source: sasmodels/sasmodels/models/core_shell_bicelle_elliptical.c @ f4f85b3

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

core shell bicelle elliptical: adapt to new orientation parameters; adjust effective shell radius in 2D model so it matches core shell bicelle when axis ratio is 1

  • Property mode set to 100644
File size: 5.5 KB
Line 
1double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length);
2double Iq(double q,
3          double radius,
4          double x_core,
5          double thick_rim,
6          double thick_face,
7          double length,
8          double core_sld,
9          double face_sld,
10          double rim_sld,
11          double solvent_sld);
12
13
14double Iqxy(double qx, double qy,
15          double radius,
16          double x_core,
17          double thick_rim,
18          double thick_face,
19          double length,
20          double core_sld,
21          double face_sld,
22          double rim_sld,
23          double solvent_sld,
24          double theta,
25          double phi,
26          double psi);
27
28// NOTE that "length" here is the full height of the core!
29double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length)
30{
31    return M_PI*(radius+thick_rim)*(radius*x_core+thick_rim)*(length+2.0*thick_face);
32}
33
34double Iq(double q,
35          double rad,
36          double x_core,
37          double radthick,
38          double facthick,
39          double length,
40          double rhoc,
41          double rhoh,
42          double rhor,
43          double rhosolv)
44{
45    double si1,si2,be1,be2;
46     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle
47     // tested against limiting cases of cylinder, elliptical_cylinder and core_shell_bicelle
48     //    const double uplim = M_PI_4;
49    const double halfheight = 0.5*length;
50    //const double va = 0.0;
51    //const double vb = 1.0;
52    // inner integral limits
53    //const double vaj=0.0;
54    //const double vbj=M_PI;
55
56    const double radius_major = rad * x_core;
57    const double rA = 0.5*(square(radius_major) + square(rad));
58    const double rB = 0.5*(square(radius_major) - square(rad));
59    const double dr1 = (rhoc-rhoh)   *M_PI*rad*radius_major*(2.0*halfheight);;
60    const double dr2 = (rhor-rhosolv)*M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
61    const double dr3 = (rhoh-rhor)   *M_PI*rad*radius_major*2.0*(halfheight+facthick);
62    //const double vol1 = M_PI*rad*radius_major*(2.0*halfheight);
63    //const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
64    //const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick);
65
66    //initialize integral
67    double outer_sum = 0.0;
68    for(int i=0;i<76;i++) {
69        //setup inner integral over the ellipsoidal cross-section
70        // since we generate these lots of times, why not store them somewhere?
71        //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
72        const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0;
73        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
74        double inner_sum=0;
75        double sinarg1 = q*halfheight*cos_alpha;
76        double sinarg2 = q*(halfheight+facthick)*cos_alpha;
77        si1 = sas_sinx_x(sinarg1);
78        si2 = sas_sinx_x(sinarg2);
79        for(int j=0;j<76;j++) {
80            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
81            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
82            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2;
83            const double rr = sqrt(rA - rB*cos(beta));
84            double besarg1 = q*rr*sin_alpha;
85            double besarg2 = q*(rr+radthick)*sin_alpha;
86            be1 = sas_2J1x_x(besarg1);
87            be2 = sas_2J1x_x(besarg2);
88            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 +
89                                              dr2*si2*be2 +
90                                              dr3*si2*be1);
91        }
92        //now calculate outer integral
93        outer_sum += Gauss76Wt[i] * inner_sum;
94    }
95
96    return outer_sum*2.5e-05;
97}
98
99double 
100Iqxy(double qx, double qy,
101          double rad,
102          double x_core,
103          double radthick,
104          double facthick,
105          double length,
106          double rhoc,
107          double rhoh,
108          double rhor,
109          double rhosolv,
110          double theta,
111          double phi,
112          double psi)
113{
114       // THIS NEEDS TESTING
115    double q, xhat, yhat, zhat;
116    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
117    const double dr1 = rhoc-rhoh;
118    const double dr2 = rhor-rhosolv;
119    const double dr3 = rhoh-rhor;
120    const double radius_major = rad*x_core;
121    const double halfheight = 0.5*length;
122    const double vol1 = M_PI*rad*radius_major*length;
123    const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
124    const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick);
125
126    // Compute:  r = sqrt((radius_major*zhat)^2 + (radius_minor*yhat)^2)
127    // Given:    radius_major = r_ratio * radius_minor 
128    // ASSUME the sin_alpha is included in the separate integration over orientation of rod angle
129    const double rad_minor = rad;
130    const double rad_major = rad*x_core;
131    const double r_hat = sqrt(square(rad_major*xhat) + square(rad_minor*yhat));
132    const double rshell_hat = sqrt(square((rad_major+radthick)*xhat)
133                                   + square((rad_minor+radthick)*yhat));
134    const double be1 = sas_2J1x_x( q*r_hat );
135    const double be2 = sas_2J1x_x( q*rshell_hat );
136    const double si1 = sas_sinx_x( q*halfheight*zhat );
137    const double si2 = sas_sinx_x( q*(halfheight + facthick)*zhat );
138    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1);
139    //const double vol = form_volume(radius_minor, r_ratio, length);
140    return 1.0e-4 * Aq;
141}
142
Note: See TracBrowser for help on using the repository browser.