# Calculating the reflectivity

## Functions

double calcReflectivity (double q, double *a, Parratt *par)
Wrapper for calculating the reflectivity with background, instrumental resolution, and surface coverage.
double calcParratt (double q, double *a, Parratt *par)
Actually calculates the reflectivity.
double calcNoFilm (double q, double *a, Parratt *par)
Calculates the reflectivity of the bulk/substrate interface.
double lookupResolution (double q, Parratt *par)
Returns the instrumental resolution at a particular q value.

## Function Documentation

 double calcReflectivity ( double q, double * a, Parratt * par )

Parameters:
 q The value to evaluate . a The set of all adjustable parameters in one array. p The Parratt struct containing all the other information necessary for calculation of .
Returns:
The reflectivity .
This function is a wrapper to add the background, instrument resolution and surface coverage. The actual reflectivity is calculated by calls to the function calcParratt() .

The instrumental resolution is determined by calculating the reflectivity averaged over several points with Gaussian weights. It's a method that was borrowed and made more general from Parratt32. It works by the following, assuming a point average,

In the first equation, we can specify the instrumental resolution present in the data as in the same manner as it is typically recorded a Chalk River, . Otherwise, the resolution can be specified as . The Gaussian averaging of points has been pre-calculated for 9 point average. The resolution is held in the struct reservd for unchanging global parameterspar->g. If no resolution smearing is required, the method is skipped for speed.

If the surface coverage of the film beign studied is other than 1, than a linear approximation of the coverage can be included in the fit. First, the reflectivity of the bulk/substrate interface is calculated in the same manner as the total film, with the same instrument resolution. If is the film surface coverage, then the reflectivity is calculated by

There is no check that the value of lies between 0 and 1.

00024 {
00025     double tmp, tmp1;
00026
00027     tmp = calcParratt(q, a, par);
00028
00029     /* determine resolution for this q */
00030     double resolution = par->g.rsn;
00031     if (par->g.rl == ARRAY)
00032         resolution = lookupResolution(q, par);
00033
00034     /* instrument resolution, skip if zero for speed */
00035     if (resolution != 0.0)
00036     {
00037
00038         if (par->g.rtp == CONSTANT)
00039         {
00040             tmp += 0.882*calcParratt(q+(0.25*resolution), a, par);
00041             tmp += 0.882*calcParratt(q-(0.25*resolution), a, par);
00042             tmp += 0.607*calcParratt(q+(0.50*resolution), a, par);
00043             tmp += 0.607*calcParratt(q-(0.50*resolution), a, par);
00044             tmp += 0.325*calcParratt(q+(0.75*resolution), a, par);
00045             tmp += 0.325*calcParratt(q-(0.75*resolution), a, par);
00046             tmp += 0.135*calcParratt(q+(1.00*resolution), a, par);
00047             tmp += 0.135*calcParratt(q-(1.00*resolution), a, par);
00048         }
00049         else if (par->g.rtp == RELATIVE)
00050         {
00051             tmp += 0.882*calcParratt(q*(1.0+(0.25*resolution)), a, par);
00052             tmp += 0.882*calcParratt(q*(1.0-(0.25*resolution)), a, par);
00053             tmp += 0.607*calcParratt(q*(1.0+(0.50*resolution)), a, par);
00054             tmp += 0.607*calcParratt(q*(1.0-(0.50*resolution)), a, par);
00055             tmp += 0.325*calcParratt(q*(1.0+(0.75*resolution)), a, par);
00056             tmp += 0.325*calcParratt(q*(1.0-(0.75*resolution)), a, par);
00057             tmp += 0.135*calcParratt(q*(1.0+(1.00*resolution)), a, par);
00058             tmp += 0.135*calcParratt(q*(1.0-(1.00*resolution)), a, par);
00059         }
00060         tmp /= (4.898);
00061     }
00062
00063     /* coverage */
00064     if (a[1] != 1.0)
00065     {
00066         tmp1 = calcNoFilm(q, a, par);
00067
00068         /* instrument resolution, skip if zero for speed */
00069         if (resolution != 0.0)
00070         {
00071             if (par->g.rtp == CONSTANT)
00072             {
00073                 tmp1 += 0.882*calcNoFilm(q+(0.25*resolution), a, par);
00074                 tmp1 += 0.882*calcNoFilm(q-(0.25*resolution), a, par);
00075                 tmp1 += 0.607*calcNoFilm(q+(0.50*resolution), a, par);
00076                 tmp1 += 0.607*calcNoFilm(q-(0.50*resolution), a, par);
00077                 tmp1 += 0.325*calcNoFilm(q+(0.75*resolution), a, par);
00078                 tmp1 += 0.325*calcNoFilm(q-(0.75*resolution), a, par);
00079                 tmp1 += 0.135*calcNoFilm(q+(1.00*resolution), a, par);
00080                 tmp1 += 0.135*calcNoFilm(q-(1.00*resolution), a, par);
00081             }
00082             else if (par->g.rtp == RELATIVE)
00083             {
00084                 tmp1 += 0.882*calcNoFilm(q*(1.0+(0.25*resolution)), a, par);
00085                 tmp1 += 0.882*calcNoFilm(q*(1.0-(0.25*resolution)), a, par);
00086                 tmp1 += 0.607*calcNoFilm(q*(1.0+(0.50*resolution)), a, par);
00087                 tmp1 += 0.607*calcNoFilm(q*(1.0-(0.50*resolution)), a, par);
00088                 tmp1 += 0.325*calcNoFilm(q*(1.0+(0.75*resolution)), a, par);
00089                 tmp1 += 0.325*calcNoFilm(q*(1.0-(0.75*resolution)), a, par);
00090                 tmp1 += 0.135*calcNoFilm(q*(1.0+(1.00*resolution)), a, par);
00091                 tmp1 += 0.135*calcNoFilm(q*(1.0-(1.00*resolution)), a, par);
00092             }
00093             tmp1 /= (4.898);
00094         }
00095
00096         /* coverage */
00097         tmp = (1.0-a[1])*tmp1 + a[1]*tmp;
00098     }
00099
00100     /* background */
00101     tmp += (a[0] * 1E-6);
00102
00103     return tmp;
00104 }


 double calcParratt ( double q, double * a, Parratt * par )

