Changeset 7ffa8196 in sasview


Ignore:
Timestamp:
Jan 5, 2012 8:25:27 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:
2d1b700
Parents:
c637521
Message:

refactor RPA model

Location:
sansmodels/src
Files:
1 deleted
4 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/c_extensions/rpa.h

    r67424cd r7ffa8196  
    11#if !defined(o_h) 
    22#define rpa_h 
     3#include "parameters.hh" 
    34 
    45/** 
    56 * Structure definition for sphere parameters 
    67 */ 
    7  //[PYTHONCLASS] = RPAModel 
    8  //[DISP_PARAMS] = background 
    9  //[DESCRIPTION] =<text>  THIS FORMALISM APPLIES TO MULTICOMPONENT POLYMER MIXTURES IN THE 
    10  //                                                     HOMOGENEOUS (MIXED) PHASE REGION ONLY.; 
    11  //                                                     CASE 0: C/D BINARY MIXTURE OF HOMOPOLYMERS 
    12  //                                                     CASE 1: C-D DIBLOCK COPOLYMER 
    13  //                                                     CASE 2: B/C/D TERNARY MIXTURE OF HOMOPOLYMERS 
    14  //                                                     CASE 3: B/C-D MIXTURE OF HOMOPOLYMER B AND 
    15  //                                                             DIBLOCK COPOLYMER C-D 
    16  //                                                     CASE 4: B-C-D TRIBLOCK COPOLYMER 
    17  //                                                     CASE 5: A/B/C/D QUATERNARY MIXTURE OF HOMOPOLYMERS 
    18  //                                                     CASE 6: A/B/C-D MIXTURE OF TWO HOMOPOLYMERS A/B 
    19  //                                                             AND A DIBLOCK C-D 
    20  //                                                     CASE 7: A/B-C-D MIXTURE OF A HOMOPOLYMER A AND A 
    21  //                                                             TRIBLOCK B-C-D 
    22  //                                                     CASE 8: A-B/C-D MIXTURE OF TWO DIBLOCK COPOLYMERS 
    23  //                                                             A-B AND C-D 
    24  //                                                     CASE 9: A-B-C-D FOUR-BLOCK COPOLYMER 
    25  //                                                     See details in the model function help 
    26  //             </text> 
    27  //[FIXED]=  <text></text> 
    28  //[NON_FITTABLE_PARAMS]= <text> lcase_n; Na; Phia; va; La; Nb; Phib; vb; Lb;Nc; Phic; vc; Lc;Nd; Phid; vd; Ld; </text> 
    29  //[ORIENTATION_PARAMS]= <text> </text> 
     8//[PYTHONCLASS] = RPAModel 
     9//[DISP_PARAMS] = background 
     10//[DESCRIPTION] =<text>  THIS FORMALISM APPLIES TO MULTICOMPONENT POLYMER MIXTURES IN THE 
     11//                                                      HOMOGENEOUS (MIXED) PHASE REGION ONLY.; 
     12//                                                      CASE 0: C/D BINARY MIXTURE OF HOMOPOLYMERS 
     13//                                                      CASE 1: C-D DIBLOCK COPOLYMER 
     14//                                                      CASE 2: B/C/D TERNARY MIXTURE OF HOMOPOLYMERS 
     15//                                                      CASE 3: B/C-D MIXTURE OF HOMOPOLYMER B AND 
     16//                                                              DIBLOCK COPOLYMER C-D 
     17//                                                      CASE 4: B-C-D TRIBLOCK COPOLYMER 
     18//                                                      CASE 5: A/B/C/D QUATERNARY MIXTURE OF HOMOPOLYMERS 
     19//                                                      CASE 6: A/B/C-D MIXTURE OF TWO HOMOPOLYMERS A/B 
     20//                                                              AND A DIBLOCK C-D 
     21//                                                      CASE 7: A/B-C-D MIXTURE OF A HOMOPOLYMER A AND A 
     22//                                                              TRIBLOCK B-C-D 
     23//                                                      CASE 8: A-B/C-D MIXTURE OF TWO DIBLOCK COPOLYMERS 
     24//                                                              A-B AND C-D 
     25//                                                      CASE 9: A-B-C-D FOUR-BLOCK COPOLYMER 
     26//                                                      See details in the model function help 
     27//              </text> 
     28//[FIXED]=  <text></text> 
     29//[NON_FITTABLE_PARAMS]= <text> lcase_n; Na; Phia; va; La; Nb; Phib; vb; Lb;Nc; Phic; vc; Lc;Nd; Phid; vd; Ld; </text> 
     30//[ORIENTATION_PARAMS]= <text> </text> 
    3031 
    31 typedef struct { 
    32         /// The Case number 
    33         //  [DEFAULT]=lcase_n=0 
    34         int lcase_n; 
     32class RPAModel{ 
     33public: 
     34  // Model parameters 
     35  /// The Case number 
     36  //  [DEFAULT]=lcase_n=0 
     37  Parameter lcase_n; 
    3538 
    36     /// Segment Length ba 
    37     //  [DEFAULT]=ba= 5.0 
    38         double ba; 
    39         /// Segment Length bb 
    40     //  [DEFAULT]=bb=5.0 
    41         double bb; 
    42         /// Segment Length bc 
    43         //  [DEFAULT]=bc= 5.0 
    44         double bc; 
    45         /// Segment Length bd 
    46         //  [DEFAULT]=bd= 5.0 
    47         double bd; 
     39  /// Segment Length ba 
     40  //  [DEFAULT]=ba= 5.0 
     41  Parameter ba; 
     42  /// Segment Length bb 
     43  //  [DEFAULT]=bb=5.0 
     44  Parameter bb; 
     45  /// Segment Length bc 
     46  //  [DEFAULT]=bc= 5.0 
     47  Parameter bc; 
     48  /// Segment Length bd 
     49  //  [DEFAULT]=bd= 5.0 
     50  Parameter bd; 
    4851 
    49         /// Chi Param ab 
    50         //  [DEFAULT]=Kab=-0.0004 
    51         double Kab; 
    52         /// Chi Param ac 
    53         //  [DEFAULT]=Kac=-0.0004 
    54     double Kac; 
    55     /// Chi Param ad 
    56     //  [DEFAULT]=Kad=-0.0004 
    57     double Kad; 
    58     /// Chi Param bc 
    59     //  [DEFAULT]=Kbc=-0.0004 
    60     double Kbc; 
    61     /// Chi Param bd 
    62     //  [DEFAULT]=Kbd=-0.0004 
    63     double Kbd; 
    64     /// Chi Param cd 
    65     //  [DEFAULT]=Kcd=-0.0004 
    66     double Kcd; 
     52  /// Chi Param ab 
     53  //  [DEFAULT]=Kab=-0.0004 
     54  Parameter Kab; 
     55  /// Chi Param ac 
     56  //  [DEFAULT]=Kac=-0.0004 
     57  Parameter Kac; 
     58  /// Chi Param ad 
     59  //  [DEFAULT]=Kad=-0.0004 
     60  Parameter Kad; 
     61  /// Chi Param bc 
     62  //  [DEFAULT]=Kbc=-0.0004 
     63  Parameter Kbc; 
     64  /// Chi Param bd 
     65  //  [DEFAULT]=Kbd=-0.0004 
     66  Parameter Kbd; 
     67  /// Chi Param cd 
     68  //  [DEFAULT]=Kcd=-0.0004 
     69  Parameter Kcd; 
    6770 
    68     /// Scale factor 
    69     //  [DEFAULT]=scale= 1.0 
    70     double scale; 
    71         /// Incoherent Background [1/cm] 
    72         //  [DEFAULT]=background=0 [1/cm] 
    73     double background; 
     71  /// Scale factor 
     72  //  [DEFAULT]=scale= 1.0 
     73  Parameter scale; 
     74  /// Incoherent Background [1/cm] 
     75  //  [DEFAULT]=background=0 [1/cm] 
     76  Parameter background; 
    7477 
    75     /// Degree OF POLYMERIZATION of a 
    76     //  [DEFAULT]=Na=1000.0 
    77     double Na; 
    78     /// Volume Fraction of a 
    79     //  [DEFAULT]=Phia=0.25 
    80     double Phia; 
    81     /// Specific Volume of a 
    82     //  [DEFAULT]=va=100.0 
    83     double va; 
    84     /// Scattering Length of a 
    85     //  [DEFAULT]=La=1.0e-012 
    86     double La; 
     78  /// Degree OF POLYMERIZATION of a 
     79  //  [DEFAULT]=Na=1000.0 
     80  Parameter Na; 
     81  /// Volume Fraction of a 
     82  //  [DEFAULT]=Phia=0.25 
     83  Parameter Phia; 
     84  /// Specific Volume of a 
     85  //  [DEFAULT]=va=100.0 
     86  Parameter va; 
     87  /// Scattering Length of a 
     88  //  [DEFAULT]=La=1.0e-012 
     89  Parameter La; 
    8790 
    88     /// Degree OF POLYMERIZATION of b 
    89     //  [DEFAULT]=Nb=1000.0 
    90     double Nb; 
    91     /// Volume Fraction of b 
    92     //  [DEFAULT]=Phib=0.25 
    93     double Phib; 
    94     /// Specific Volume of b 
    95     //  [DEFAULT]=vb=100.0 
    96     double vb; 
    97     /// Scattering Length of b 
    98     //  [DEFAULT]=Lb=1.0e-012 
    99     double Lb; 
     91  /// Degree OF POLYMERIZATION of b 
     92  //  [DEFAULT]=Nb=1000.0 
     93  Parameter Nb; 
     94  /// Volume Fraction of b 
     95  //  [DEFAULT]=Phib=0.25 
     96  Parameter Phib; 
     97  /// Specific Volume of b 
     98  //  [DEFAULT]=vb=100.0 
     99  Parameter vb; 
     100  /// Scattering Length of b 
     101  //  [DEFAULT]=Lb=1.0e-012 
     102  Parameter Lb; 
    100103 
    101     /// Degree OF POLYMERIZATION of c 
    102     //  [DEFAULT]=Nc=1000.0 
    103     double Nc; 
    104     /// Volume Fraction of c 
    105     //  [DEFAULT]=Phic=0.25 
    106     double Phic; 
    107     /// Specific Volume of c 
    108     //  [DEFAULT]=vc=100.0 
    109     double vc; 
    110     /// Scattering Length of c 
    111     //  [DEFAULT]=Lc=1.0e-012 
    112     double Lc; 
     104  /// Degree OF POLYMERIZATION of c 
     105  //  [DEFAULT]=Nc=1000.0 
     106  Parameter Nc; 
     107  /// Volume Fraction of c 
     108  //  [DEFAULT]=Phic=0.25 
     109  Parameter Phic; 
     110  /// Specific Volume of c 
     111  //  [DEFAULT]=vc=100.0 
     112  Parameter vc; 
     113  /// Scattering Length of c 
     114  //  [DEFAULT]=Lc=1.0e-012 
     115  Parameter Lc; 
    113116 
    114     /// Degree OF POLYMERIZATION of d 
    115     //  [DEFAULT]=Nd=1000.0 
    116     double Nd; 
    117     /// Volume Fraction of d 
    118     //  [DEFAULT]=Phid=0.25 
    119     double Phid; 
    120     /// Specific Volume of d 
    121     //  [DEFAULT]=vd=100.0 
    122     double vd; 
    123     /// Scattering Length of d 
    124     //  [DEFAULT]=Ld=0.0e-012 
    125     double Ld; 
     117  /// Degree OF POLYMERIZATION of d 
     118  //  [DEFAULT]=Nd=1000.0 
     119  Parameter Nd; 
     120  /// Volume Fraction of d 
     121  //  [DEFAULT]=Phid=0.25 
     122  Parameter Phid; 
     123  /// Specific Volume of d 
     124  //  [DEFAULT]=vd=100.0 
     125  Parameter vd; 
     126  /// Scattering Length of d 
     127  //  [DEFAULT]=Ld=0.0e-012 
     128  Parameter Ld; 
    126129 
    127 } RPAParameters; 
     130  // Constructor 
     131  RPAModel(); 
    128132 
    129 double rpa_kernel(double dq[], double q); 
    130  
    131 /// 1D scattering function 
    132 double rpa_analytical_1D(RPAParameters *pars, double q); 
    133  
    134 /// 2D scattering function 
    135 double rpa_analytical_2D(RPAParameters *pars, double q, double phi); 
    136 double rpa_analytical_2DXY(RPAParameters *pars, double qx, double qy); 
    137  
     133  // Operators to get I(Q) 
     134  double operator()(double q); 
     135  double operator()(double qx, double qy); 
     136  double calculate_ER(); 
     137  double evaluate_rphi(double q, double phi); 
     138}; 
    138139#endif 

  • sansmodels/src/c_models/models.hh

    rc637521 r7ffa8196  
    2626 
    2727using namespace std; 
    28  
    29  
    30  
    31  
    32  
    33 class RPAModel{ 
    34 public: 
    35         // Model parameters 
    36         Parameter lcase_n; 
    37         Parameter ba; 
    38         Parameter bb; 
    39         Parameter bc; 
    40         Parameter bd; 
    41  
    42         Parameter Kab; 
    43         Parameter Kac; 
    44         Parameter Kad; 
    45         Parameter Kbc; 
    46         Parameter Kbd; 
    47         Parameter Kcd; 
    48  
    49         Parameter scale; 
    50         Parameter background; 
    51  
    52         Parameter Na; 
    53         Parameter Phia; 
    54         Parameter va; 
    55         Parameter La; 
    56  
    57         Parameter Nb; 
    58         Parameter Phib; 
    59         Parameter vb; 
    60         Parameter Lb; 
    61  
    62         Parameter Nc; 
    63         Parameter Phic; 
    64         Parameter vc; 
    65         Parameter Lc; 
    66  
    67         Parameter Nd; 
    68         Parameter Phid; 
    69         Parameter vd; 
    70         Parameter Ld; 
    71  
    72         // Constructor 
    73         RPAModel(); 
    74  
    75         // Operators to get I(Q) 
    76         double operator()(double q); 
    77         double operator()(double qx, double qy); 
    78         double calculate_ER(); 
    79         double evaluate_rphi(double q, double phi); 
    80 }; 
    81  
    8228 
    8329class ReflModel{ 

  • sansmodels/src/c_models/rpa.cpp

    r67424cd r7ffa8196  
    2121 
    2222#include <math.h> 
    23 #include "models.hh" 
    2423#include "parameters.hh" 
    2524#include <stdio.h> 
     25#include "rpa.h" 
     26 
    2627using namespace std; 
    2728 
    28 extern "C" { 
    29         #include "rpa.h" 
     29static double rpa_kernel(double dp[], double q) { 
     30  int lCASE = dp[0]; 
     31 
     32  double Na,Nb,Nc,Nd,Nab,Nac,Nad,Nbc,Nbd,Ncd; 
     33  double Phia,Phib,Phic,Phid,Phiab,Phiac,Phiad; 
     34  double Phibc,Phibd,Phicd; 
     35  double va,vb,vc,vd,vab,vac,vad,vbc,vbd,vcd; 
     36  double m; 
     37  double ba,bb,bc,bd; 
     38  double Xa,Xb,Xc,Xd; 
     39  double Paa,S0aa,Pab,S0ab,Pac,S0ac,Pad,S0ad; 
     40  double S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd; 
     41  double S0ca,S0cb,Pcc,S0cc,Pcd,S0cd; 
     42  double S0da,S0db,S0dc,Pdd,S0dd; 
     43  double Kaa,Kab,Kac,Kad,Kba,Kbb,Kbc,Kbd; 
     44  double Kca,Kcb,Kcc,Kcd,Kda,Kdb,Kdc,Kdd; 
     45  double Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc; 
     46  double DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33; 
     47  double Y1,Y2,Y3,X11,X12,X13,X21,X22,X23,X31,X32,X33; 
     48  double ZZ,DenQ1,DenQ2,DenQ3,DenQ,Q11,Q12,Q13,Q21,Q22,Q23,Q31,Q32,Q33; 
     49  double N11,N12,N13,N21,N22,N23,N31,N32,N33; 
     50  double M11,M12,M13,M21,M22,M23,M31,M32,M33; 
     51  double S11,S12,S13,S14,S21,S22,S23,S24; 
     52  double S31,S32,S33,S34,S41,S42,S43,S44; 
     53  double La,Lb,Lc,Ld,Lad,Lbd,Lcd,Nav,Intg; 
     54  double scale, background; 
     55 
     56  int ii=13;  //dp[ii<13] = fittable, else fixed 
     57 
     58  Na=1000.0; 
     59  Nb=1000.0; 
     60  Nc=1000.0; 
     61  Nd=1000.0;  //DEGREE OF POLYMERIZATION 
     62  Phia=0.25; 
     63  Phib=0.25; 
     64  Phic=0.25; 
     65  Phid=0.25 ; //VOL FRACTION 
     66  Kab=-0.0004; 
     67  Kac=-0.0004; 
     68  Kad=-0.0004;  //CHI PARAM 
     69  Kbc=-0.0004; 
     70  Kbd=-0.0004; 
     71  Kcd=-0.0004; 
     72  La=1.0e-12; 
     73  Lb=1.0e-12; 
     74  Lc=1.0e-12; 
     75  Ld=0.0e-12;   //SCATT. LENGTH 
     76  va=100.0; 
     77  vb=100.0; 
     78  vc=100.0; 
     79  vd=100.0  ; //SPECIFIC VOLUME 
     80  ba=5.0; 
     81  bb=5.0; 
     82  bc=5.0; 
     83  bd=5.0;   //SEGMENT LENGTH 
     84 
     85  //lCASE was shifted to N-1 from the original code 
     86  if (lCASE <= 1){ 
     87    Phia=0.0000001; 
     88    Phib=0.0000001; 
     89    Phic=0.5; 
     90    Phid=0.5; 
     91    Nc=dp[ii+8]; 
     92    Phic=dp[ii+9]; 
     93    vc=dp[ii+10]; 
     94    Lc=dp[ii+11]; 
     95    bc=dp[3]; 
     96    Nd=dp[ii+12]; 
     97    Phid=dp[ii+13]; 
     98    vd=dp[ii+14]; 
     99    Ld=dp[ii+15]; 
     100    bd=dp[4]; 
     101    Kcd=dp[10]; 
     102    scale=dp[11]; 
     103    background=dp[12]; 
     104  } 
     105  else if ((lCASE > 1) && (lCASE <= 4)){ 
     106    Phia=0.0000001; 
     107    Phib=0.333333; 
     108    Phic=0.333333; 
     109    Phid=0.333333; 
     110    Nb=dp[ii+4]; 
     111    Phib=dp[ii+5]; 
     112    vb=dp[ii+6]; 
     113    Lb=dp[ii+7]; 
     114    bb=dp[2]; 
     115    Nc=dp[ii+8]; 
     116    Phic=dp[ii+9]; 
     117    vc=dp[ii+10]; 
     118    Lc=dp[ii+11]; 
     119    bc=dp[3]; 
     120    Nd=dp[ii+12]; 
     121    Phid=dp[ii+13]; 
     122    vd=dp[ii+14]; 
     123    Ld=dp[ii+15]; 
     124    bd=dp[4]; 
     125    Kbc=dp[8]; 
     126    Kbd=dp[9]; 
     127    Kcd=dp[10]; 
     128    scale=dp[11]; 
     129    background=dp[12]; 
     130  } 
     131  else { 
     132    Phia=0.25; 
     133    Phib=0.25; 
     134    Phic=0.25; 
     135    Phid=0.25; 
     136    Na=dp[ii+0]; 
     137    Phia=dp[ii+1]; 
     138    va=dp[ii+2]; 
     139    La=dp[ii+3]; 
     140    ba=dp[1]; 
     141    Nb=dp[ii+4]; 
     142    Phib=dp[ii+5]; 
     143    vb=dp[ii+6]; 
     144    Lb=dp[ii+7]; 
     145    bb=dp[2]; 
     146    Nc=dp[ii+8]; 
     147    Phic=dp[ii+9]; 
     148    vc=dp[ii+10]; 
     149    Lc=dp[ii+11]; 
     150    bc=dp[3]; 
     151    Nd=dp[ii+12]; 
     152    Phid=dp[ii+13]; 
     153    vd=dp[ii+14]; 
     154    Ld=dp[ii+15]; 
     155    bd=dp[4]; 
     156    Kab=dp[5]; 
     157    Kac=dp[6]; 
     158    Kad=dp[7]; 
     159    Kbc=dp[8]; 
     160    Kbd=dp[9]; 
     161    Kcd=dp[10]; 
     162    scale=dp[11]; 
     163    background=dp[12]; 
     164  } 
     165 
     166  Nab=pow((Na*Nb),(0.5)); 
     167  Nac=pow((Na*Nc),(0.5)); 
     168  Nad=pow((Na*Nd),(0.5)); 
     169  Nbc=pow((Nb*Nc),(0.5)); 
     170  Nbd=pow((Nb*Nd),(0.5)); 
     171  Ncd=pow((Nc*Nd),(0.5)); 
     172 
     173  vab=pow((va*vb),(0.5)); 
     174  vac=pow((va*vc),(0.5)); 
     175  vad=pow((va*vd),(0.5)); 
     176  vbc=pow((vb*vc),(0.5)); 
     177  vbd=pow((vb*vd),(0.5)); 
     178  vcd=pow((vc*vd),(0.5)); 
     179 
     180  Phiab=pow((Phia*Phib),(0.5)); 
     181  Phiac=pow((Phia*Phic),(0.5)); 
     182  Phiad=pow((Phia*Phid),(0.5)); 
     183  Phibc=pow((Phib*Phic),(0.5)); 
     184  Phibd=pow((Phib*Phid),(0.5)); 
     185  Phicd=pow((Phic*Phid),(0.5)); 
     186 
     187  Xa=q*q*ba*ba*Na/6.0; 
     188  Xb=q*q*bb*bb*Nb/6.0; 
     189  Xc=q*q*bc*bb*Nc/6.0; 
     190  Xd=q*q*bd*bb*Nd/6.0; 
     191 
     192  Paa=2.0*(exp(-Xa)-1.0+Xa)/pow(Xa,2.0); 
     193  S0aa=Na*Phia*va*Paa; 
     194  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); 
     195  S0ab=(Phiab*vab*Nab)*Pab; 
     196  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 
     197  S0ac=(Phiac*vac*Nac)*Pac; 
     198  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 
     199  S0ad=(Phiad*vad*Nad)*Pad; 
     200 
     201  S0ba=S0ab; 
     202  Pbb=2.0*(exp(-Xb)-1.0+Xb)/pow(Xb,2.0); 
     203  S0bb=Nb*Phib*vb*Pbb; 
     204  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 
     205  S0bc=(Phibc*vbc*Nbc)*Pbc; 
     206  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 
     207  S0bd=(Phibd*vbd*Nbd)*Pbd; 
     208 
     209  S0ca=S0ac; 
     210  S0cb=S0bc; 
     211  Pcc=2.0*(exp(-Xc)-1.0+Xc)/pow(Xc,2.0); 
     212  S0cc=Nc*Phic*vc*Pcc; 
     213  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 
     214  S0cd=(Phicd*vcd*Ncd)*Pcd; 
     215 
     216  S0da=S0ad; 
     217  S0db=S0bd; 
     218  S0dc=S0cd; 
     219  Pdd=2.0*(exp(-Xd)-1.0+Xd)/pow(Xd,2.0); 
     220  S0dd=Nd*Phid*vd*Pdd; 
     221  //lCASE was shifted to N-1 from the original code 
     222  switch(lCASE){ 
     223  case 0: 
     224    S0aa=0.000001; 
     225    S0ab=0.000002; 
     226    S0ac=0.000003; 
     227    S0ad=0.000004; 
     228    S0bb=0.000005; 
     229    S0bc=0.000006; 
     230    S0bd=0.000007; 
     231    S0cd=0.000008; 
     232    break; 
     233  case 1: 
     234    S0aa=0.000001; 
     235    S0ab=0.000002; 
     236    S0ac=0.000003; 
     237    S0ad=0.000004; 
     238    S0bb=0.000005; 
     239    S0bc=0.000006; 
     240    S0bd=0.000007; 
     241    break; 
     242  case 2: 
     243    S0aa=0.000001; 
     244    S0ab=0.000002; 
     245    S0ac=0.000003; 
     246    S0ad=0.000004; 
     247    S0bc=0.000005; 
     248    S0bd=0.000006; 
     249    S0cd=0.000007; 
     250    break; 
     251  case 3: 
     252    S0aa=0.000001; 
     253    S0ab=0.000002; 
     254    S0ac=0.000003; 
     255    S0ad=0.000004; 
     256    S0bc=0.000005; 
     257    S0bd=0.000006; 
     258    break; 
     259  case 4: 
     260    S0aa=0.000001; 
     261    S0ab=0.000002; 
     262    S0ac=0.000003; 
     263    S0ad=0.000004; 
     264    break; 
     265  case 5: 
     266    S0ab=0.000001; 
     267    S0ac=0.000002; 
     268    S0ad=0.000003; 
     269    S0bc=0.000004; 
     270    S0bd=0.000005; 
     271    S0cd=0.000006; 
     272    break; 
     273  case 6: 
     274    S0ab=0.000001; 
     275    S0ac=0.000002; 
     276    S0ad=0.000003; 
     277    S0bc=0.000004; 
     278    S0bd=0.000005; 
     279    break; 
     280  case 7: 
     281    S0ab=0.000001; 
     282    S0ac=0.000002; 
     283    S0ad=0.000003; 
     284    break; 
     285  case 8: 
     286    S0ac=0.000001; 
     287    S0ad=0.000002; 
     288    S0bc=0.000003; 
     289    S0bd=0.000004; 
     290    break; 
     291  default : //case 9: 
     292    break; 
     293  } 
     294  S0ba=S0ab; 
     295  S0ca=S0ac; 
     296  S0cb=S0bc; 
     297  S0da=S0ad; 
     298  S0db=S0bd; 
     299  S0dc=S0cd; 
     300 
     301  Kaa=0.0; 
     302  Kbb=0.0; 
     303  Kcc=0.0; 
     304  Kdd=0.0; 
     305 
     306  Kba=Kab; 
     307  Kca=Kac; 
     308  Kcb=Kbc; 
     309  Kda=Kad; 
     310  Kdb=Kbd; 
     311  Kdc=Kcd; 
     312 
     313  Zaa=Kaa-Kad-Kad; 
     314  Zab=Kab-Kad-Kbd; 
     315  Zac=Kac-Kad-Kcd; 
     316  Zba=Kba-Kbd-Kad; 
     317  Zbb=Kbb-Kbd-Kbd; 
     318  Zbc=Kbc-Kbd-Kcd; 
     319  Zca=Kca-Kcd-Kad; 
     320  Zcb=Kcb-Kcd-Kbd; 
     321  Zcc=Kcc-Kcd-Kcd; 
     322 
     323  DenT=(-(S0ac*S0bb*S0ca) + S0ab*S0bc*S0ca + S0ac*S0ba*S0cb - S0aa*S0bc*S0cb - S0ab*S0ba*S0cc + S0aa*S0bb*S0cc); 
     324 
     325  T11= (-(S0bc*S0cb) + S0bb*S0cc)/DenT; 
     326  T12= (S0ac*S0cb - S0ab*S0cc)/DenT; 
     327  T13= (-(S0ac*S0bb) + S0ab*S0bc)/DenT; 
     328  T21= (S0bc*S0ca - S0ba*S0cc)/DenT; 
     329  T22= (-(S0ac*S0ca) + S0aa*S0cc)/DenT; 
     330  T23= (S0ac*S0ba - S0aa*S0bc)/DenT; 
     331  T31= (-(S0bb*S0ca) + S0ba*S0cb)/DenT; 
     332  T32= (S0ab*S0ca - S0aa*S0cb)/DenT; 
     333  T33= (-(S0ab*S0ba) + S0aa*S0bb)/DenT; 
     334 
     335  Y1=T11*S0ad+T12*S0bd+T13*S0cd+1.0; 
     336  Y2=T21*S0ad+T22*S0bd+T23*S0cd+1.0; 
     337  Y3=T31*S0ad+T32*S0bd+T33*S0cd+1.0; 
     338 
     339  X11=Y1*Y1; 
     340  X12=Y1*Y2; 
     341  X13=Y1*Y3; 
     342  X21=Y2*Y1; 
     343  X22=Y2*Y2; 
     344  X23=Y2*Y3; 
     345  X31=Y3*Y1; 
     346  X32=Y3*Y2; 
     347  X33=Y3*Y3; 
     348 
     349  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 
     350 
     351  m=1.0/(S0dd-ZZ); 
     352 
     353  N11=m*X11+Zaa; 
     354  N12=m*X12+Zab; 
     355  N13=m*X13+Zac; 
     356  N21=m*X21+Zba; 
     357  N22=m*X22+Zbb; 
     358  N23=m*X23+Zbc; 
     359  N31=m*X31+Zca; 
     360  N32=m*X32+Zcb; 
     361  N33=m*X33+Zcc; 
     362 
     363  M11= N11*S0aa + N12*S0ab + N13*S0ac; 
     364  M12= N11*S0ab + N12*S0bb + N13*S0bc; 
     365  M13= N11*S0ac + N12*S0bc + N13*S0cc; 
     366  M21= N21*S0aa + N22*S0ab + N23*S0ac; 
     367  M22= N21*S0ab + N22*S0bb + N23*S0bc; 
     368  M23= N21*S0ac + N22*S0bc + N23*S0cc; 
     369  M31= N31*S0aa + N32*S0ab + N33*S0ac; 
     370  M32= N31*S0ab + N32*S0bb + N33*S0bc; 
     371  M33= N31*S0ac + N32*S0bc + N33*S0cc; 
     372 
     373  DenQ1=1.0+M11-M12*M21+M22+M11*M22-M13*M31-M13*M22*M31; 
     374  DenQ2=  M12*M23*M31+M13*M21*M32-M23*M32-M11*M23*M32+M33+M11*M33; 
     375  DenQ3=  -M12*M21*M33+M22*M33+M11*M22*M33; 
     376  DenQ=DenQ1+DenQ2+DenQ3; 
     377 
     378  Q11= (1.0 + M22-M23*M32 + M33 + M22*M33)/DenQ; 
     379  Q12= (-M12 + M13*M32 - M12*M33)/DenQ; 
     380  Q13= (-M13 - M13*M22 + M12*M23)/DenQ; 
     381  Q21= (-M21 + M23*M31 - M21*M33)/DenQ; 
     382  Q22= (1.0 + M11 - M13*M31 + M33 + M11*M33)/DenQ; 
     383  Q23= (M13*M21 - M23 - M11*M23)/DenQ; 
     384  Q31= (-M31 - M22*M31 + M21*M32)/DenQ; 
     385  Q32= (M12*M31 - M32 - M11*M32)/DenQ; 
     386  Q33= (1.0 + M11 - M12*M21 + M22 + M11*M22)/DenQ; 
     387 
     388  S11= Q11*S0aa + Q21*S0ab + Q31*S0ac; 
     389  S12= Q12*S0aa + Q22*S0ab + Q32*S0ac; 
     390  S13= Q13*S0aa + Q23*S0ab + Q33*S0ac; 
     391  S14=-S11-S12-S13; 
     392  S21= Q11*S0ba + Q21*S0bb + Q31*S0bc; 
     393  S22= Q12*S0ba + Q22*S0bb + Q32*S0bc; 
     394  S23= Q13*S0ba + Q23*S0bb + Q33*S0bc; 
     395  S24=-S21-S22-S23; 
     396  S31= Q11*S0ca + Q21*S0cb + Q31*S0cc; 
     397  S32= Q12*S0ca + Q22*S0cb + Q32*S0cc; 
     398  S33= Q13*S0ca + Q23*S0cb + Q33*S0cc; 
     399  S34=-S31-S32-S33; 
     400  S41=S14; 
     401  S42=S24; 
     402  S43=S34; 
     403  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 
     404 
     405  Nav=6.022045e+23; 
     406  Lad=(La/va-Ld/vd)*sqrt(Nav); 
     407  Lbd=(Lb/vb-Ld/vd)*sqrt(Nav); 
     408  Lcd=(Lc/vc-Ld/vd)*sqrt(Nav); 
     409 
     410  Intg=pow(Lad,2.0)*S11+pow(Lbd,2.0)*S22+pow(Lcd,2.0)*S33+2.0*Lad*Lbd*S12+2.0*Lbd*Lcd*S23+2.0*Lad*Lcd*S13; 
     411 
     412  Intg *= scale; 
     413  Intg += background; 
     414 
     415  return Intg; 
     416 
    30417} 
    31418 
    32419RPAModel :: RPAModel() { 
    33         lcase_n = Parameter(0.0); 
    34         ba = Parameter(5.0); 
    35         bb = Parameter(5.0); 
    36         bc = Parameter(5.0); 
    37         bd = Parameter(5.0); 
    38  
    39     Kab = Parameter(-0.0004); 
    40     Kac = Parameter(-0.0004); 
    41     Kad = Parameter(-0.0004); 
    42     Kbc = Parameter(-0.0004); 
    43     Kbd = Parameter(-0.0004); 
    44     Kcd = Parameter(-0.0004); 
    45  
    46     scale = Parameter(1.0); 
    47     background = Parameter(0.0); 
    48  
    49     Na = Parameter(1000.0); 
    50     Phia = Parameter(0.25); 
    51     va = Parameter(100.0); 
    52     La = Parameter(1.0e-012); 
    53  
    54     Nb = Parameter(1000.0); 
    55     Phib = Parameter(0.25); 
    56     vb = Parameter(100.0); 
    57     Lb = Parameter(1.0e-012); 
    58  
    59     Nc = Parameter(1000.0); 
    60     Phic = Parameter(0.25); 
    61     vc = Parameter(100.0); 
    62     Lc = Parameter(1.0e-012); 
    63  
    64     Nd = Parameter(1000.0); 
    65     Phid = Parameter(0.25); 
    66     vd = Parameter(100.0); 
    67     Ld = Parameter(0.0e-012); 
     420  lcase_n = Parameter(0.0); 
     421  ba = Parameter(5.0); 
     422  bb = Parameter(5.0); 
     423  bc = Parameter(5.0); 
     424  bd = Parameter(5.0); 
     425 
     426  Kab = Parameter(-0.0004); 
     427  Kac = Parameter(-0.0004); 
     428  Kad = Parameter(-0.0004); 
     429  Kbc = Parameter(-0.0004); 
     430  Kbd = Parameter(-0.0004); 
     431  Kcd = Parameter(-0.0004); 
     432 
     433  scale = Parameter(1.0); 
     434  background = Parameter(0.0); 
     435 
     436  Na = Parameter(1000.0); 
     437  Phia = Parameter(0.25); 
     438  va = Parameter(100.0); 
     439  La = Parameter(1.0e-012); 
     440 
     441  Nb = Parameter(1000.0); 
     442  Phib = Parameter(0.25); 
     443  vb = Parameter(100.0); 
     444  Lb = Parameter(1.0e-012); 
     445 
     446  Nc = Parameter(1000.0); 
     447  Phic = Parameter(0.25); 
     448  vc = Parameter(100.0); 
     449  Lc = Parameter(1.0e-012); 
     450 
     451  Nd = Parameter(1000.0); 
     452  Phid = Parameter(0.25); 
     453  vd = Parameter(100.0); 
     454  Ld = Parameter(0.0e-012); 
    68455 
    69456} 
    … …  
    76463 */ 
    77464double RPAModel :: operator()(double q) { 
    78         double dp[29]; 
    79         // Fill parameter array for IGOR library 
    80         // Add the background after averaging 
    81         //Fittable parameters 
    82         dp[0] = lcase_n(); 
    83  
    84         dp[1] = ba(); 
    85         dp[2] = bb(); 
    86         dp[3] = bc(); 
    87         dp[4] = bd(); 
    88  
    89         dp[5] = Kab(); 
    90         dp[6] = Kac(); 
    91         dp[7] = Kad(); 
    92         dp[8] = Kbc(); 
    93         dp[9] = Kbd(); 
    94         dp[10] = Kcd(); 
    95  
    96         dp[11] = scale(); 
    97         dp[12] = background(); 
    98  
    99         //Fixed parameters 
    100         dp[13] = Na(); 
    101         dp[14] = Phia(); 
    102         dp[15] = va(); 
    103         dp[16] = La(); 
    104  
    105         dp[17] = Nb(); 
    106         dp[18] = Phib(); 
    107         dp[19] = vb(); 
    108         dp[20] = Lb(); 
    109  
    110         dp[21] = Nc(); 
    111         dp[22] = Phic(); 
    112         dp[23] = vc(); 
    113         dp[24] = Lc(); 
    114  
    115         dp[25] = Nd(); 
    116         dp[26] = Phid(); 
    117         dp[27] = vd(); 
    118         dp[28] = Ld(); 
    119  
    120         return rpa_kernel(dp,q); 
     465  double dp[29]; 
     466  // Fill parameter array for IGOR library 
     467  // Add the background after averaging 
     468  //Fittable parameters 
     469  dp[0] = lcase_n(); 
     470 
     471  dp[1] = ba(); 
     472  dp[2] = bb(); 
     473  dp[3] = bc(); 
     474  dp[4] = bd(); 
     475 
     476  dp[5] = Kab(); 
     477  dp[6] = Kac(); 
     478  dp[7] = Kad(); 
     479  dp[8] = Kbc(); 
     480  dp[9] = Kbd(); 
     481  dp[10] = Kcd(); 
     482 
     483  dp[11] = scale(); 
     484  dp[12] = background(); 
     485 
     486  //Fixed parameters 
     487  dp[13] = Na(); 
     488  dp[14] = Phia(); 
     489  dp[15] = va(); 
     490  dp[16] = La(); 
     491 
     492  dp[17] = Nb(); 
     493  dp[18] = Phib(); 
     494  dp[19] = vb(); 
     495  dp[20] = Lb(); 
     496 
     497  dp[21] = Nc(); 
     498  dp[22] = Phic(); 
     499  dp[23] = vc(); 
     500  dp[24] = Lc(); 
     501 
     502  dp[25] = Nd(); 
     503  dp[26] = Phid(); 
     504  dp[27] = vd(); 
     505  dp[28] = Ld(); 
     506 
     507  return rpa_kernel(dp,q); 
    121508} 
    122509 
    … …  
    128515 */ 
    129516double RPAModel :: operator()(double qx, double qy) { 
    130         double q = sqrt(qx*qx + qy*qy); 
    131         return (*this).operator()(q); 
     517  double q = sqrt(qx*qx + qy*qy); 
     518  return (*this).operator()(q); 
    132519} 
    133520 
    … …  
    140527 */ 
    141528double RPAModel :: evaluate_rphi(double q, double phi) { 
    142         return (*this).operator()(q); 
     529  return (*this).operator()(q); 
    143530} 
    144531 
    … …  
    148535 */ 
    149536double RPAModel :: calculate_ER() { 
    150 //NOT implemented!!! 
    151         return 0.0; 
     537  //NOT implemented!!! 
     538  return 0.0; 
    152539} 

  • sansmodels/src/python_wrapper/CRPAModel.cpp

    r67424cd r7ffa8196  
    1818 * 
    1919 * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY 
    20  *          DO NOT MODIFY THIS FILE, MODIFY rpa.h 
     20 *          DO NOT MODIFY THIS FILE, MODIFY ../c_extensions/rpa.h 
    2121 *          AND RE-RUN THE GENERATOR SCRIPT 
    2222 * 
    … …  
    3333#include <math.h> 
    3434#include <time.h> 
     35 
     36} 
     37 
    3538#include "rpa.h" 
    36 } 
    37  
    38 #include "models.hh" 
    3939#include "dispersion_visitor.hh" 
    4040 
Note: See TracChangeset for help on using the changeset viewer.