source: sasmodels/sasmodels/models/core_shell_parallelepiped.c @ e077231

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

correct internal comments in core shell parallelepiped

  • Property mode set to 100644
File size: 5.1 KB
RevLine 
[fc0b7aa]1// Set OVERLAPPING to 1 in order to fill in the edges of the box, with
2// c endcaps and b overlapping a.  With the proper choice of parameters,
3// (setting rim slds to sld, core sld to solvent, rim thickness to thickness
4// and subtracting 2*thickness from length, this should match the hollow
5// rectangular prism.)  Set it to 0 for the documented behaviour.
6#define OVERLAPPING 0
[becded3]7static double
8form_volume(double length_a, double length_b, double length_c,
9    double thick_rim_a, double thick_rim_b, double thick_rim_c)
[44bd2be]10{
[fc0b7aa]11    return
12#if OVERLAPPING
13        // Hollow rectangular prism only includes the volume of the shell
14        // so uncomment the next line when comparing.  Solid rectangular
15        // prism, or parallelepiped want filled cores, so comment when
16        // comparing.
17        //-length_a * length_b * length_c +
18        (length_a + 2.0*thick_rim_a) *
19        (length_b + 2.0*thick_rim_b) *
20        (length_c + 2.0*thick_rim_c);
21#else
22        length_a * length_b * length_c +
23        2.0 * thick_rim_a * length_b * length_c +
24        2.0 * length_a * thick_rim_b * length_c +
25        2.0 * length_a * length_b * thick_rim_c;
26#endif
[44bd2be]27}
28
[becded3]29static double
30Iq(double q,
[44bd2be]31    double core_sld,
32    double arim_sld,
33    double brim_sld,
34    double crim_sld,
35    double solvent_sld,
[2222134]36    double length_a,
37    double length_b,
38    double length_c,
39    double thick_rim_a,
40    double thick_rim_b,
41    double thick_rim_c)
[44bd2be]42{
[3a1fc7d]43    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
[44bd2be]44    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
[e077231]45    // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker)
46    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK)
[2a0b2b1]47
[4493288]48    const double half_q = 0.5*q;
[2a0b2b1]49
[4493288]50    const double tA = length_a + 2.0*thick_rim_a;
51    const double tB = length_b + 2.0*thick_rim_b;
52    const double tC = length_c + 2.0*thick_rim_c;
[14838a3]53
[4493288]54    // Scale factors
[3a1fc7d]55    const double dr0 = (core_sld-solvent_sld);
56    const double drA = (arim_sld-solvent_sld);
57    const double drB = (brim_sld-solvent_sld);
58    const double drC = (crim_sld-solvent_sld);
[14838a3]59
60    // outer integral (with gauss points), integration limits = 0, 1
[3a1fc7d]61    double outer_sum = 0; //initialize integral
[74768cb]62    for( int i=0; i<GAUSS_N; i++) {
[a261a83]63        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 );
[4493288]64        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha);
[14838a3]65
[3a1fc7d]66        // inner integral (with gauss points), integration limits = 0, pi/2
[4493288]67        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q);
68        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q);
[3a1fc7d]69        double inner_sum = 0.0;
[74768cb]70        for(int j=0; j<GAUSS_N; j++) {
[a261a83]71            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 );
[4493288]72            double sin_beta, cos_beta;
73            SINCOS(M_PI_2*beta, sin_beta, cos_beta);
74            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta);
75            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta);
76            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta);
77            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta);
[3a1fc7d]78
[fc0b7aa]79#if OVERLAPPING
[4493288]80            const double f = dr0*siA*siB*siC
81                + drA*(siAt-siA)*siB*siC
82                + drB*siAt*(siBt-siB)*siC
83                + drC*siAt*siBt*(siCt-siC);
[fc0b7aa]84#else
[4493288]85            const double f = dr0*siA*siB*siC
86                + drA*(siAt-siA)*siB*siC
87                + drB*siA*(siBt-siB)*siC
88                + drC*siA*siB*(siCt-siC);
[fc0b7aa]89#endif
90
[74768cb]91            inner_sum += GAUSS_W[j] * f * f;
[44bd2be]92        }
[3a1fc7d]93        inner_sum *= 0.5;
[14838a3]94        // now sum up the outer integral
[74768cb]95        outer_sum += GAUSS_W[i] * inner_sum;
[44bd2be]96    }
[3a1fc7d]97    outer_sum *= 0.5;
[44bd2be]98
[14838a3]99    //convert from [1e-12 A-1] to [cm-1]
[3a1fc7d]100    return 1.0e-4 * outer_sum;
[44bd2be]101}
102
[becded3]103static double
[108e70e]104Iqabc(double qa, double qb, double qc,
[44bd2be]105    double core_sld,
106    double arim_sld,
107    double brim_sld,
108    double crim_sld,
109    double solvent_sld,
[2222134]110    double length_a,
111    double length_b,
112    double length_c,
113    double thick_rim_a,
114    double thick_rim_b,
[becded3]115    double thick_rim_c)
[44bd2be]116{
[14838a3]117    // cspkernel in csparallelepiped recoded here
118    const double dr0 = core_sld-solvent_sld;
119    const double drA = arim_sld-solvent_sld;
120    const double drB = brim_sld-solvent_sld;
121    const double drC = crim_sld-solvent_sld;
122
[fc0b7aa]123    const double tA = length_a + 2.0*thick_rim_a;
124    const double tB = length_b + 2.0*thick_rim_b;
125    const double tC = length_c + 2.0*thick_rim_c;
[4493288]126    const double siA = length_a*sas_sinx_x(0.5*length_a*qa);
127    const double siB = length_b*sas_sinx_x(0.5*length_b*qb);
128    const double siC = length_c*sas_sinx_x(0.5*length_c*qc);
129    const double siAt = tA*sas_sinx_x(0.5*tA*qa);
130    const double siBt = tB*sas_sinx_x(0.5*tB*qb);
131    const double siCt = tC*sas_sinx_x(0.5*tC*qc);
[fc0b7aa]132
133#if OVERLAPPING
[4493288]134    const double f = dr0*siA*siB*siC
135        + drA*(siAt-siA)*siB*siC
136        + drB*siAt*(siBt-siB)*siC
137        + drC*siAt*siBt*(siCt-siC);
[fc0b7aa]138#else
[4493288]139    const double f = dr0*siA*siB*siC
140        + drA*(siAt-siA)*siB*siC
141        + drB*siA*(siBt-siB)*siC
142        + drC*siA*siB*(siCt-siC);
[fc0b7aa]143#endif
[2a0b2b1]144
[44bd2be]145    return 1.0e-4 * f * f;
146}
Note: See TracBrowser for help on using the repository browser.