source: sasmodels/sasmodels/models/core_shell_bicelle_elliptical.c @ 1e7b0db0

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 1e7b0db0 was 1e7b0db0, checked in by wojciech, 8 years ago

sinc tranfered to sas_sinx_x

  • Property mode set to 100644
File size: 5.4 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 
35                Iq(double qq,
36                   double rad,
37                   double x_core,
38                   double radthick,
39                   double facthick,
40                   double length,
41                   double rhoc,
42                   double rhoh,
43                   double rhor,
44                   double rhosolv)
45{
46    double si1,si2,be1,be2;
47     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle
48     // tested against limiting cases of cylinder, elliptical_cylinder and core_shell_bicelle
49     //    const double uplim = M_PI_4;
50    const double halfheight = 0.5*length;
51    //const double va = 0.0;
52    //const double vb = 1.0;
53    // inner integral limits
54    //const double vaj=0.0;
55    //const double vbj=M_PI;
56
57    const double radius_major = rad * x_core;
58    const double rA = 0.5*(square(radius_major) + square(rad));
59    const double rB = 0.5*(square(radius_major) - square(rad));
60    const double dr1 = (rhoc-rhoh)   *M_PI*rad*radius_major*(2.0*halfheight);;
61    const double dr2 = (rhor-rhosolv)*M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
62    const double dr3 = (rhoh-rhor)   *M_PI*rad*radius_major*2.0*(halfheight+facthick);
63    //const double vol1 = M_PI*rad*radius_major*(2.0*halfheight);
64    //const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
65    //const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick);
66
67    //initialize integral
68    double outer_sum = 0.0;
69    for(int i=0;i<76;i++) {
70        //setup inner integral over the ellipsoidal cross-section
71        // since we generate these lots of times, why not store them somewhere?
72        //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
73        const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0;
74        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
75        double inner_sum=0;
76        double sinarg1 = qq*halfheight*cos_alpha;
77        double sinarg2 = qq*(halfheight+facthick)*cos_alpha;
78        si1 = sas_sinx_x(sinarg1);
79        si2 = sas_sinx_x(sinarg2);
80        for(int j=0;j<76;j++) {
81            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
82            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
83            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2;
84            const double rr = sqrt(rA - rB*cos(beta));
85            double besarg1 = qq*rr*sin_alpha;
86            double besarg2 = qq*(rr+radthick)*sin_alpha;
87            be1 = sas_J1c(besarg1);
88            be2 = sas_J1c(besarg2);
89            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 +
90                                              dr2*si2*be2 +
91                                              dr3*si2*be1);
92        }
93        //now calculate outer integral
94        outer_sum += Gauss76Wt[i] * inner_sum;
95    }
96
97    return outer_sum*2.5e-05;
98}
99
100double 
101Iqxy(double qx, double qy,
102          double rad,
103          double x_core,
104          double radthick,
105          double facthick,
106          double length,
107          double rhoc,
108          double rhoh,
109          double rhor,
110          double rhosolv,
111          double theta,
112          double phi,
113          double psi)
114{
115       // THIS NEEDS TESTING
116    double qq, cos_val, cos_mu, cos_nu;
117    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, qq, cos_val, cos_mu, cos_nu);
118    const double dr1 = rhoc-rhoh;
119    const double dr2 = rhor-rhosolv;
120    const double dr3 = rhoh-rhor;
121    const double radius_major = rad*x_core;
122    const double halfheight = 0.5*length;
123    const double vol1 = M_PI*rad*radius_major*length;
124    const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
125    const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick);
126
127    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2)
128    // Given:    radius_major = r_ratio * radius_minor 
129    // ASSUME the sin_alpha is included in the separate integration over orientation of rod angle
130    const double r = rad*sqrt(square(x_core*cos_nu) + cos_mu*cos_mu);
131    const double be1 = sas_J1c(qq*r);
132    const double be2 = sas_J1c( qq*(r + radthick ) );
133    const double si1 = sas_sinx_x( qq*halfheight*cos_val );
134    const double si2 = sas_sinx_x( qq*(halfheight + facthick)*cos_val );
135    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1);
136    //const double vol = form_volume(radius_minor, r_ratio, length);
137    return 1.0e-4 * Aq;
138}
139
Note: See TracBrowser for help on using the repository browser.