// The original code, of which work was not DANSE funded, // was provided by J. Cho. // And modified to fit sansmodels/sansview: JC #include #include "librefl.h" #include #include #if defined(_MSC_VER) #define NEED_ERF #endif #if defined(NEED_ERF) /* erf.c - public domain implementation of error function erf(3m) reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten (New Algorithm handbook in C language) (Gijyutsu hyouron sha, Tokyo, 1991) p.227 [in Japanese] */ #ifdef _WIN32 # include # if !defined __MINGW32__ || defined __NO_ISOCEXT # ifndef isnan # define isnan(x) _isnan(x) # endif # ifndef isinf # define isinf(x) (!_finite(x) && !_isnan(x)) # endif # ifndef finite # define finite(x) _finite(x) # endif # endif #endif static double q_gamma(double, double, double); /* Incomplete gamma function 1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt */ static double p_gamma(double a, double x, double loggamma_a) { int k; double result, term, previous; if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a); if (x == 0) return 0; result = term = exp(a * log(x) - x - loggamma_a) / a; for (k = 1; k < 1000; k++) { term *= x / (a + k); previous = result; result += term; if (result == previous) return result; } fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__); return result; } /* Incomplete gamma function 1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt */ static double q_gamma(double a, double x, double loggamma_a) { int k; double result, w, temp, previous; double la = 1, lb = 1 + x - a; /* Laguerre polynomial */ if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a); w = exp(a * log(x) - x - loggamma_a); result = w / lb; for (k = 2; k < 1000; k++) { temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k; la = lb; lb = temp; w *= (k - 1 - a) / k; temp = w / (la * lb); previous = result; result += temp; if (result == previous) return result; } fprintf(stderr, "erf.c:%d:q_gamma() could not converge.", __LINE__); return result; } #define LOG_PI_OVER_2 0.572364942924700087071713675675 /* log_e(PI)/2 */ double erf(double x) { if (!finite(x)) { if (isnan(x)) return x; /* erf(NaN) = NaN */ return (x>0 ? 1.0 : -1.0); /* erf(+-inf) = +-1.0 */ } if (x >= 0) return p_gamma(0.5, x * x, LOG_PI_OVER_2); else return - p_gamma(0.5, x * x, LOG_PI_OVER_2); } double erfc(double x) { if (!finite(x)) { if (isnan(x)) return x; /* erfc(NaN) = NaN */ return (x>0 ? 0.0 : 2.0); /* erfc(+-inf) = 0.0, 2.0 */ } if (x >= 0) return q_gamma(0.5, x * x, LOG_PI_OVER_2); else return 1 + p_gamma(0.5, x * x, LOG_PI_OVER_2); } #endif // NEED_ERF complex cassign(real, imag) double real, imag; { complex x; x.re = real; x.im = imag; return x; } complex cplx_add(x,y) complex x,y; { complex z; z.re = x.re + y.re; z.im = x.im + y.im; return z; } complex rcmult(x,y) double x; complex y; { complex z; z.re = x*y.re; z.im = x*y.im; return z; } complex cplx_sub(x,y) complex x,y; { complex z; z.re = x.re - y.re; z.im = x.im - y.im; return z; } complex cplx_mult(x,y) complex x,y; { complex z; z.re = x.re*y.re - x.im*y.im; z.im = x.re*y.im + x.im*y.re; return z; } complex cplx_div(x,y) complex x,y; { complex z; z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im); z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im); return z; } complex cplx_exp(b) complex b; { complex z; double br,bi; br=b.re; bi=b.im; z.re = exp(br)*cos(bi); z.im = exp(br)*sin(bi); return z; } complex cplx_sqrt(z) //see Schaum`s Math Handbook p. 22, 6.6 and 6.10 complex z; { complex c; double zr,zi,x,y,r,w; zr=z.re; zi=z.im; if (zr==0.0 && zi==0.0) { c.re=0.0; c.im=0.0; return c; } else { x=fabs(zr); y=fabs(zi); if (x>y) { r=y/x; w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r))); } else { r=x/y; w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r))); } if (zr >=0.0) { c.re=w; c.im=zi/(2.0*w); } else { c.im=(zi >= 0) ? w : -w; c.re=zi/(2.0*c.im); } return c; } } complex cplx_cos(b) complex b; { complex zero,two,z,i,bi,negbi; zero = cassign(0.0,0.0); two = cassign(2.0,0.0); i = cassign(0.0,1.0); bi = cplx_mult(b,i); negbi = cplx_sub(zero,bi); z = cplx_div(cplx_add(cplx_exp(bi),cplx_exp(negbi)),two); return z; } // normalized and modified erf // | // 1 + __ - - - - // | _ // | _ // | __ // 0 + - - - // |-------------+------------+-- // 0 center n_sub ---> // ind // // n_sub = total no. of bins(or sublayers) // ind = x position: 0 to max // nu = max x to integration double err_mod_func(double n_sub, double ind, double nu) { double center, func; if (nu == 0.0) nu = 1e-14; if (n_sub == 0.0) n_sub = 1.0; //ind = (n_sub-1.0)/2.0-1.0 +ind; center = n_sub/2.0; // transform it so that min(ind) = 0 ind -= center; // normalize by max limit ind /= center; // divide by sqrt(2) to get Gaussian func nu /= sqrt(2.0); ind *= nu; // re-scale and normalize it so that max(erf)=1, min(erf)=0 func = erf(ind)/erf(nu)/2.0; // shift it by +0.5 in y-direction so that min(erf) = 0 func += 0.5; return func; } double linearfunc(double n_sub, double ind, double nu) { double bin_size, func; if (n_sub == 0.0) n_sub = 1.0; bin_size = 1.0/n_sub; //size of each sub-layer // rescale ind *= bin_size; func = ind; return func; } // use the right hand side from the center of power func double power_r(double n_sub, double ind, double nu) { double bin_size,func; if (nu == 0.0) nu = 1e-14; if (n_sub == 0.0) n_sub = 1.0; bin_size = 1.0/n_sub; //size of each sub-layer // rescale ind *= bin_size; func = pow(ind, nu); return func; } // use the left hand side from the center of power func double power_l(double n_sub, double ind, double nu) { double bin_size, func; if (nu == 0.0) nu = 1e-14; if (n_sub == 0.0) n_sub = 1.0; bin_size = 1.0/n_sub; //size of each sub-layer // rescale ind *= bin_size; func = 1.0-pow((1.0-ind),nu); return func; } // use 1-exp func from x=0 to x=1 double exp_r(double n_sub, double ind, double nu) { double bin_size, func; if (nu == 0.0) nu = 1e-14; if (n_sub == 0.0) n_sub = 1.0; bin_size = 1.0/n_sub; //size of each sub-layer // rescale ind *= bin_size; // modify func so that func(0) =0 and func(max)=1 func = 1.0-exp(-nu*ind); // normalize by its max func /= (1.0-exp(-nu)); return func; } // use the left hand side mirror image of exp func double exp_l(double n_sub, double ind, double nu) { double bin_size, func; if (nu == 0.0) nu = 1e-14; if (n_sub == 0.0) n_sub = 1.0; bin_size = 1.0/n_sub; //size of each sub-layer // rescale ind *= bin_size; // modify func func = exp(-nu*(1.0-ind))-exp(-nu); // normalize by its max func /= (1.0-exp(-nu)); return func; } // To select function called // At nu = 0 (singular point), call line function double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r) { double sld_i, func; // this condition protects an error from the singular point if (nu == 0.0){ nu = 1e-13; } // select func switch(fun_type){ case 1 : func = power_r(n_sub, i, nu); break; case 2 : func = power_l(n_sub, i, nu); break; case 3 : func = exp_r(n_sub, i, nu); break; case 4 : func = exp_l(n_sub, i, nu); break; case 5 : func = linearfunc(n_sub, i, nu); break; default: func = err_mod_func(n_sub, i, nu); break; } // compute sld if (sld_r>sld_l){ sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick); } else if (sld_rsld_l){ sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick); } else if (sld_r