source: sasmodels/sasmodels/models/core_shell_parallelepiped.c @ 3a1fc7d

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 3a1fc7d was 3a1fc7d, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

core_shell_parallelepiped: rewrite so that integrated 2D matches 1D

  • Property mode set to 100644
File size: 5.9 KB
RevLine 
[becded3]1static double
2form_volume(double length_a, double length_b, double length_c,
3    double thick_rim_a, double thick_rim_b, double thick_rim_c)
[44bd2be]4{
[2a0b2b1]5    return length_a * length_b * length_c +
6           2.0 * thick_rim_a * length_b * length_c +
[2222134]7           2.0 * thick_rim_b * length_a * length_c +
8           2.0 * thick_rim_c * length_a * length_b;
[44bd2be]9}
10
[becded3]11static double
12Iq(double q,
[44bd2be]13    double core_sld,
14    double arim_sld,
15    double brim_sld,
16    double crim_sld,
17    double solvent_sld,
[2222134]18    double length_a,
19    double length_b,
20    double length_c,
21    double thick_rim_a,
22    double thick_rim_b,
23    double thick_rim_c)
[44bd2be]24{
[3a1fc7d]25    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
[44bd2be]26    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
[c69d6d6]27    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
[2a0b2b1]28
[14838a3]29    const double mu = 0.5 * q * length_b;
[2a0b2b1]30
[44bd2be]31    // Scale sides by B
[3a1fc7d]32    const double a_over_b = length_a / length_b;
33    const double c_over_b = length_c / length_b;
[14838a3]34
[3a1fc7d]35    double tA_over_b = a_over_b + 2.0*thick_rim_a/length_b;
36    double tB_over_b = 1+ 2.0*thick_rim_b/length_b;
37    double tC_over_b = c_over_b + 2.0*thick_rim_c/length_b;
[14838a3]38
39    double Vin = length_a * length_b * length_c;
[3a1fc7d]40    double VtA = (2.0 * thick_rim_a * length_b * length_c);
41    double VtB = (2.0 * length_a * thick_rim_b * length_c);
42    double VtC = (2.0 * length_a * length_b * thick_rim_c);
[14838a3]43
44    // Scale factors (note that drC is not used later)
[3a1fc7d]45    const double dr0 = (core_sld-solvent_sld);
46    const double drA = (arim_sld-solvent_sld);
47    const double drB = (brim_sld-solvent_sld);
48    const double drC = (crim_sld-solvent_sld);
[14838a3]49
50    // Precompute scale factors for combining cross terms from the shape
[3a1fc7d]51    const double dr0_Vin = dr0*Vin;
52    const double drA_VtA = drA*VtA;
53    const double drB_VtB = drB*VtB;
54    const double drC_VtC = drC*VtC;
55    const double drV_delta = dr0_Vin - drA_VtA - drB_VtB - drC_VtC;
56
57    /*  *************** algorithm description ******************
58
59    // Rewrite f as x*siC + y*siCt to move the siC/siCt calculation out
60    // of the inner loop.  That is:
61
62    f = (di-da-db-dc) sa sb sc + da sa' sb sc + db sa sb' sc + dc sa sb sc'
63      =  [ (di-da-db-dc) sa sb + da sa' sb + db sa sb' ] sc  + [dc sa sb] sc'
64      = x sc + y sc'
65
66    // where:
67    di = delta rho_core V_core
68    da = delta rho_rimA V_rimA
69    db = delta rho_rimB V_rimB
70    dc = delta rho_rimC V_rimC
71    sa = j_0 (q_a a/2)    // siA the code
72    sb = j_0 (q_b b/2)
73    sc = j_0 (q_c c/2)
74    sa' = j_0(q_a a_rim/2)  // siAt the code
75    sb' = j_0(q_b b_rim/2)
76    sc' = j_0(q_c c_rim/2)
77
78    // qa, qb, and qc are generated using polar coordinates, with the
79    // outer loop integrating over [0,1] after the u-substitution
80    //    sigma = cos(theta), sqrt(1-sigma^2) = sin(theta)
81    // and inner loop integrating over [0, pi/2] as
82    //    uu = phi
83
84    ************************************************************  */
[14838a3]85
86    // outer integral (with gauss points), integration limits = 0, 1
[3a1fc7d]87    double outer_sum = 0; //initialize integral
[14838a3]88    for( int i=0; i<76; i++) {
89        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
90        double mu_proj = mu * sqrt(1.0-sigma*sigma);
91
[3a1fc7d]92        // inner integral (with gauss points), integration limits = 0, pi/2
93        const double siC = sas_sinx_x(mu * sigma * c_over_b);
94        const double siCt = sas_sinx_x(mu * sigma * tC_over_b);
95        double inner_sum = 0.0;
[14838a3]96        for(int j=0; j<76; j++) {
97            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
98            double sin_uu, cos_uu;
99            SINCOS(M_PI_2*uu, sin_uu, cos_uu);
[3a1fc7d]100            const double siA = sas_sinx_x(mu_proj * sin_uu * a_over_b);
101            const double siB = sas_sinx_x(mu_proj * cos_uu );
102            const double siAt = sas_sinx_x(mu_proj * sin_uu * tA_over_b);
103            const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b);
104
105            const double x = drV_delta*siA*siB + drA_VtA*siB*siAt + drB_VtB*siA*siBt;
106            const double form = x*siC + drC_VtC*siA*siB*siCt;
107
108            inner_sum += Gauss76Wt[j] * form * form;
[44bd2be]109        }
[3a1fc7d]110        inner_sum *= 0.5;
[14838a3]111        // now sum up the outer integral
[3a1fc7d]112        outer_sum += Gauss76Wt[i] * inner_sum;
[44bd2be]113    }
[3a1fc7d]114    outer_sum *= 0.5;
[44bd2be]115
[14838a3]116    //convert from [1e-12 A-1] to [cm-1]
[3a1fc7d]117    return 1.0e-4 * outer_sum;
[44bd2be]118}
119
[becded3]120static double
121Iqxy(double qa, double qb, double qc,
[44bd2be]122    double core_sld,
123    double arim_sld,
124    double brim_sld,
125    double crim_sld,
126    double solvent_sld,
[2222134]127    double length_a,
128    double length_b,
129    double length_c,
130    double thick_rim_a,
131    double thick_rim_b,
[becded3]132    double thick_rim_c)
[44bd2be]133{
[14838a3]134    // cspkernel in csparallelepiped recoded here
135    const double dr0 = core_sld-solvent_sld;
136    const double drA = arim_sld-solvent_sld;
137    const double drB = brim_sld-solvent_sld;
138    const double drC = crim_sld-solvent_sld;
139
140    double Vin = length_a * length_b * length_c;
[3a1fc7d]141    double VtA = 2.0 * thick_rim_a * length_b * length_c;
142    double VtB = 2.0 * length_a * thick_rim_b * length_c;
143    double VtC = 2.0 * length_a * length_b * thick_rim_c;
[14838a3]144
[44bd2be]145    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
[c69d6d6]146    // the scaling by B.
[3a1fc7d]147    double tA = length_a + 2.0*thick_rim_a;
148    double tB = length_b + 2.0*thick_rim_b;
149    double tC = length_c + 2.0*thick_rim_c;
[44bd2be]150    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
[2a0b2b1]151    double siA = sas_sinx_x(0.5*length_a*qa);
152    double siB = sas_sinx_x(0.5*length_b*qb);
153    double siC = sas_sinx_x(0.5*length_c*qc);
[3a1fc7d]154    double siAt = sas_sinx_x(0.5*tA*qa);
155    double siBt = sas_sinx_x(0.5*tB*qb);
156    double siCt = sas_sinx_x(0.5*tC*qc);
[2a0b2b1]157
[44bd2be]158
159    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
160    // in the 1D code, but should be checked!
[3a1fc7d]161    double f = ( dr0*Vin*siA*siB*siC
162               + drA*VtA*(siAt-siA)*siB*siC
163               + drB*VtB*siA*(siBt-siB)*siC
164               + drC*VtC*siA*siB*(siCt-siC));
[2a0b2b1]165
[44bd2be]166    return 1.0e-4 * f * f;
167}
Note: See TracBrowser for help on using the repository browser.