source: sasmodels/sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c @ 108e70e

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

use Iqac/Iqabc? for the new orientation interface but Iqxy for the old

  • Property mode set to 100644
File size: 4.6 KB
Line 
1// NOTE that "length" here is the full height of the core!
2static double
3form_volume(double r_minor,
4        double x_core,
5        double thick_rim,
6        double thick_face,
7        double length)
8{
9    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length +
10                 square(r_minor)*x_core*2.0*thick_face  );
11}
12
13static double
14Iq(double q,
15        double r_minor,
16        double x_core,
17        double thick_rim,
18        double thick_face,
19        double length,
20        double rhoc,
21        double rhoh,
22        double rhor,
23        double rhosolv,
24        double sigma)
25{
26    double si1,si2,be1,be2;
27     // core_shell_bicelle_elliptical_belt, RKH 5th Oct 2017, core_shell_bicelle_elliptical
28     // tested briefly against limiting cases of cylinder, hollow cylinder & elliptical cylinder models
29     //    const double uplim = M_PI_4;
30    const double halfheight = 0.5*length;
31    //const double va = 0.0;
32    //const double vb = 1.0;
33    // inner integral limits
34    //const double vaj=0.0;
35    //const double vbj=M_PI;
36
37    const double r_major = r_minor * x_core;
38    const double r2A = 0.5*(square(r_major) + square(r_minor));
39    const double r2B = 0.5*(square(r_major) - square(r_minor));
40    // dr1,2,3 are now for Vcore, Vcore+rim, Vcore+face,
41    const double dr1 = (-rhor - rhoh + rhoc + rhosolv) *M_PI*r_minor*r_major*
42          2.0*halfheight;
43    const double dr2 = (rhor-rhosolv) *M_PI*(r_minor+thick_rim)*(
44         r_major+thick_rim)* 2.0*halfheight;
45    const double dr3 = (rhoh-rhosolv) *M_PI*r_minor*r_major*
46         2.0*(halfheight+thick_face);
47    //initialize integral
48    double outer_sum = 0.0;
49    for(int i=0;i<GAUSS_N;i++) {
50        //setup inner integral over the ellipsoidal cross-section
51        // since we generate these lots of times, why not store them somewhere?
52        //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
53        const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0;
54        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
55        double inner_sum=0;
56        double sinarg1 = q*halfheight*cos_alpha;
57        double sinarg2 = q*(halfheight+thick_face)*cos_alpha;
58        si1 = sas_sinx_x(sinarg1);
59        si2 = sas_sinx_x(sinarg2);
60        for(int j=0;j<GAUSS_N;j++) {
61            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
62            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
63            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2;
64            const double rr = sqrt(r2A - r2B*cos(beta));
65            double besarg1 = q*rr*sin_alpha;
66            double besarg2 = q*(rr+thick_rim)*sin_alpha;
67            be1 = sas_2J1x_x(besarg1);
68            be2 = sas_2J1x_x(besarg2);
69            inner_sum += GAUSS_W[j] *square(dr1*si1*be1 +
70                                              dr2*si1*be2 +
71                                              dr3*si2*be1);
72        }
73        //now calculate outer integral
74        outer_sum += GAUSS_W[i] * inner_sum;
75    }
76
77    return outer_sum*2.5e-05*exp(-0.5*square(q*sigma));
78}
79
80static double
81Iqabc(double qa, double qb, double qc,
82          double r_minor,
83          double x_core,
84          double thick_rim,
85          double thick_face,
86          double length,
87          double rhoc,
88          double rhoh,
89          double rhor,
90          double rhosolv,
91          double sigma)
92{
93    // integrated 2d seems to match 1d reasonably well, except perhaps at very high Q
94    // Vol1,2,3 and dr1,2,3 are now for Vcore, Vcore+rim, Vcore+face,
95    const double dr1 = -rhor - rhoh + rhoc + rhosolv;
96    const double dr2 = rhor-rhosolv;
97    const double dr3 = rhoh-rhosolv;
98    const double r_major = r_minor*x_core;
99    const double halfheight = 0.5*length;
100    const double vol1 = M_PI*r_minor*r_major*length;
101    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*halfheight;
102    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
103
104    // Compute effective radius in rotated coordinates
105    const double qr_hat = sqrt(square(r_major*qb) + square(r_minor*qa));
106    // does this need to be changed for the "missing corners" where there there is no "belt" ?
107    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qb)
108                                   + square((r_minor+thick_rim)*qa));
109    const double be1 = sas_2J1x_x( qr_hat );
110    const double be2 = sas_2J1x_x( qrshell_hat );
111    const double si1 = sas_sinx_x( halfheight*qc );
112    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc );
113    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1);
114    return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma));
115}
Note: See TracBrowser for help on using the repository browser.