Changeset 2d1b700 in sasview for sansmodels/src/c_models/refl.cpp


Ignore:
Timestamp:
Jan 5, 2012 8:52:13 AM (12 years ago)
Author:
Mathieu Doucet <doucetm@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
dd60b45
Parents:
7ffa8196
Message:

refactor refl model and auto-generate c++ wrapper at compile time.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/c_models/refl.cpp

    r67424cd r2d1b700  
    11 
    22#include <math.h> 
    3 #include "models.hh" 
    43#include "parameters.hh" 
    54#include <stdio.h> 
     5#include <stdlib.h> 
     6#include "refl.h" 
    67using namespace std; 
    78 
    89extern "C" { 
    9         #include "refl.h" 
    10 } 
    11  
     10#include "libmultifunc/librefl.h" 
     11} 
     12 
     13#define lamda 4.62 
     14 
     15double re_kernel(double dp[], double q) { 
     16  int n = dp[0]; 
     17  int i,j; 
     18 
     19  double scale = dp[1]; 
     20  double thick_inter_sub = dp[2]; 
     21  double sld_sub = dp[4]; 
     22  double sld_super = dp[5]; 
     23  double background = dp[6]; 
     24 
     25  double total_thick=0.0; 
     26  double nsl=21.0; //nsl = Num_sub_layer: 
     27  int n_s; 
     28  double sld_i,dz,phi,R,ko2; 
     29  double fun; 
     30  double pi; 
     31 
     32  double* sld; 
     33  double* thick_inter; 
     34  double* thick; 
     35  int*fun_type; 
     36  complex  phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,two,n_sub,n_sup,knp1,Xnp1; 
     37 
     38  sld = (double*)malloc((n+2)*sizeof(double)); 
     39  thick_inter = (double*)malloc((n+2)*sizeof(double)); 
     40  thick = (double*)malloc((n+2)*sizeof(double)); 
     41  fun_type = (int*)malloc((n+2)*sizeof(int)); 
     42 
     43  fun_type[0] = dp[3]; 
     44  for (i =1; i<=n; i++){ 
     45    sld[i] = dp[i+6]; 
     46    thick_inter[i]= dp[i+16]; 
     47    thick[i] = dp[i+26]; 
     48    fun_type[i] = dp[i+36]; 
     49    total_thick += thick[i] + thick_inter[i]; 
     50  } 
     51  sld[0] = sld_sub; 
     52  sld[n+1] = sld_super; 
     53 
     54  thick[0] = total_thick/5.0; 
     55  thick[n+1] = total_thick/5.0; 
     56  thick_inter[0] = thick_inter_sub; 
     57  thick_inter[n+1] = 0.0; 
     58 
     59  pi = 4.0*atan(1.0); 
     60  Xn = cassign(0.0,0.0); 
     61  one = cassign(1.0,0.0); 
     62  two = cassign(0.0,-2.0); 
     63 
     64  //Checking if floor is available. 
     65  //no imaginary sld inputs in this function yet 
     66  n_sub=cassign(1.0-sld_sub*pow(lamda,2.0)/(2.0*pi),0.0); 
     67  n_sup=cassign(1.0-sld_super*pow(lamda,2.0)/(2.0*pi),0.0); 
     68  ko2 = pow(2.0*pi/lamda,2.0); 
     69 
     70  phi = asin(lamda*q/(4.0*pi)); 
     71  phi1 = cplx_div(rcmult(phi,one),n_sup); 
     72  alpha = cplx_mult(n_sup,cplx_cos(phi1)); 
     73  alpha2 = cplx_mult(alpha,alpha); 
     74 
     75  nnp1=n_sub; 
     76  knp1=cplx_sqrt(rcmult(ko2,cplx_sub(cplx_mult(nnp1,nnp1),alpha2)));  //nnp1*ko*sin(phinp1) 
     77  Xnp1=cassign(0.0,0.0); 
     78  dz = 0.0; 
     79  // iteration for # of layers +sub from the top 
     80  for (i=1;i<=n+1; i++){ 
     81    if (fun_type[i-1]==1) 
     82      fun = 5; 
     83    else 
     84      fun = 0; 
     85    //iteration for 9 sub-layers 
     86    for (j=0;j<2;j++){ 
     87      for (n_s=0;n_s<nsl; n_s++){ 
     88        if (j==1){ 
     89          if (i==n+1) 
     90            break; 
     91          dz = thick[i]; 
     92          sld_i = sld[i]; 
     93        } 
     94        else{ 
     95          dz = thick_inter[i-1]/nsl;//nsl; 
     96          if (sld[i-1] == sld[i]){ 
     97            sld_i = sld[i]; 
     98          } 
     99          else{ 
     100            sld_i = intersldfunc(fun,nsl, n_s+0.5, 2.5, sld[i-1], sld[i]); 
     101          } 
     102        } 
     103        nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi),0.0); 
     104        nn2=cplx_mult(nn,nn); 
     105 
     106        kn=cplx_sqrt(rcmult(ko2,cplx_sub(nn2,alpha2)));        //nn*ko*sin(phin) 
     107        an=cplx_exp(rcmult(dz,cplx_mult(two,kn))); 
     108 
     109        fnm=cplx_sub(kn,knp1); 
     110        fnp=cplx_add(kn,knp1); 
     111        rn=cplx_div(fnm,fnp); 
     112        Xn=cplx_mult(an,cplx_div(cplx_add(rn,Xnp1),cplx_add(one,cplx_mult(rn,Xnp1))));    //Xn=an*((rn+Xnp1*anp1)/(1+rn*Xnp1*anp1)) 
     113 
     114        Xnp1=Xn; 
     115        knp1=kn; 
     116 
     117        if (j==1) 
     118          break; 
     119      } 
     120    } 
     121  } 
     122  R=pow(Xn.re,2.0)+pow(Xn.im,2.0); 
     123  // This temperarily fixes the total reflection for Rfunction and linear. 
     124  // ToDo: Show why it happens that Xn.re=0 and Xn.im >1! 
     125  if (Xn.im == 0.0){ 
     126    R=1.0; 
     127  } 
     128  R *= scale; 
     129  R += background; 
     130 
     131  free(sld); 
     132  free(thick_inter); 
     133  free(thick); 
     134  free(fun_type); 
     135 
     136  return R; 
     137 
     138} 
    12139ReflModel :: ReflModel() { 
    13         n_layers = Parameter(1.0); 
    14         scale = Parameter(1.0); 
    15         thick_inter0 = Parameter(1.0); 
    16         func_inter0 = Parameter(0); 
    17         sld_bottom0 = Parameter(2.07e-06); 
    18         sld_medium = Parameter(1.0e-06); 
    19     background = Parameter(0.0); 
    20  
    21  
    22     sld_flat1 = Parameter(3.0e-06); 
    23     sld_flat2 = Parameter(3.5e-06); 
    24     sld_flat3 = Parameter(4.0e-06); 
    25     sld_flat4 = Parameter(3.5e-06); 
    26     sld_flat5 = Parameter(4.0e-06); 
    27     sld_flat6 = Parameter(3.5e-06); 
    28     sld_flat7 = Parameter(4.0e-06); 
    29     sld_flat8 = Parameter(3.5e-06); 
    30     sld_flat9 = Parameter(4.0e-06); 
    31     sld_flat10 = Parameter(3.5e-06); 
    32  
    33  
    34     thick_inter1 = Parameter(1); 
    35     thick_inter2 = Parameter(1); 
    36     thick_inter3 = Parameter(1); 
    37     thick_inter4 = Parameter(1); 
    38     thick_inter5 = Parameter(1); 
    39     thick_inter6 = Parameter(1); 
    40     thick_inter7 = Parameter(1); 
    41     thick_inter8 = Parameter(1); 
    42     thick_inter9 = Parameter(1); 
    43     thick_inter10 = Parameter(1); 
    44  
    45  
    46     thick_flat1 = Parameter(15); 
    47     thick_flat2 = Parameter(100); 
    48     thick_flat3 = Parameter(100); 
    49     thick_flat4 = Parameter(100); 
    50     thick_flat5 = Parameter(100); 
    51     thick_flat6 = Parameter(100); 
    52     thick_flat7 = Parameter(100); 
    53     thick_flat8 = Parameter(100); 
    54     thick_flat9 = Parameter(100); 
    55     thick_flat10 = Parameter(100); 
    56  
    57  
    58     func_inter1 = Parameter(0); 
    59     func_inter2 = Parameter(0); 
    60     func_inter3 = Parameter(0); 
    61     func_inter4 = Parameter(0); 
    62     func_inter5 = Parameter(0); 
    63     func_inter6 = Parameter(0); 
    64     func_inter7 = Parameter(0); 
    65     func_inter8 = Parameter(0); 
    66     func_inter9 = Parameter(0); 
    67     func_inter10 = Parameter(0); 
     140  n_layers = Parameter(1.0); 
     141  scale = Parameter(1.0); 
     142  thick_inter0 = Parameter(1.0); 
     143  func_inter0 = Parameter(0); 
     144  sld_bottom0 = Parameter(2.07e-06); 
     145  sld_medium = Parameter(1.0e-06); 
     146  background = Parameter(0.0); 
     147 
     148 
     149  sld_flat1 = Parameter(3.0e-06); 
     150  sld_flat2 = Parameter(3.5e-06); 
     151  sld_flat3 = Parameter(4.0e-06); 
     152  sld_flat4 = Parameter(3.5e-06); 
     153  sld_flat5 = Parameter(4.0e-06); 
     154  sld_flat6 = Parameter(3.5e-06); 
     155  sld_flat7 = Parameter(4.0e-06); 
     156  sld_flat8 = Parameter(3.5e-06); 
     157  sld_flat9 = Parameter(4.0e-06); 
     158  sld_flat10 = Parameter(3.5e-06); 
     159 
     160 
     161  thick_inter1 = Parameter(1); 
     162  thick_inter2 = Parameter(1); 
     163  thick_inter3 = Parameter(1); 
     164  thick_inter4 = Parameter(1); 
     165  thick_inter5 = Parameter(1); 
     166  thick_inter6 = Parameter(1); 
     167  thick_inter7 = Parameter(1); 
     168  thick_inter8 = Parameter(1); 
     169  thick_inter9 = Parameter(1); 
     170  thick_inter10 = Parameter(1); 
     171 
     172 
     173  thick_flat1 = Parameter(15); 
     174  thick_flat2 = Parameter(100); 
     175  thick_flat3 = Parameter(100); 
     176  thick_flat4 = Parameter(100); 
     177  thick_flat5 = Parameter(100); 
     178  thick_flat6 = Parameter(100); 
     179  thick_flat7 = Parameter(100); 
     180  thick_flat8 = Parameter(100); 
     181  thick_flat9 = Parameter(100); 
     182  thick_flat10 = Parameter(100); 
     183 
     184 
     185  func_inter1 = Parameter(0); 
     186  func_inter2 = Parameter(0); 
     187  func_inter3 = Parameter(0); 
     188  func_inter4 = Parameter(0); 
     189  func_inter5 = Parameter(0); 
     190  func_inter6 = Parameter(0); 
     191  func_inter7 = Parameter(0); 
     192  func_inter8 = Parameter(0); 
     193  func_inter9 = Parameter(0); 
     194  func_inter10 = Parameter(0); 
    68195 
    69196} 
    … …  
    75202 */ 
    76203double ReflModel :: operator()(double q) { 
    77         double dp[47]; 
    78         // Fill parameter array for IGOR library 
    79         // Add the background after averaging 
    80         dp[0] = n_layers(); 
    81         dp[1] = scale(); 
    82         dp[2] = thick_inter0(); 
    83         dp[3] = func_inter0(); 
    84         dp[4] = sld_bottom0(); 
    85         dp[5] = sld_medium(); 
    86         dp[6] = background(); 
    87  
    88         dp[7] = sld_flat1(); 
    89         dp[8] = sld_flat2(); 
    90         dp[9] = sld_flat3(); 
    91         dp[10] = sld_flat4(); 
    92         dp[11] = sld_flat5(); 
    93         dp[12] = sld_flat6(); 
    94         dp[13] = sld_flat7(); 
    95         dp[14] = sld_flat8(); 
    96         dp[15] = sld_flat9(); 
    97         dp[16] = sld_flat10(); 
    98  
    99         dp[17] = thick_inter1(); 
    100         dp[18] = thick_inter2(); 
    101         dp[19] = thick_inter3(); 
    102         dp[20] = thick_inter4(); 
    103         dp[21] = thick_inter5(); 
    104         dp[22] = thick_inter6(); 
    105         dp[23] = thick_inter7(); 
    106         dp[24] = thick_inter8(); 
    107         dp[25] = thick_inter9(); 
    108         dp[26] = thick_inter10(); 
    109  
    110         dp[27] = thick_flat1(); 
    111         dp[28] = thick_flat2(); 
    112         dp[29] = thick_flat3(); 
    113         dp[30] = thick_flat4(); 
    114         dp[31] = thick_flat5(); 
    115         dp[32] = thick_flat6(); 
    116         dp[33] = thick_flat7(); 
    117         dp[34] = thick_flat8(); 
    118         dp[35] = thick_flat9(); 
    119         dp[36] = thick_flat10(); 
    120  
    121         dp[37] = func_inter1(); 
    122         dp[38] = func_inter2(); 
    123         dp[39] = func_inter3(); 
    124         dp[40] = func_inter4(); 
    125         dp[41] = func_inter5(); 
    126         dp[42] = func_inter6(); 
    127         dp[43] = func_inter7(); 
    128         dp[44] = func_inter8(); 
    129         dp[45] = func_inter9(); 
    130         dp[46] = func_inter10(); 
    131  
    132         // Get the dispersion points for the radius 
    133         //vector<WeightPoint> weights_thick; 
    134         //thick_inter0.get_weights(weights_thick_inter0); 
    135  
    136  
    137         return re_kernel(dp,q); 
     204  double dp[47]; 
     205  // Fill parameter array for IGOR library 
     206  // Add the background after averaging 
     207  dp[0] = n_layers(); 
     208  dp[1] = scale(); 
     209  dp[2] = thick_inter0(); 
     210  dp[3] = func_inter0(); 
     211  dp[4] = sld_bottom0(); 
     212  dp[5] = sld_medium(); 
     213  dp[6] = background(); 
     214 
     215  dp[7] = sld_flat1(); 
     216  dp[8] = sld_flat2(); 
     217  dp[9] = sld_flat3(); 
     218  dp[10] = sld_flat4(); 
     219  dp[11] = sld_flat5(); 
     220  dp[12] = sld_flat6(); 
     221  dp[13] = sld_flat7(); 
     222  dp[14] = sld_flat8(); 
     223  dp[15] = sld_flat9(); 
     224  dp[16] = sld_flat10(); 
     225 
     226  dp[17] = thick_inter1(); 
     227  dp[18] = thick_inter2(); 
     228  dp[19] = thick_inter3(); 
     229  dp[20] = thick_inter4(); 
     230  dp[21] = thick_inter5(); 
     231  dp[22] = thick_inter6(); 
     232  dp[23] = thick_inter7(); 
     233  dp[24] = thick_inter8(); 
     234  dp[25] = thick_inter9(); 
     235  dp[26] = thick_inter10(); 
     236 
     237  dp[27] = thick_flat1(); 
     238  dp[28] = thick_flat2(); 
     239  dp[29] = thick_flat3(); 
     240  dp[30] = thick_flat4(); 
     241  dp[31] = thick_flat5(); 
     242  dp[32] = thick_flat6(); 
     243  dp[33] = thick_flat7(); 
     244  dp[34] = thick_flat8(); 
     245  dp[35] = thick_flat9(); 
     246  dp[36] = thick_flat10(); 
     247 
     248  dp[37] = func_inter1(); 
     249  dp[38] = func_inter2(); 
     250  dp[39] = func_inter3(); 
     251  dp[40] = func_inter4(); 
     252  dp[41] = func_inter5(); 
     253  dp[42] = func_inter6(); 
     254  dp[43] = func_inter7(); 
     255  dp[44] = func_inter8(); 
     256  dp[45] = func_inter9(); 
     257  dp[46] = func_inter10(); 
     258 
     259  // Get the dispersion points for the radius 
     260  //vector<WeightPoint> weights_thick; 
     261  //thick_inter0.get_weights(weights_thick_inter0); 
     262 
     263 
     264  return re_kernel(dp,q); 
    138265} 
    139266 
    … …  
    145272 */ 
    146273double ReflModel :: operator()(double qx, double qy) { 
    147         // For 2D set qy as q, ignoring qx. 
    148         double q = qy;//sqrt(qx*qx + qy*qy); 
    149         if (q < 0.0){ 
    150                 return 0.0; 
    151         } 
    152         return (*this).operator()(q); 
     274  // For 2D set qy as q, ignoring qx. 
     275  double q = qy;//sqrt(qx*qx + qy*qy); 
     276  if (q < 0.0){ 
     277    return 0.0; 
     278  } 
     279  return (*this).operator()(q); 
    153280} 
    154281 
    … …  
    161288 */ 
    162289double ReflModel :: evaluate_rphi(double q, double phi) { 
    163         return (*this).operator()(q); 
     290  return (*this).operator()(q); 
    164291} 
    165292 
    … …  
    169296 */ 
    170297double ReflModel :: calculate_ER() { 
    171         //NOT implemented yet!!! 
    172         return 0.0; 
    173 } 
     298  //NOT implemented yet!!! 
     299  return 0.0; 
     300} 
Note: See TracChangeset for help on using the changeset viewer.