source: sasmodels/sasmodels/models/core_shell_parallelepiped.c @ 58c3367

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 58c3367 was abdd01c, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

floating point constants should include decimal point

  • Property mode set to 100644
File size: 10.2 KB
Line 
1double form_volume(double a_side, double b_side, double c_side, 
2                   double arim_thickness, double brim_thickness, double crim_thickness);
3double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld,
4          double solvent_sld, double a_side, double b_side, double c_side,
5          double arim_thickness, double brim_thickness, double crim_thickness);
6double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld,
7            double crim_sld, double solvent_sld, double a_side, double b_side,
8            double c_side, double arim_thickness, double brim_thickness,
9            double crim_thickness, double theta, double phi, double psi);
10
11double form_volume(double a_side, double b_side, double c_side, 
12                   double arim_thickness, double brim_thickness, double crim_thickness)
13{
14    //return a_side * b_side * c_side;
15    return a_side * b_side * c_side + 
16           2.0 * arim_thickness * b_side * c_side + 
17           2.0 * brim_thickness * a_side * c_side +
18           2.0 * crim_thickness * a_side * b_side;
19}
20
21double Iq(double q,
22    double core_sld,
23    double arim_sld,
24    double brim_sld,
25    double crim_sld,
26    double solvent_sld,
27    double a_side,
28    double b_side,
29    double c_side,
30    double arim_thickness,
31    double brim_thickness,
32    double crim_thickness)
33{
34    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
35    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
36   
37    double t1, t2, t3, t4, tmp, answer;   
38    double mu = q * b_side;
39   
40    //calculate volume before rescaling (in original code, but not used)
41    //double vol = form_volume(a_side, b_side, c_side, arim_thickness, brim_thickness, crim_thickness);         
42    //double vol = a_side * b_side * c_side +
43    //       2.0 * arim_thickness * b_side * c_side +
44    //       2.0 * brim_thickness * a_side * c_side +
45    //       2.0 * crim_thickness * a_side * b_side;
46   
47    // Scale sides by B
48    double a_scaled = a_side / b_side;
49    double c_scaled = c_side / b_side;
50
51    // DelRho values (note that drC is not used later)       
52        double dr0 = core_sld-solvent_sld;
53        double drA = arim_sld-solvent_sld;
54        double drB = brim_sld-solvent_sld;
55        //double drC = crim_sld-solvent_sld;
56
57    //Order of integration
58    int nordi=76;                               
59    int nordj=76;
60   
61    // outer integral (with nordi gauss points), integration limits = 0, 1
62    double summ = 0; //initialize integral
63
64    for( int i=0; i<nordi; i++) {
65               
66        // inner integral (with nordj gauss points), integration limits = 0, 1
67       
68        double summj = 0.0;
69            double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );         
70               
71            for(int j=0; j<nordj; j++) {
72
73            double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
74            double mudum = mu * sqrt(1.0-sigma*sigma);
75
76                double Vin = a_side * b_side * c_side;
77                //double Vot = (a_side * b_side * c_side +
78            //            2.0 * arim_thickness * b_side * c_side +
79            //            2.0 * a_side * brim_thickness * c_side +
80            //            2.0 * a_side * b_side * crim_thickness);
81                double V1 = (2.0 * arim_thickness * b_side * c_side);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
82                double V2 = (2.0 * a_side * brim_thickness * c_side);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
83       
84            // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG)
85            // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c,
86            // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions
87            // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!)           
88            double ta = (a_scaled+2.0*arim_thickness)/b_side; 
89            double tb = (a_scaled+2.0*brim_thickness)/b_side;
90   
91                double arg1 = (0.5*mudum*a_scaled) * sin(0.5*M_PI*uu);
92                double arg2 = (0.5*mudum) * cos(0.5*M_PI*uu);
93                double arg3=  (0.5*mudum*ta) * sin(0.5*M_PI*uu);
94                double arg4=  (0.5*mudum*tb) * cos(0.5*M_PI*uu);
95
96                if(arg1==0.0){
97                        t1 = 1.0;
98                } else {
99                        t1 = (sin(arg1)/arg1);                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)   
100                }
101                if(arg2==0.0){
102                        t2 = 1.0;
103                } else {
104                        t2 = (sin(arg2)/arg2);           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)   
105                }       
106                if(arg3==0.0){
107                        t3 = 1.0;
108                } else {
109                        t3 = sin(arg3)/arg3;
110                }
111                if(arg4==0.0){
112                        t4 = 1.0;
113                } else {
114                        t4 = sin(arg4)/arg4;
115                }
116           
117            // Expression in libCylinder.c (neither drC nor Vot are used)
118                tmp =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 );   //  correct FF : square of sum of phase factors           
119           
120            // To note also that in csparallelepiped.cpp, there is a function called
121            // cspkernel, which seems to make more sense and has the following comment:
122            //   This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010.
123            //   tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)*
124            //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
125            // This is the function called by csparallelepiped_analytical_2D_scaled,
126            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c       
127           
128            summj += Gauss76Wt[j] * tmp;
129
130        }
131               
132        // value of the inner integral
133        answer = 0.5 * summj;
134
135                // finish the outer integral
136                double arg = 0.5 * mu* c_scaled *sigma;
137                if ( arg == 0.0 ) {
138                        answer *= 1.0;
139                } else {
140                answer *= sin(arg)*sin(arg)/arg/arg;
141                }
142               
143                // now sum up the outer integral
144                summ += Gauss76Wt[i] * answer;
145
146    }
147   
148        answer = 0.5 * summ;
149
150        //convert from [1e-12 A-1] to [cm-1]
151        answer *= 1.0e-4;
152       
153        return answer;
154}
155
156double Iqxy(double qx, double qy,
157    double core_sld,
158    double arim_sld,
159    double brim_sld,
160    double crim_sld,
161    double solvent_sld,
162    double a_side,
163    double b_side,
164    double c_side,
165    double arim_thickness,
166    double brim_thickness,
167    double crim_thickness,
168    double theta,
169    double phi,
170    double psi)
171{
172    double tmp1, tmp2, tmp3, tmpt1, tmpt2, tmpt3;   
173
174    double q = sqrt(qx*qx+qy*qy);
175    double qx_scaled = qx/q;
176    double qy_scaled = qy/q;
177
178    // Convert angles given in degrees to radians
179    theta *= M_PI_180;
180    phi   *= M_PI_180;
181    psi   *= M_PI_180;
182   
183    // Parallelepiped c axis orientation
184    double cparallel_x = cos(theta) * cos(phi);
185    double cparallel_y = sin(theta);
186   
187    // Compute angle between q and parallelepiped axis
188    double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz;
189
190    // Parallelepiped a axis orientation
191    double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
192    double parallel_y = sin(psi)*cos(theta);
193    double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled;
194
195    // Parallelepiped b axis orientation
196    double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
197    double bparallel_y = cos(theta)*cos(psi);
198    double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled;
199
200    // The following tests should always pass
201    if (fabs(cos_val_c)>1.0) {
202      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
203      cos_val_c = 1.0;
204    }
205    if (fabs(cos_val_a)>1.0) {
206      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
207      cos_val_a = 1.0;
208    }
209    if (fabs(cos_val_b)>1.0) {
210      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
211      cos_val_b = 1.0;
212    }
213   
214    // cspkernel in csparallelepiped recoded here
215    double dr0 = core_sld-solvent_sld;
216        double drA = arim_sld-solvent_sld;
217        double drB = brim_sld-solvent_sld;
218        double drC = crim_sld-solvent_sld;
219        double Vin = a_side * b_side * c_side;
220    // As for the 1D case, Vot is not used
221        //double Vot = (a_side * b_side * c_side +
222    //              2.0 * arim_thickness * b_side * c_side +
223    //              2.0 * a_side * brim_thickness * c_side +
224    //              2.0 * a_side * b_side * crim_thickness);
225        double V1 = (2.0 * arim_thickness * b_side * c_side);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
226        double V2 = (2.0 * a_side * brim_thickness * c_side);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
227    double V3 = (2.0 * a_side * b_side * crim_thickness);
228    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
229    // the scaling by B. The use of a_side for the 3 of them seems clearly a mistake to me,
230    // but for the moment I let it like this until understanding better the code.
231        double ta = a_side + 2.0*arim_thickness;
232    double tb = a_side + 2.0*brim_thickness;
233    double tc = a_side + 2.0*crim_thickness;
234    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
235    double argA = 0.5*q*a_side*cos_val_a;
236    double argB = 0.5*q*b_side*cos_val_b;
237    double argC = 0.5*q*c_side*cos_val_c;
238    double argtA = 0.5*q*ta*cos_val_a;
239    double argtB = 0.5*q*tb*cos_val_b;
240    double argtC = 0.5*q*tc*cos_val_c;
241   
242    if(argA==0.0) {
243        tmp1 = 1.0;
244    } else {
245        tmp1 = sin(argA)/argA;
246    }
247    if (argB==0.0) {
248        tmp2 = 1.0;
249    } else {
250        tmp2 = sin(argB)/argB;
251    }
252    if (argC==0.0) {
253        tmp3 = 1.0;
254    } else {
255        tmp3 = sin(argC)/argC;
256    }
257    if(argtA==0.0) {
258        tmpt1 = 1.0;
259    } else {
260        tmpt1 = sin(argtA)/argtA;
261    }
262    if (argtB==0.0) {
263        tmpt2 = 1.0;
264    } else {
265        tmpt2 = sin(argtB)/argtB;
266    }
267    if (argtC==0.0) {
268        tmpt3 = 1.0;
269    } else {
270        tmpt3 = sin(argtC)*sin(argtC)/argtC/argtC;
271    }
272
273    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
274    // in the 1D code, but should be checked!
275    double f = ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1 + 
276               drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);
277   
278    return 1.0e-4 * f * f;
279}
Note: See TracBrowser for help on using the repository browser.