Parameters:
 q The value to evaluate . a The set of all adjustable parameters in one array. p The Parratt struct containing all the other information necessary for calculation of .
Returns:
The reflectivity .
For a set of layers, that run from the front layer that the neutrons enter, (), to the back, bulk layer, (), reflectivity can be calculated by

where the wave vector in layer is , and the Fresnel reflectivity from the interface between layers is . The reflectivity of the entire system is given by the reflectivity from the top surface, and is shown in the last line above.

In this way the reflectivity is calculated from the back surface to the front. Here, and are real valued, and everything else, including , is complex.

As with the profile calculation, the total number of layers, which is given by the passed parameter n, includes the substrate. For practical purposes of C arrays, the layer numbering runs as . Therefore the parameters (6 for each layer), can be found in the following manner:

for each layer i = n-2 to 0
index = idx[i]
par[index+0] = Thickness
par[index+1] = Roughness
par[index+2] = Real part of nuclear SLD
par[index+3] = Imaginary part of nuclear SLD
par[index+4] = Real part of magnetic SLD
par[index+5] = Imaginary part of magnetic SLD


The main loop runs from layer n-2 to 0. The function relies on the GNU scientific library complex math functions.

00107 {
00108     int          l;
00109     unsigned int ix1, ix2;
00110     double       k0z, tmpr, tmpi, rsld, isld, result;
00111     gsl_complex  kzn, kznp1;
00112     gsl_complex  fr, r, R;
00113     gsl_complex  tmp1, tmp2, tmp3;
00114     gsl_complex  argm, numtr, denom;
00115
00116     /* k0z */
00117     ix1 = par->l.idx[0];
00118     if (par->g.plz == POLARIZED_UP)
00119     {
00120         rsld = (a[ix1+2] + a[ix1+4]);
00121         isld = (a[ix1+3] + a[ix1+5]);
00122     }
00123     else
00124     {
00125         rsld = (a[ix1+2] - a[ix1+4]);
00126         isld = (a[ix1+3] - a[ix1+5]);
00127     }
00128     tmpr = (q * q / 4.0) + (4.0 * M_PI * rsld * 1E-6);
00129     tmpi = 4.0 * M_PI * isld * 1E-6;
00130     GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
00131     k0z  = gsl_complex_abs(gsl_complex_sqrt(tmp1));
00132
00133     GSL_SET_COMPLEX(&R, 0.0, 0.0);
00134     GSL_SET_COMPLEX(&r, 0.0, 0.0);
00135
00136     for (l = par->l.num-2; l>=0; l--)
00137     {
00138         ix1 = par->l.idx[l];
00139         ix2 = par->l.idx[l+1];
00140
00141         /* k_n */
00142         if (par->g.plz == POLARIZED_UP)
00143         {
00144             rsld = (a[ix1+2] + a[ix1+4]);
00145             isld = (a[ix1+3] + a[ix1+5]);
00146         }
00147         else
00148         {
00149             rsld = (a[ix1+2] - a[ix1+4]);
00150             isld = (a[ix1+3] - a[ix1+5]);
00151         }
00152         tmpr = (k0z * k0z)-(4.0 * M_PI * rsld * 1E-6);
00153         tmpi = 4.0 * M_PI * isld * 1E-6;
00154         GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
00155         kzn   = gsl_complex_sqrt(tmp1);
00156
00157         /* k_n+1 */
00158         if (par->g.plz == POLARIZED_UP)
00159         {
00160             rsld = (a[ix2+2] + a[ix2+4]);
00161             isld = (a[ix2+3] + a[ix2+5]);
00162         }
00163         else
00164         {
00165             rsld = (a[ix2+2] - a[ix2+4]);
00166             isld = (a[ix2+3] - a[ix2+5]);
00167         }
00168         tmpr = (k0z * k0z)-(4.0 * M_PI * rsld * 1E-6);
00169         tmpi = 4.0 * M_PI * isld * 1E-6;
00170         GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
00171         kznp1 = gsl_complex_sqrt(tmp1);
00172
00173         /* Fresnel */
00174         tmp1 = gsl_complex_sub(kzn, kznp1);
00176         fr   = gsl_complex_div(tmp1, tmp2);
00177
00178         /* r = Fresnel * roughness */
00179         tmp1 = gsl_complex_mul(kzn,kznp1);
00180         tmpr = -2.0*a[ix2+1]*a[ix2+1];
00181         tmp2 = gsl_complex_mul_real(tmp1, tmpr);
00182         tmp3 = gsl_complex_exp(tmp2);
00183         r    = gsl_complex_mul(fr, tmp3);
00184
00185         /* argm = X_n+1 * exp(2 i k_n+1 d_n+1) */
00186         tmp1 = gsl_complex_mul_imag(kznp1, 2.0*a[ix2+0]);
00187         tmp2 = gsl_complex_exp(tmp1);
00188         argm = gsl_complex_mul(R, tmp2);
00189
00190         /* numtr = r_n,n+1 + X_n+1 * exp(2 i d_n+1) */
00192
00193         /* denom = 1 + r_n,n+1 * X_n+1 * exp(2 i k_n+1 d_n+1) */
00194         denom = gsl_complex_add_real(gsl_complex_mul(r, argm), 1.0);
00195
00196         R = gsl_complex_div(numtr, denom);
00197     }
00198
00199     result = gsl_complex_abs2(R);
00200
00201     return result;
00202 }


 double calcNoFilm ( double q, double * a, Parratt * par )

