Changeset 7ffa8196 in sasview
- Timestamp:
- Jan 5, 2012 10:25:27 AM (13 years ago)
- 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
- Location:
- sansmodels/src
- Files:
-
- 1 deleted
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/c_extensions/rpa.h
r67424cd r7ffa8196 1 1 #if !defined(o_h) 2 2 #define rpa_h 3 #include "parameters.hh" 3 4 4 5 /** 5 6 * Structure definition for sphere parameters 6 7 */ 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 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> 30 31 31 typedef struct { 32 /// The Case number 33 // [DEFAULT]=lcase_n=0 34 int lcase_n; 32 class RPAModel{ 33 public: 34 // Model parameters 35 /// The Case number 36 // [DEFAULT]=lcase_n=0 37 Parameter lcase_n; 35 38 36 37 38 doubleba;39 40 41 doublebb;42 43 44 doublebc;45 46 47 doublebd;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; 48 51 49 50 51 doubleKab;52 53 54 doubleKac;55 56 57 doubleKad;58 59 60 doubleKbc;61 62 63 doubleKbd;64 65 66 doubleKcd;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; 67 70 68 69 70 doublescale;71 72 73 doublebackground;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; 74 77 75 76 77 doubleNa;78 79 80 doublePhia;81 82 83 doubleva;84 85 86 doubleLa;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; 87 90 88 89 90 doubleNb;91 92 93 doublePhib;94 95 96 doublevb;97 98 99 doubleLb;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; 100 103 101 102 103 doubleNc;104 105 106 doublePhic;107 108 109 doublevc;110 111 112 doubleLc;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; 113 116 114 115 116 doubleNd;117 118 119 doublePhid;120 121 122 doublevd;123 124 125 doubleLd;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; 126 129 127 } RPAParameters; 130 // Constructor 131 RPAModel(); 128 132 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 }; 138 139 #endif -
sansmodels/src/c_models/models.hh
rc637521 r7ffa8196 26 26 27 27 using namespace std; 28 29 30 31 32 33 class RPAModel{34 public:35 // Model parameters36 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 // Constructor73 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 82 28 83 29 class ReflModel{ -
sansmodels/src/c_models/rpa.cpp
r67424cd r7ffa8196 21 21 22 22 #include <math.h> 23 #include "models.hh"24 23 #include "parameters.hh" 25 24 #include <stdio.h> 25 #include "rpa.h" 26 26 27 using namespace std; 27 28 28 extern "C" { 29 #include "rpa.h" 29 static 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 30 417 } 31 418 32 419 RPAModel :: RPAModel() { 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 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); 68 455 69 456 } … … 76 463 */ 77 464 double RPAModel :: operator()(double q) { 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 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); 121 508 } 122 509 … … 128 515 */ 129 516 double RPAModel :: operator()(double qx, double qy) { 130 131 517 double q = sqrt(qx*qx + qy*qy); 518 return (*this).operator()(q); 132 519 } 133 520 … … 140 527 */ 141 528 double RPAModel :: evaluate_rphi(double q, double phi) { 142 529 return (*this).operator()(q); 143 530 } 144 531 … … 148 535 */ 149 536 double RPAModel :: calculate_ER() { 150 //NOT implemented!!!151 537 //NOT implemented!!! 538 return 0.0; 152 539 } -
sansmodels/src/python_wrapper/CRPAModel.cpp
r67424cd r7ffa8196 18 18 * 19 19 * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY 20 * DO NOT MODIFY THIS FILE, MODIFY rpa.h20 * DO NOT MODIFY THIS FILE, MODIFY ../c_extensions/rpa.h 21 21 * AND RE-RUN THE GENERATOR SCRIPT 22 22 * … … 33 33 #include <math.h> 34 34 #include <time.h> 35 36 } 37 35 38 #include "rpa.h" 36 }37 38 #include "models.hh"39 39 #include "dispersion_visitor.hh" 40 40
Note: See TracChangeset
for help on using the changeset viewer.