Changeset dd60b45 in sasview for sansmodels/src


Ignore:
Timestamp:
Jan 5, 2012 11:01:26 AM (13 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:
a8eab1c
Parents:
2d1b700
Message:

refl_adv model refactor

Location:
sansmodels/src
Files:
2 deleted
3 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/c_extensions/refl_adv.h

    r67424cd rdd60b45  
    33#if !defined(o_h) 
    44#define refl_adv_h 
     5#include "parameters.hh" 
    56 
    67/** 
    78 * Structure definition for sphere parameters 
    89 */ 
    9  //[PYTHONCLASS] = ReflAdvModel 
    10  //[DISP_PARAMS] = thick_inter0 
    11  //[DESCRIPTION] =<text>Calculate neutron reflectivity using the Parratt iterative formula 
    12  //                             Parameters: 
    13  //                             background:background 
    14  //                             scale: scale factor 
    15  //                             sld_bottom0: the SLD of the substrate 
    16  //                             sld_medium: the SLD of the incident medium 
    17  //                                     or superstrate 
    18  //                             sld_flatN: the SLD of the flat region of 
    19  //                                     the N'th layer 
    20  //                             thick_flatN: the thickness of the flat 
    21  //                                     region of the N'th layer 
    22  //                             func_interN: the function used to describe 
    23  //                                     the interface of the N'th layer 
    24  //                             nu_interN: the coefficient for the func_interN 
    25  //                             thick_interN: the thickness of the interface 
    26  //                                     of the N'th layer 
    27  //                             Note: the layer number starts to increase 
    28  //                                     from the bottom (substrate) to the top. 
    29  //             </text> 
    30  //[FIXED]=  <text></text> 
    31  //[NON_FITTABLE_PARAMS]= <text>n_layers;func_inter0;func_inter1;func_inter2;func_inter3;func_inter4;func_inter5;func_inter5;func_inter7;func_inter8;func_inter9;func_inter10 </text> 
    32  //[ORIENTATION_PARAMS]= <text> </text> 
    33  
    34 typedef struct { 
    35         /// number of layers 
    36         //  [DEFAULT]=n_layers=1 
    37         int n_layers; 
    38     /// Scale factor 
    39     //  [DEFAULT]=scale= 1.0 
    40         double scale; 
    41     /// thick_inter0 [A] 
    42     //  [DEFAULT]=thick_inter0=50.0 [A] 
    43         double thick_inter0; 
    44         ///     func_inter0 
    45         //  [DEFAULT]=func_inter0= 0 
    46         double func_inter0; 
    47         ///     sld_bottom0 [1/A^(2)] 
    48         //  [DEFAULT]=sld_bottom0= 2.07e-6 [1/A^(2)] 
    49         double sld_bottom0; 
    50         ///     sld_medium [1/A^(2)] 
    51         //  [DEFAULT]=sld_medium= 1.0e-6 [1/A^(2)] 
    52         double sld_medium; 
    53         /// Background 
    54         //  [DEFAULT]=background=0 
    55         double background; 
    56  
    57     //  [DEFAULT]=sld_flat1=4.0e-06 [1/A^(2)] 
    58     double sld_flat1; 
    59     //  [DEFAULT]=sld_flat2=3.5e-06 [1/A^(2)] 
    60     double sld_flat2; 
    61     //  [DEFAULT]=sld_flat3=4.0e-06 [1/A^(2)] 
    62     double sld_flat3; 
    63     //  [DEFAULT]=sld_flat4=3.5e-06 [1/A^(2)] 
    64     double sld_flat4; 
    65     //  [DEFAULT]=sld_flat5=4.0e-06 [1/A^(2)] 
    66     double sld_flat5; 
    67     //  [DEFAULT]=sld_flat6=3.5e-06 [1/A^(2)] 
    68     double sld_flat6; 
    69     //  [DEFAULT]=sld_flat7=4.0e-06 [1/A^(2)] 
    70     double sld_flat7; 
    71     //  [DEFAULT]=sld_flat8=3.5e-06 [1/A^(2)] 
    72     double sld_flat8; 
    73     //  [DEFAULT]=sld_flat9=4.0e-06 [1/A^(2)] 
    74     double sld_flat9; 
    75     //  [DEFAULT]=sld_flat10=3.5e-06 [1/A^(2)] 
    76     double sld_flat10; 
    77  
    78     //  [DEFAULT]=thick_inter1=50 [A] 
    79     double thick_inter1; 
    80     //  [DEFAULT]=thick_inter2=50 [A] 
    81     double thick_inter2; 
    82     //  [DEFAULT]=thick_inter3=50 [A] 
    83     double thick_inter3; 
    84     //  [DEFAULT]=thick_inter4=50 [A] 
    85     double thick_inter4; 
    86     //  [DEFAULT]=thick_inter5=50 [A] 
    87     double thick_inter5; 
    88     //  [DEFAULT]=thick_inter6=50 [A] 
    89     double thick_inter6; 
    90     //  [DEFAULT]=thick_inter7=50 [A] 
    91     double thick_inter7; 
    92     //  [DEFAULT]=thick_inter8=50 [A] 
    93     double thick_inter8; 
    94     //  [DEFAULT]=thick_inter9=50 [A] 
    95     double thick_inter9; 
    96     //  [DEFAULT]=thick_inter10=50 [A] 
    97     double thick_inter10; 
    98  
    99     //  [DEFAULT]=thick_flat1=100 [A] 
    100     double thick_flat1; 
    101     //  [DEFAULT]=thick_flat2=100 [A] 
    102     double thick_flat2; 
    103     //  [DEFAULT]=thick_flat3=100 [A] 
    104     double thick_flat3; 
    105     //  [DEFAULT]=thick_flat4=100 [A] 
    106     double thick_flat4; 
    107     //  [DEFAULT]=thick_flat5=100 [A] 
    108     double thick_flat5; 
    109     //  [DEFAULT]=thick_flat6=100 [A] 
    110     double thick_flat6; 
    111     //  [DEFAULT]=thick_flat7=100 [A] 
    112     double thick_flat7; 
    113     //  [DEFAULT]=thick_flat8=100 [A] 
    114     double thick_flat8; 
    115     //  [DEFAULT]=thick_flat9=100 [A] 
    116     double thick_flat9; 
    117     //  [DEFAULT]=thick_flat10=100 [A] 
    118     double thick_flat10; 
    119  
    120     //  [DEFAULT]=func_inter1=0 
    121     double func_inter1; 
    122     //  [DEFAULT]=func_inter2=0 
    123     double func_inter2; 
    124     //  [DEFAULT]=func_inter3=0 
    125     double func_inter3; 
    126     //  [DEFAULT]=func_inter4=0 
    127     double func_inter4; 
    128     //  [DEFAULT]=func_inter5=0 
    129     double func_inter5; 
    130     //  [DEFAULT]=func_inter6=0 
    131     double func_inter6; 
    132     //  [DEFAULT]=func_inter7=0 
    133     double func_inter7; 
    134     //  [DEFAULT]=func_inter8=0 
    135     double func_inter8; 
    136     //  [DEFAULT]=func_inter9=0 
    137     double func_inter9; 
    138     //  [DEFAULT]=func_inter10=0 
    139     double func_inter10; 
    140  
    141     //  [DEFAULT]=sldIM_flat1=0 [1/A^(2)] 
    142     double sldIM_flat1; 
    143     //  [DEFAULT]=sldIM_flat2=0 [1/A^(2)] 
    144     double sldIM_flat2; 
    145     //  [DEFAULT]=sldIM_flat3=0 [1/A^(2)] 
    146     double sldIM_flat3; 
    147     //  [DEFAULT]=sldIM_flat4=0 [1/A^(2)] 
    148     double sldIM_flat4; 
    149     //  [DEFAULT]=sldIM_flat5=0 [1/A^(2)] 
    150     double sldIM_flat5; 
    151     //  [DEFAULT]=sldIM_flat6=0 [1/A^(2)] 
    152     double sldIM_flat6; 
    153     //  [DEFAULT]=sldIM_flat7=0 [1/A^(2)] 
    154     double sldIM_flat7; 
    155     //  [DEFAULT]=sldIM_flat8=0 [1/A^(2)] 
    156     double sldIM_flat8; 
    157     //  [DEFAULT]=sldIM_flat9=0 [1/A^(2)] 
    158     double sldIM_flat9; 
    159     //  [DEFAULT]=sldIM_flat10=0 [1/A^(2)] 
    160     double sldIM_flat10; 
    161  
    162     //  [DEFAULT]=nu_inter1=2.5 
    163     double nu_inter1; 
    164     //  [DEFAULT]=nu_inter2=2.5 
    165     double nu_inter2; 
    166     //  [DEFAULT]=nu_inter3=2.5 
    167     double nu_inter3; 
    168     //  [DEFAULT]=nu_inter4=2.5 
    169     double nu_inter4; 
    170     //  [DEFAULT]=nu_inter5=2.5 
    171     double nu_inter5; 
    172     //  [DEFAULT]=nu_inter6=2.5 
    173     double nu_inter6; 
    174     //  [DEFAULT]=nu_inter7=2.5 
    175     double nu_inter7; 
    176     //  [DEFAULT]=nu_inter8=2.5 
    177     double nu_inter8; 
    178     //  [DEFAULT]=nu_inter9=2.5 
    179     double nu_inter9; 
    180     //  [DEFAULT]=nu_inter10=2.5 
    181     double nu_inter10; 
    182  
    183  
    184     //  [DEFAULT]=sldIM_sub0=0 
    185     double sldIM_sub0; 
    186     //  [DEFAULT]=sldIM_medium=0 
    187     double sldIM_medium; 
    188     //  [DEFAULT]=npts_inter=21.0 
    189     double npts_inter; 
    190     //  [DEFAULT]=nu_inter0=2.5 
    191     double nu_inter0; 
    192 } ReflAdvParameters; 
    193  
    194 double re_adv_kernel(double dq[], double q); 
    195  
    196 /// 1D scattering function 
    197 double refl_adv_analytical_1D(ReflAdvParameters *pars, double q); 
    198  
    199 /// 2D scattering function 
    200 double refl_adv_analytical_2D(ReflAdvParameters *pars, double q, double phi); 
    201 double refl_adv_analytical_2DXY(ReflAdvParameters *pars, double qx, double qy); 
     10//[PYTHONCLASS] = ReflAdvModel 
     11//[DISP_PARAMS] = thick_inter0 
     12//[DESCRIPTION] =<text>Calculate neutron reflectivity using the Parratt iterative formula 
     13//                              Parameters: 
     14//                              background:background 
     15//                              scale: scale factor 
     16//                              sld_bottom0: the SLD of the substrate 
     17//                              sld_medium: the SLD of the incident medium 
     18//                                      or superstrate 
     19//                              sld_flatN: the SLD of the flat region of 
     20//                                      the N'th layer 
     21//                              thick_flatN: the thickness of the flat 
     22//                                      region of the N'th layer 
     23//                              func_interN: the function used to describe 
     24//                                      the interface of the N'th layer 
     25//                              nu_interN: the coefficient for the func_interN 
     26//                              thick_interN: the thickness of the interface 
     27//                                      of the N'th layer 
     28//                              Note: the layer number starts to increase 
     29//                                      from the bottom (substrate) to the top. 
     30//              </text> 
     31//[FIXED]=  <text></text> 
     32//[NON_FITTABLE_PARAMS]= <text>n_layers;func_inter0;func_inter1;func_inter2;func_inter3;func_inter4;func_inter5;func_inter5;func_inter7;func_inter8;func_inter9;func_inter10 </text> 
     33//[ORIENTATION_PARAMS]= <text> </text> 
     34 
     35class ReflAdvModel{ 
     36public: 
     37  // Model parameters 
     38  /// number of layers 
     39  //  [DEFAULT]=n_layers=1 
     40  Parameter n_layers; 
     41  /// Scale factor 
     42  //  [DEFAULT]=scale= 1.0 
     43  Parameter scale; 
     44  /// thick_inter0 [A] 
     45  //  [DEFAULT]=thick_inter0=50.0 [A] 
     46  Parameter thick_inter0; 
     47  /// func_inter0 
     48  //  [DEFAULT]=func_inter0= 0 
     49  Parameter func_inter0; 
     50  /// sld_bottom0 [1/A^(2)] 
     51  //  [DEFAULT]=sld_bottom0= 2.07e-6 [1/A^(2)] 
     52  Parameter sld_bottom0; 
     53  /// sld_medium [1/A^(2)] 
     54  //  [DEFAULT]=sld_medium= 1.0e-6 [1/A^(2)] 
     55  Parameter sld_medium; 
     56  /// Background 
     57  //  [DEFAULT]=background=0 
     58  Parameter background; 
     59 
     60  //  [DEFAULT]=sld_flat1=4.0e-06 [1/A^(2)] 
     61  Parameter sld_flat1; 
     62  //  [DEFAULT]=sld_flat2=3.5e-06 [1/A^(2)] 
     63  Parameter sld_flat2; 
     64  //  [DEFAULT]=sld_flat3=4.0e-06 [1/A^(2)] 
     65  Parameter sld_flat3; 
     66  //  [DEFAULT]=sld_flat4=3.5e-06 [1/A^(2)] 
     67  Parameter sld_flat4; 
     68  //  [DEFAULT]=sld_flat5=4.0e-06 [1/A^(2)] 
     69  Parameter sld_flat5; 
     70  //  [DEFAULT]=sld_flat6=3.5e-06 [1/A^(2)] 
     71  Parameter sld_flat6; 
     72  //  [DEFAULT]=sld_flat7=4.0e-06 [1/A^(2)] 
     73  Parameter sld_flat7; 
     74  //  [DEFAULT]=sld_flat8=3.5e-06 [1/A^(2)] 
     75  Parameter sld_flat8; 
     76  //  [DEFAULT]=sld_flat9=4.0e-06 [1/A^(2)] 
     77  Parameter sld_flat9; 
     78  //  [DEFAULT]=sld_flat10=3.5e-06 [1/A^(2)] 
     79  Parameter sld_flat10; 
     80 
     81  //  [DEFAULT]=thick_inter1=50 [A] 
     82  Parameter thick_inter1; 
     83  //  [DEFAULT]=thick_inter2=50 [A] 
     84  Parameter thick_inter2; 
     85  //  [DEFAULT]=thick_inter3=50 [A] 
     86  Parameter thick_inter3; 
     87  //  [DEFAULT]=thick_inter4=50 [A] 
     88  Parameter thick_inter4; 
     89  //  [DEFAULT]=thick_inter5=50 [A] 
     90  Parameter thick_inter5; 
     91  //  [DEFAULT]=thick_inter6=50 [A] 
     92  Parameter thick_inter6; 
     93  //  [DEFAULT]=thick_inter7=50 [A] 
     94  Parameter thick_inter7; 
     95  //  [DEFAULT]=thick_inter8=50 [A] 
     96  Parameter thick_inter8; 
     97  //  [DEFAULT]=thick_inter9=50 [A] 
     98  Parameter thick_inter9; 
     99  //  [DEFAULT]=thick_inter10=50 [A] 
     100  Parameter thick_inter10; 
     101 
     102  //  [DEFAULT]=thick_flat1=100 [A] 
     103  Parameter thick_flat1; 
     104  //  [DEFAULT]=thick_flat2=100 [A] 
     105  Parameter thick_flat2; 
     106  //  [DEFAULT]=thick_flat3=100 [A] 
     107  Parameter thick_flat3; 
     108  //  [DEFAULT]=thick_flat4=100 [A] 
     109  Parameter thick_flat4; 
     110  //  [DEFAULT]=thick_flat5=100 [A] 
     111  Parameter thick_flat5; 
     112  //  [DEFAULT]=thick_flat6=100 [A] 
     113  Parameter thick_flat6; 
     114  //  [DEFAULT]=thick_flat7=100 [A] 
     115  Parameter thick_flat7; 
     116  //  [DEFAULT]=thick_flat8=100 [A] 
     117  Parameter thick_flat8; 
     118  //  [DEFAULT]=thick_flat9=100 [A] 
     119  Parameter thick_flat9; 
     120  //  [DEFAULT]=thick_flat10=100 [A] 
     121  Parameter thick_flat10; 
     122 
     123  //  [DEFAULT]=func_inter1=0 
     124  Parameter func_inter1; 
     125  //  [DEFAULT]=func_inter2=0 
     126  Parameter func_inter2; 
     127  //  [DEFAULT]=func_inter3=0 
     128  Parameter func_inter3; 
     129  //  [DEFAULT]=func_inter4=0 
     130  Parameter func_inter4; 
     131  //  [DEFAULT]=func_inter5=0 
     132  Parameter func_inter5; 
     133  //  [DEFAULT]=func_inter6=0 
     134  Parameter func_inter6; 
     135  //  [DEFAULT]=func_inter7=0 
     136  Parameter func_inter7; 
     137  //  [DEFAULT]=func_inter8=0 
     138  Parameter func_inter8; 
     139  //  [DEFAULT]=func_inter9=0 
     140  Parameter func_inter9; 
     141  //  [DEFAULT]=func_inter10=0 
     142  Parameter func_inter10; 
     143 
     144  //  [DEFAULT]=sldIM_flat1=0 [1/A^(2)] 
     145  Parameter sldIM_flat1; 
     146  //  [DEFAULT]=sldIM_flat2=0 [1/A^(2)] 
     147  Parameter sldIM_flat2; 
     148  //  [DEFAULT]=sldIM_flat3=0 [1/A^(2)] 
     149  Parameter sldIM_flat3; 
     150  //  [DEFAULT]=sldIM_flat4=0 [1/A^(2)] 
     151  Parameter sldIM_flat4; 
     152  //  [DEFAULT]=sldIM_flat5=0 [1/A^(2)] 
     153  Parameter sldIM_flat5; 
     154  //  [DEFAULT]=sldIM_flat6=0 [1/A^(2)] 
     155  Parameter sldIM_flat6; 
     156  //  [DEFAULT]=sldIM_flat7=0 [1/A^(2)] 
     157  Parameter sldIM_flat7; 
     158  //  [DEFAULT]=sldIM_flat8=0 [1/A^(2)] 
     159  Parameter sldIM_flat8; 
     160  //  [DEFAULT]=sldIM_flat9=0 [1/A^(2)] 
     161  Parameter sldIM_flat9; 
     162  //  [DEFAULT]=sldIM_flat10=0 [1/A^(2)] 
     163  Parameter sldIM_flat10; 
     164 
     165  //  [DEFAULT]=nu_inter1=2.5 
     166  Parameter nu_inter1; 
     167  //  [DEFAULT]=nu_inter2=2.5 
     168  Parameter nu_inter2; 
     169  //  [DEFAULT]=nu_inter3=2.5 
     170  Parameter nu_inter3; 
     171  //  [DEFAULT]=nu_inter4=2.5 
     172  Parameter nu_inter4; 
     173  //  [DEFAULT]=nu_inter5=2.5 
     174  Parameter nu_inter5; 
     175  //  [DEFAULT]=nu_inter6=2.5 
     176  Parameter nu_inter6; 
     177  //  [DEFAULT]=nu_inter7=2.5 
     178  Parameter nu_inter7; 
     179  //  [DEFAULT]=nu_inter8=2.5 
     180  Parameter nu_inter8; 
     181  //  [DEFAULT]=nu_inter9=2.5 
     182  Parameter nu_inter9; 
     183  //  [DEFAULT]=nu_inter10=2.5 
     184  Parameter nu_inter10; 
     185 
     186 
     187  //  [DEFAULT]=sldIM_sub0=0 
     188  Parameter sldIM_sub0; 
     189  //  [DEFAULT]=sldIM_medium=0 
     190  Parameter sldIM_medium; 
     191  //  [DEFAULT]=npts_inter=21.0 
     192  Parameter npts_inter; 
     193  //  [DEFAULT]=nu_inter0=2.5 
     194  Parameter nu_inter0; 
     195 
     196  // Constructor 
     197  ReflAdvModel(); 
     198 
     199  // Operators to get I(Q) 
     200  double operator()(double q); 
     201  double operator()(double qx, double qy); 
     202  double calculate_ER(); 
     203  double evaluate_rphi(double q, double phi); 
     204}; 
    202205 
    203206#endif 

  • sansmodels/src/c_models/models.hh

    r2d1b700 rdd60b45  
    2727using namespace std; 
    2828 
    29 class ReflAdvModel{ 
    30 public: 
    31         // Model parameters 
    32         Parameter n_layers; 
    33         Parameter scale; 
    34         Parameter thick_inter0; 
    35         Parameter func_inter0; 
    36         Parameter sld_bottom0; 
    37         Parameter sld_medium; 
    38         Parameter background; 
    39  
    40         Parameter sld_flat1; 
    41         Parameter sld_flat2; 
    42         Parameter sld_flat3; 
    43         Parameter sld_flat4; 
    44         Parameter sld_flat5; 
    45         Parameter sld_flat6; 
    46         Parameter sld_flat7; 
    47         Parameter sld_flat8; 
    48         Parameter sld_flat9; 
    49         Parameter sld_flat10; 
    50  
    51         Parameter thick_inter1; 
    52         Parameter thick_inter2; 
    53         Parameter thick_inter3; 
    54         Parameter thick_inter4; 
    55         Parameter thick_inter5; 
    56         Parameter thick_inter6; 
    57         Parameter thick_inter7; 
    58         Parameter thick_inter8; 
    59         Parameter thick_inter9; 
    60         Parameter thick_inter10; 
    61  
    62         Parameter thick_flat1; 
    63         Parameter thick_flat2; 
    64         Parameter thick_flat3; 
    65         Parameter thick_flat4; 
    66         Parameter thick_flat5; 
    67         Parameter thick_flat6; 
    68         Parameter thick_flat7; 
    69         Parameter thick_flat8; 
    70         Parameter thick_flat9; 
    71         Parameter thick_flat10; 
    72  
    73         Parameter func_inter1; 
    74         Parameter func_inter2; 
    75         Parameter func_inter3; 
    76         Parameter func_inter4; 
    77         Parameter func_inter5; 
    78         Parameter func_inter6; 
    79         Parameter func_inter7; 
    80         Parameter func_inter8; 
    81         Parameter func_inter9; 
    82         Parameter func_inter10; 
    83  
    84         Parameter sldIM_flat1; 
    85         Parameter sldIM_flat2; 
    86         Parameter sldIM_flat3; 
    87         Parameter sldIM_flat4; 
    88         Parameter sldIM_flat5; 
    89         Parameter sldIM_flat6; 
    90         Parameter sldIM_flat7; 
    91         Parameter sldIM_flat8; 
    92         Parameter sldIM_flat9; 
    93         Parameter sldIM_flat10; 
    94  
    95         Parameter nu_inter1; 
    96         Parameter nu_inter2; 
    97         Parameter nu_inter3; 
    98         Parameter nu_inter4; 
    99         Parameter nu_inter5; 
    100         Parameter nu_inter6; 
    101         Parameter nu_inter7; 
    102         Parameter nu_inter8; 
    103         Parameter nu_inter9; 
    104         Parameter nu_inter10; 
    105  
    106         Parameter sldIM_sub0; 
    107         Parameter sldIM_medium; 
    108         Parameter npts_inter; 
    109         Parameter nu_inter0; 
    110  
    111         // Constructor 
    112         ReflAdvModel(); 
    113  
    114         // Operators to get I(Q) 
    115         double operator()(double q); 
    116         double operator()(double qx, double qy); 
    117         double calculate_ER(); 
    118         double evaluate_rphi(double q, double phi); 
    119 }; 
     29 
    12030 
    12131 

  • sansmodels/src/c_models/refl_adv.cpp

    r67424cd rdd60b45  
    11 
    22#include <math.h> 
    3 #include "models.hh" 
    43#include "parameters.hh" 
    54#include <stdio.h> 
     5#include <stdlib.h> 
     6#include "refl_adv.h" 
     7 
     8extern "C" { 
     9#include "libmultifunc/librefl.h" 
     10} 
     11 
    612using namespace std; 
    713 
    8 extern "C" { 
    9         #include "refl_adv.h" 
    10 } 
    11  
     14#define lamda 4.62 
     15 
     16 
     17double re_adv_kernel(double dp[], double q) { 
     18  int n = dp[0]; 
     19  int i,j; 
     20  double nsl; 
     21 
     22  double scale = dp[1]; 
     23  double thick_inter_sub = dp[2]; 
     24  double sld_sub = dp[4]; 
     25  double sld_super = dp[5]; 
     26  double background = dp[6]; 
     27  double npts = dp[69]; //number of sub_layers in each interface 
     28 
     29  double total_thick=0.0; 
     30 
     31  int n_s; 
     32  double sld_i,sldim_i,dz,phi,R,ko2; 
     33  double pi; 
     34 
     35  int* fun_type; 
     36  double* sld; 
     37  double* sld_im; 
     38  double* thick_inter; 
     39  double* thick; 
     40  double* fun_coef; 
     41  complex  phi1,alpha,alpha2,kn,fnm,fnp,rn,Xn,nn,nn2,an,nnp1,one,two,n_sub,n_sup,knp1,Xnp1; 
     42 
     43  fun_type = (int*)malloc((n+2)*sizeof(int)); 
     44  sld = (double*)malloc((n+2)*sizeof(double)); 
     45  sld_im = (double*)malloc((n+2)*sizeof(double)); 
     46  thick_inter = (double*)malloc((n+2)*sizeof(double)); 
     47  thick = (double*)malloc((n+2)*sizeof(double)); 
     48  fun_coef = (double*)malloc((n+2)*sizeof(double)); 
     49 
     50  fun_type[0] = dp[3]; 
     51  fun_coef[0] = fabs(dp[70]); 
     52  for (i =1; i<=n; i++){ 
     53    sld[i] = dp[i+6]; 
     54    thick_inter[i]= dp[i+16]; 
     55    thick[i] = dp[i+26]; 
     56    fun_type[i] = dp[i+36]; 
     57    sld_im[i] = dp[i+46]; 
     58    fun_coef[i] = fabs(dp[i+56]); 
     59    //printf("type_func2 =%g\n",fun_coef[i]); 
     60 
     61    total_thick += thick[i] + thick_inter[i]; 
     62  } 
     63  sld[0] = sld_sub; 
     64  sld[n+1] = sld_super; 
     65  sld_im[0] = fabs(dp[1+66]); 
     66  sld_im[n+1] = fabs(dp[2+66]); 
     67  thick[0] = total_thick/5.0; 
     68  thick[n+1] = total_thick/5.0; 
     69  thick_inter[0] = thick_inter_sub; 
     70  thick_inter[n+1] = 0.0; 
     71  fun_coef[n+1] = 0.0; 
     72 
     73  nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now 
     74 
     75  pi = 4.0*atan(1.0); 
     76  one = cassign(1.0,0.0); 
     77  Xn = cassign(0.0,0.0); 
     78  two = cassign(0.0,-2.0); 
     79 
     80  //Checking if floor is available. 
     81  //no imaginary sld inputs in this function yet 
     82  n_sub=cassign(1.0-sld_sub*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sld_im[0]); 
     83  n_sup=cassign(1.0-sld_super*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sld_im[n+1]); 
     84  ko2 = pow(2.0*pi/lamda,2.0); 
     85 
     86  phi = asin(lamda*q/(4.0*pi)); 
     87  phi1 = cplx_div(rcmult(phi,one),n_sup); 
     88  alpha = cplx_mult(n_sup,cplx_cos(phi1)); 
     89  alpha2 = cplx_mult(alpha,alpha); 
     90 
     91  nnp1=n_sub; 
     92  knp1=cplx_sqrt(rcmult(ko2,cplx_sub(cplx_mult(nnp1,nnp1),alpha2)));  //nnp1*ko*sin(phinp1) 
     93  Xnp1=cassign(0.0,0.0); 
     94  dz = 0.0; 
     95  // iteration for # of layers +sub from the top 
     96  for (i=1;i<=n+1; i++){ 
     97    //if (fun_coef[i]==0.0) 
     98    //  // this condition protects an error in numerical multiplication 
     99    //  fun_coef[i] = 1e-14; 
     100    //iteration for 9 sub-layers 
     101    for (j=0;j<2;j++){ 
     102      for (n_s=0;n_s<nsl; n_s++){ 
     103        // for flat layer 
     104        if (j==1){ 
     105          if (i==n+1) 
     106            break; 
     107          dz = thick[i]; 
     108          sld_i = sld[i]; 
     109          sldim_i = sld_im[i]; 
     110        } 
     111        // for interface 
     112        else{ 
     113          dz = thick_inter[i-1]/nsl; 
     114          if (sld[i-1] == sld[i]){ 
     115            sld_i = sld[i]; 
     116          } 
     117          else{ 
     118            sld_i = intersldfunc(fun_type[i-1],nsl, n_s+0.5, fun_coef[i-1], sld[i-1], sld[i]); 
     119          } 
     120          if (sld_im[i-1] == sld_im[i]){ 
     121            sldim_i = sld_im[i]; 
     122          } 
     123          else{ 
     124            sldim_i = intersldfunc(fun_type[i-1],nsl, n_s+0.5, fun_coef[i-1], sld_im[i-1], sld_im[i]); 
     125          } 
     126        } 
     127        nn = cassign(1.0-sld_i*pow(lamda,2.0)/(2.0*pi),pow(lamda,2.0)/(2.0*pi)*sldim_i); 
     128        nn2=cplx_mult(nn,nn); 
     129 
     130        kn=cplx_sqrt(rcmult(ko2,cplx_sub(nn2,alpha2)));        //nn*ko*sin(phin) 
     131        an=cplx_exp(rcmult(dz,cplx_mult(two,kn))); 
     132 
     133        fnm=cplx_sub(kn,knp1); 
     134        fnp=cplx_add(kn,knp1); 
     135        rn=cplx_div(fnm,fnp); 
     136        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)) 
     137 
     138        Xnp1=Xn; 
     139        knp1=kn; 
     140        // no for-loop for flat layer 
     141        if (j==1) 
     142          break; 
     143      } 
     144    } 
     145  } 
     146  R=pow(Xn.re,2.0)+pow(Xn.im,2.0); 
     147  // This temperarily fixes the total reflection for Rfunction and linear. 
     148  // ToDo: Show why it happens that Xn.re=0 and Xn.im >1! 
     149  if (Xn.im == 0.0 || R > 1){ 
     150    R=1.0; 
     151  } 
     152  R *= scale; 
     153  R += background; 
     154 
     155  free(fun_type); 
     156  free(sld); 
     157  free(sld_im); 
     158  free(thick_inter); 
     159  free(thick); 
     160  free(fun_coef); 
     161 
     162  return R; 
     163 
     164} 
    12165ReflAdvModel :: ReflAdvModel() { 
    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(2.7e-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.0); 
    35     thick_inter2 = Parameter(1.0); 
    36     thick_inter3 = Parameter(1.0); 
    37     thick_inter4 = Parameter(1.0); 
    38     thick_inter5 = Parameter(1.0); 
    39     thick_inter6 = Parameter(1.0); 
    40     thick_inter7 = Parameter(1.0); 
    41     thick_inter8 = Parameter(1.0); 
    42     thick_inter9 = Parameter(1.0); 
    43     thick_inter10 = Parameter(1.0); 
    44  
    45  
    46     thick_flat1 = Parameter(15.0); 
    47     thick_flat2 = Parameter(100.0); 
    48     thick_flat3 = Parameter(100.0); 
    49     thick_flat4 = Parameter(100.0); 
    50     thick_flat5 = Parameter(100.0); 
    51     thick_flat6 = Parameter(100.0); 
    52     thick_flat7 = Parameter(100.0); 
    53     thick_flat8 = Parameter(100.0); 
    54     thick_flat9 = Parameter(100.0); 
    55     thick_flat10 = Parameter(100.0); 
    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); 
    68  
    69     sldIM_flat1 = Parameter(0); 
    70     sldIM_flat2 = Parameter(0); 
    71     sldIM_flat3 = Parameter(0); 
    72     sldIM_flat4 = Parameter(0); 
    73     sldIM_flat5 = Parameter(0); 
    74     sldIM_flat6 = Parameter(0); 
    75     sldIM_flat7 = Parameter(0); 
    76     sldIM_flat8 = Parameter(0); 
    77     sldIM_flat9 = Parameter(0); 
    78     sldIM_flat10 = Parameter(0); 
    79  
    80     nu_inter1 = Parameter(2.5); 
    81     nu_inter2 = Parameter(2.5); 
    82     nu_inter3 = Parameter(2.5); 
    83     nu_inter4 = Parameter(2.5); 
    84     nu_inter5 = Parameter(2.5); 
    85     nu_inter6 = Parameter(2.5); 
    86     nu_inter7 = Parameter(2.5); 
    87     nu_inter8 = Parameter(2.5); 
    88     nu_inter9 = Parameter(2.5); 
    89     nu_inter10 = Parameter(2.5); 
    90  
    91     sldIM_sub0 = Parameter(0.0); 
    92     sldIM_medium = Parameter(0.0); 
    93     npts_inter = Parameter(21.0); 
    94     nu_inter0 = Parameter(2.5); 
     166  n_layers = Parameter(1.0); 
     167  scale = Parameter(1.0); 
     168  thick_inter0 = Parameter(1.0); 
     169  func_inter0 = Parameter(0); 
     170  sld_bottom0 = Parameter(2.07e-06); 
     171  sld_medium = Parameter(1.0e-06); 
     172  background = Parameter(0.0); 
     173 
     174 
     175  sld_flat1 = Parameter(2.7e-06); 
     176  sld_flat2 = Parameter(3.5e-06); 
     177  sld_flat3 = Parameter(4.0e-06); 
     178  sld_flat4 = Parameter(3.5e-06); 
     179  sld_flat5 = Parameter(4.0e-06); 
     180  sld_flat6 = Parameter(3.5e-06); 
     181  sld_flat7 = Parameter(4.0e-06); 
     182  sld_flat8 = Parameter(3.5e-06); 
     183  sld_flat9 = Parameter(4.0e-06); 
     184  sld_flat10 = Parameter(3.5e-06); 
     185 
     186 
     187  thick_inter1 = Parameter(1.0); 
     188  thick_inter2 = Parameter(1.0); 
     189  thick_inter3 = Parameter(1.0); 
     190  thick_inter4 = Parameter(1.0); 
     191  thick_inter5 = Parameter(1.0); 
     192  thick_inter6 = Parameter(1.0); 
     193  thick_inter7 = Parameter(1.0); 
     194  thick_inter8 = Parameter(1.0); 
     195  thick_inter9 = Parameter(1.0); 
     196  thick_inter10 = Parameter(1.0); 
     197 
     198 
     199  thick_flat1 = Parameter(15.0); 
     200  thick_flat2 = Parameter(100.0); 
     201  thick_flat3 = Parameter(100.0); 
     202  thick_flat4 = Parameter(100.0); 
     203  thick_flat5 = Parameter(100.0); 
     204  thick_flat6 = Parameter(100.0); 
     205  thick_flat7 = Parameter(100.0); 
     206  thick_flat8 = Parameter(100.0); 
     207  thick_flat9 = Parameter(100.0); 
     208  thick_flat10 = Parameter(100.0); 
     209 
     210 
     211  func_inter1 = Parameter(0); 
     212  func_inter2 = Parameter(0); 
     213  func_inter3 = Parameter(0); 
     214  func_inter4 = Parameter(0); 
     215  func_inter5 = Parameter(0); 
     216  func_inter6 = Parameter(0); 
     217  func_inter7 = Parameter(0); 
     218  func_inter8 = Parameter(0); 
     219  func_inter9 = Parameter(0); 
     220  func_inter10 = Parameter(0); 
     221 
     222  sldIM_flat1 = Parameter(0); 
     223  sldIM_flat2 = Parameter(0); 
     224  sldIM_flat3 = Parameter(0); 
     225  sldIM_flat4 = Parameter(0); 
     226  sldIM_flat5 = Parameter(0); 
     227  sldIM_flat6 = Parameter(0); 
     228  sldIM_flat7 = Parameter(0); 
     229  sldIM_flat8 = Parameter(0); 
     230  sldIM_flat9 = Parameter(0); 
     231  sldIM_flat10 = Parameter(0); 
     232 
     233  nu_inter1 = Parameter(2.5); 
     234  nu_inter2 = Parameter(2.5); 
     235  nu_inter3 = Parameter(2.5); 
     236  nu_inter4 = Parameter(2.5); 
     237  nu_inter5 = Parameter(2.5); 
     238  nu_inter6 = Parameter(2.5); 
     239  nu_inter7 = Parameter(2.5); 
     240  nu_inter8 = Parameter(2.5); 
     241  nu_inter9 = Parameter(2.5); 
     242  nu_inter10 = Parameter(2.5); 
     243 
     244  sldIM_sub0 = Parameter(0.0); 
     245  sldIM_medium = Parameter(0.0); 
     246  npts_inter = Parameter(21.0); 
     247  nu_inter0 = Parameter(2.5); 
    95248} 
    96249 
    … …  
    101254 */ 
    102255double ReflAdvModel :: operator()(double q) { 
    103         double dp[71]; 
    104         // Fill parameter array for IGOR library 
    105         // Add the background after averaging 
    106         dp[0] = n_layers(); 
    107         dp[1] = scale(); 
    108         dp[2] = thick_inter0(); 
    109         dp[3] = func_inter0(); 
    110         dp[4] = sld_bottom0(); 
    111         dp[5] = sld_medium(); 
    112         dp[6] = background(); 
    113  
    114         dp[7] = sld_flat1(); 
    115         dp[8] = sld_flat2(); 
    116         dp[9] = sld_flat3(); 
    117         dp[10] = sld_flat4(); 
    118         dp[11] = sld_flat5(); 
    119         dp[12] = sld_flat6(); 
    120         dp[13] = sld_flat7(); 
    121         dp[14] = sld_flat8(); 
    122         dp[15] = sld_flat9(); 
    123         dp[16] = sld_flat10(); 
    124  
    125         dp[17] = thick_inter1(); 
    126         dp[18] = thick_inter2(); 
    127         dp[19] = thick_inter3(); 
    128         dp[20] = thick_inter4(); 
    129         dp[21] = thick_inter5(); 
    130         dp[22] = thick_inter6(); 
    131         dp[23] = thick_inter7(); 
    132         dp[24] = thick_inter8(); 
    133         dp[25] = thick_inter9(); 
    134         dp[26] = thick_inter10(); 
    135  
    136         dp[27] = thick_flat1(); 
    137         dp[28] = thick_flat2(); 
    138         dp[29] = thick_flat3(); 
    139         dp[30] = thick_flat4(); 
    140         dp[31] = thick_flat5(); 
    141         dp[32] = thick_flat6(); 
    142         dp[33] = thick_flat7(); 
    143         dp[34] = thick_flat8(); 
    144         dp[35] = thick_flat9(); 
    145         dp[36] = thick_flat10(); 
    146  
    147         dp[37] = func_inter1(); 
    148         dp[38] = func_inter2(); 
    149         dp[39] = func_inter3(); 
    150         dp[40] = func_inter4(); 
    151         dp[41] = func_inter5(); 
    152         dp[42] = func_inter6(); 
    153         dp[43] = func_inter7(); 
    154         dp[44] = func_inter8(); 
    155         dp[45] = func_inter9(); 
    156         dp[46] = func_inter10(); 
    157  
    158         dp[47] = sldIM_flat1(); 
    159         dp[48] = sldIM_flat2(); 
    160         dp[49] = sldIM_flat3(); 
    161         dp[50] = sldIM_flat4(); 
    162         dp[51] = sldIM_flat5(); 
    163         dp[52] = sldIM_flat6(); 
    164         dp[53] = sldIM_flat7(); 
    165         dp[54] = sldIM_flat8(); 
    166         dp[55] = sldIM_flat9(); 
    167         dp[56] = sldIM_flat10(); 
    168  
    169         dp[57] = nu_inter1(); 
    170         dp[58] = nu_inter2(); 
    171         dp[59] = nu_inter3(); 
    172         dp[60] = nu_inter4(); 
    173         dp[61] = nu_inter5(); 
    174         dp[62] = nu_inter6(); 
    175         dp[63] = nu_inter7(); 
    176         dp[64] = nu_inter8(); 
    177         dp[65] = nu_inter9(); 
    178         dp[66] = nu_inter10(); 
    179  
    180         dp[67] = sldIM_sub0(); 
    181         dp[68] = sldIM_medium(); 
    182         dp[69] = npts_inter(); 
    183         dp[70] = nu_inter0(); 
    184         // Get the dispersion points for the radius 
    185         //vector<WeightPoint> weights_thick; 
    186         //thick_inter0.get_weights(weights_thick_inter0); 
    187  
    188  
    189         return re_adv_kernel(dp,q); 
     256  double dp[71]; 
     257  // Fill parameter array for IGOR library 
     258  // Add the background after averaging 
     259  dp[0] = n_layers(); 
     260  dp[1] = scale(); 
     261  dp[2] = thick_inter0(); 
     262  dp[3] = func_inter0(); 
     263  dp[4] = sld_bottom0(); 
     264  dp[5] = sld_medium(); 
     265  dp[6] = background(); 
     266 
     267  dp[7] = sld_flat1(); 
     268  dp[8] = sld_flat2(); 
     269  dp[9] = sld_flat3(); 
     270  dp[10] = sld_flat4(); 
     271  dp[11] = sld_flat5(); 
     272  dp[12] = sld_flat6(); 
     273  dp[13] = sld_flat7(); 
     274  dp[14] = sld_flat8(); 
     275  dp[15] = sld_flat9(); 
     276  dp[16] = sld_flat10(); 
     277 
     278  dp[17] = thick_inter1(); 
     279  dp[18] = thick_inter2(); 
     280  dp[19] = thick_inter3(); 
     281  dp[20] = thick_inter4(); 
     282  dp[21] = thick_inter5(); 
     283  dp[22] = thick_inter6(); 
     284  dp[23] = thick_inter7(); 
     285  dp[24] = thick_inter8(); 
     286  dp[25] = thick_inter9(); 
     287  dp[26] = thick_inter10(); 
     288 
     289  dp[27] = thick_flat1(); 
     290  dp[28] = thick_flat2(); 
     291  dp[29] = thick_flat3(); 
     292  dp[30] = thick_flat4(); 
     293  dp[31] = thick_flat5(); 
     294  dp[32] = thick_flat6(); 
     295  dp[33] = thick_flat7(); 
     296  dp[34] = thick_flat8(); 
     297  dp[35] = thick_flat9(); 
     298  dp[36] = thick_flat10(); 
     299 
     300  dp[37] = func_inter1(); 
     301  dp[38] = func_inter2(); 
     302  dp[39] = func_inter3(); 
     303  dp[40] = func_inter4(); 
     304  dp[41] = func_inter5(); 
     305  dp[42] = func_inter6(); 
     306  dp[43] = func_inter7(); 
     307  dp[44] = func_inter8(); 
     308  dp[45] = func_inter9(); 
     309  dp[46] = func_inter10(); 
     310 
     311  dp[47] = sldIM_flat1(); 
     312  dp[48] = sldIM_flat2(); 
     313  dp[49] = sldIM_flat3(); 
     314  dp[50] = sldIM_flat4(); 
     315  dp[51] = sldIM_flat5(); 
     316  dp[52] = sldIM_flat6(); 
     317  dp[53] = sldIM_flat7(); 
     318  dp[54] = sldIM_flat8(); 
     319  dp[55] = sldIM_flat9(); 
     320  dp[56] = sldIM_flat10(); 
     321 
     322  dp[57] = nu_inter1(); 
     323  dp[58] = nu_inter2(); 
     324  dp[59] = nu_inter3(); 
     325  dp[60] = nu_inter4(); 
     326  dp[61] = nu_inter5(); 
     327  dp[62] = nu_inter6(); 
     328  dp[63] = nu_inter7(); 
     329  dp[64] = nu_inter8(); 
     330  dp[65] = nu_inter9(); 
     331  dp[66] = nu_inter10(); 
     332 
     333  dp[67] = sldIM_sub0(); 
     334  dp[68] = sldIM_medium(); 
     335  dp[69] = npts_inter(); 
     336  dp[70] = nu_inter0(); 
     337  // Get the dispersion points for the radius 
     338  //vector<WeightPoint> weights_thick; 
     339  //thick_inter0.get_weights(weights_thick_inter0); 
     340 
     341 
     342  return re_adv_kernel(dp,q); 
    190343} 
    191344 
    … …  
    197350 */ 
    198351double ReflAdvModel :: operator()(double qx, double qy) { 
    199         // For 2D set qy as q, ignoring qx. 
    200         double q = qy;//sqrt(qx*qx + qy*qy); 
    201         if (q < 0.0){ 
    202                 return 0.0; 
    203         } 
    204         return (*this).operator()(q); 
     352  // For 2D set qy as q, ignoring qx. 
     353  double q = qy;//sqrt(qx*qx + qy*qy); 
     354  if (q < 0.0){ 
     355    return 0.0; 
     356  } 
     357  return (*this).operator()(q); 
    205358} 
    206359 
    … …  
    213366 */ 
    214367double ReflAdvModel :: evaluate_rphi(double q, double phi) { 
    215         return (*this).operator()(q); 
     368  return (*this).operator()(q); 
    216369} 
    217370 
    … …  
    221374 */ 
    222375double ReflAdvModel :: calculate_ER() { 
    223         //NOT implemented yet!!! 
    224         return 0.0; 
    225 } 
     376  //NOT implemented yet!!! 
     377  return 0.0; 
     378} 
Note: See TracChangeset for help on using the changeset viewer.