Functions | |
double | calcProfile (double z, double *a, Parratt *par) |
Calculate the profile. |
double calcProfile | ( | double | z, | |
double * | a, | |||
Parratt * | par | |||
) |
z | The value of 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 . |
where is the thickness, is the real part of the scattering length density, the roughness, each of layer . The location of the error function, , is calculated by a running total of the thicknesses.
In actuality, the total number of layers, which is given by num
, includes the substrate. For practical purposes of C arrays, the layer numbering runs as num-1
. Therefore the parameters (6 for each layer), can be found in the following manner:
for each layer i = 0 to num-1 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
When the roughness is zero, the error function becomes a Heavyside step function.
00023 { 00024 int j, indx, indx1; 00025 double thk, tmp, tmp1, tmp2; 00026 00027 /******************** 00028 * The first layer is the layer 00029 * the neutrons enter from. 00030 */ 00031 indx = par->l.idx[0]; 00032 if (par->g.plz == POLARIZED_UP) tmp = a[indx+2] + a[indx+4]; 00033 else tmp = a[indx+2] - a[indx+4]; 00034 00035 /******************** 00036 * Loop over the layers 00037 */ 00038 thk = 0.0; 00039 for (j=1; j<par->l.num; j++) 00040 { 00041 indx = par->l.idx[j]; 00042 indx1 = par->l.idx[j-1]; 00043 00044 if (par->g.plz == POLARIZED_UP) 00045 { 00046 tmp2 = 0.5*((a[indx+2] + a[indx+4]) - (a[indx1+2] + a[indx1+4])); 00047 } 00048 else 00049 { 00050 tmp2 = 0.5*((a[indx+2] - a[indx+4]) - (a[indx1+2] - a[indx1+4])); 00051 } 00052 00053 /* zero roughness : erf becomes step function */ 00054 if (a[indx+1] == 0) 00055 { 00056 if (z>=thk) tmp += 2.0 * tmp2; 00057 } 00058 else 00059 { 00060 tmp1 = (z-thk)/(M_SQRT2*fabs(a[indx+1])); 00061 tmp += tmp2 * (1.0 + gsl_sf_erf(tmp1)); 00062 } 00063 thk += a[indx+0]; 00064 } 00065 00066 return tmp; 00067 }