This function calculates the reflectivity assuming no film on the substrate. Equation are the same as those found in calcParratt().

00205 {
00206     unsigned int ix1, ix2;
00207     double       k0z, tmpr, tmpi, rsld, isld, result;
00208     gsl_complex  kzn, kznp1;
00209     gsl_complex  fr, r, R;
00210     gsl_complex  tmp1, tmp2, tmp3;
00211     gsl_complex  argm, numtr, denom;
00212
00213     /* k0z */
00214     ix1 = par->l.idx[0];
00215     if (par->g.plz == POLARIZED_UP)
00216     {
00217         rsld = (a[ix1+2] + a[ix1+4]);
00218         isld = (a[ix1+3] + a[ix1+5]);
00219     }
00220     else
00221     {
00222         rsld = (a[ix1+2] - a[ix1+4]);
00223         isld = (a[ix1+3] - a[ix1+5]);
00224     }
00225     tmpr = (q * q / 4.0) + (4.0 * M_PI * rsld * 1E-6);
00226     tmpi = 4.0 * M_PI * isld * 1E-6;
00227     GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
00228     k0z  = gsl_complex_abs(gsl_complex_sqrt(tmp1));
00229
00230     GSL_SET_COMPLEX(&R, 0.0, 0.0);
00231     GSL_SET_COMPLEX(&r, 0.0, 0.0);
00232
00233     /* Index first and last layers only.*/
00234     ix1 = par->l.idx[0];
00235     ix2 = par->l.idx[par->l.num-1];
00236
00237     /* k_n */
00238     if (par->g.plz == POLARIZED_UP)
00239     {
00240         rsld = (a[ix1+2] + a[ix1+4]);
00241         isld = (a[ix1+3] + a[ix1+5]);
00242     }
00243     else
00244     {
00245         rsld = (a[ix1+2] - a[ix1+4]);
00246         isld = (a[ix1+3] - a[ix1+5]);
00247     }
00248     tmpr = (k0z * k0z)-(4.0 * M_PI * rsld * 1E-6);
00249     tmpi = 4.0 * M_PI * isld * 1E-6;
00250     GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
00251     kzn   = gsl_complex_sqrt(tmp1);
00252
00253     /* k_n+1 */
00254     if (par->g.plz == POLARIZED_UP)
00255     {
00256         rsld = (a[ix2+2] + a[ix2+4]);
00257         isld = (a[ix2+3] + a[ix2+5]);
00258     }
00259     else
00260     {
00261         rsld = (a[ix2+2] - a[ix2+4]);
00262         isld = (a[ix2+3] - a[ix2+5]);
00263     }
00264     tmpr = (k0z * k0z)-(4.0 * M_PI * rsld * 1E-6);
00265     tmpi = 4.0 * M_PI * isld * 1E-6;
00266     GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
00267     kznp1 = gsl_complex_sqrt(tmp1);
00268
00269     /* Fresnel */
00270     tmp1 = gsl_complex_sub(kzn, kznp1);
00272     fr   = gsl_complex_div(tmp1, tmp2);
00273
00274     /* r = Fresnel * roughness */
00275     tmp1 = gsl_complex_mul(kzn,kznp1);
00276     tmpr = -2.0*a[ix2+1]*a[ix2+1];
00277     tmp2 = gsl_complex_mul_real(tmp1, tmpr);
00278     tmp3 = gsl_complex_exp(tmp2);
00279     r    = gsl_complex_mul(fr, tmp3);
00280
00281     /* argm = X_n+1 * exp(2 i k_n+1 d_n+1) */
00282     tmp1 = gsl_complex_mul_imag(kznp1, 2.0*a[ix2+0]);
00283     tmp2 = gsl_complex_exp(tmp1);
00284     argm = gsl_complex_mul(R, tmp2);
00285
00286     /* numtr = r_n,n+1 + X_n+1 * exp(2 i d_n+1) */
00288
00289     /* denom = 1 + r_n,n+1 * X_n+1 * exp(2 i k_n+1 d_n+1) */
00290     denom = gsl_complex_add_real(gsl_complex_mul(r, argm), 1.0);
00291
00292     R = gsl_complex_div(numtr, denom);
00293
00294     result = gsl_complex_abs2(R);
00295
00296     return result;
00297 }


 double lookupResolution ( double q, Parratt * par )

This function returns the value of the instrumental resolution at a particular q value, based on the user-supplied list of resolution values. If the exact q-value does not exist in the user-supplied listing, then the nearest smaller q is assumed.

00300 {
00301
00302     int i=1;
00303     /* Walk through the array of q-values until you walk past
00304        the target value */
00305     /* The first part of the while check prevents us from
00306        ever accessing past the end of the array */
00307     while( (i<(par->g.rsNum)) && (par->g.rsArrayq[i]<=q) )
00308     {
00309         i++;
00310     }
00311     /* Return the previous element (since we just stepped
00312        past the target */
00313     return par->g.rsArrayval[i-1];
00314
00315 }


Generated on Wed Mar 14 13:24:56 2007 for Parratt by  1.4.7