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 $q$ value to evaluate $R(q)$.
a The set of all adjustable parameters in one array.
p The Parratt struct containing all the other information necessary for calculation of $R(q)$.
Returns:
The reflectivity $R(q)$.
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 $n$ point average,

\begin{eqnarray*} R(q) &=& {\sum_{a=-n}^{n} w_a R(q[1+\frac{a}{n}\frac{\delta q}{q}])}/ {\sum_{a=-n}^{n} w_a} \\ \mathrm{or} &=& {\sum_{a=-n}^{n} w_a R(q+\frac{a}{n}\delta q)}/ {\sum_{a=-n}^{n} w_a} \\ w_a &=& \exp\left(-2\left(\frac{a}{n}\right)^2\right) \end{eqnarray*}

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, $\delta q/q=\mathrm{const}$. Otherwise, the resolution can be specified as $\delta q=\mathrm{const}$. 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 $\alpha$ is the film surface coverage, then the reflectivity is calculated by

\[ R(q) = (1-\alpha )R_\mathrm{substrate}(q) + \alpha R_\mathrm{film}(q) \]

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

Finally, the background is added. 

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 $q$ value to evaluate $R(q)$.
a The set of all adjustable parameters in one array.
p The Parratt struct containing all the other information necessary for calculation of $R(q)$.
Returns:
The reflectivity $R(q)$.
For a set of $N$ layers, that run from the front layer that the neutrons enter, ($i=0$), to the back, bulk layer, ($i=N+1$), reflectivity can be calculated by

\begin{eqnarray*} k_0 &=& \sqrt{q^2/4+4\pi\rho_0} \\ k_i &=& \sqrt{k_0^2-4\pi\rho_i} \\ r_{i,i+1} &=& \frac{k_i - k_{i+1}} {k_i - k_{i+1}} e^{-2\sigma_{i+1}^2k_ik_{i+1}} \\ R_{N+1} &=& 0 \\ R_{N,N+1} &=& r_{N,N+1} \\ R_{i,i+1} &=& \frac{r_{i,i+1}+R_{i+1,i+2}e^{2id_{i+1}k_{i+1}}} {1+r_{i,i+1}R_{i+1,i+2}e^{2id_{i+1}k_{i+1}}} \\ R(q) &=& |R_{0,1}|^2 \end{eqnarray*}

where the wave vector in layer $i$ is $k_i$, and the Fresnel reflectivity from the interface between layers $i,i+1$ is $r_{i,i+1}$. 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, $q$ and $k_0$ are real valued, and everything else, including $\rho$, 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 $i=0\rightarrow n-1$. 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);
00175         tmp2 = gsl_complex_add(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) */
00191         numtr = gsl_complex_add(r, argm);
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);
00271     tmp2 = gsl_complex_add(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) */
00287     numtr = gsl_complex_add(r, argm);
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  doxygen 1.4.7