source: sasmodels/sasmodels/models/core_shell_bicelle.c @ 43b7eea

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 43b7eea was 43b7eea, checked in by wojciech, 8 years ago

Initial clean-up of Bessel functions

  • Property mode set to 100644
File size: 4.8 KB
Line 
1double form_volume(double radius, double rim_thickness, double face_thickness, double length);
2double Iq(double q,
3          double radius,
4          double rim_thickness,
5          double face_thickness,
6          double length,
7          double core_sld,
8          double face_sld,
9          double rim_sld,
10          double solvent_sld);
11
12
13double Iqxy(double qx, double qy,
14          double radius,
15          double rim_thickness,
16          double face_thickness,
17          double length,
18          double core_sld,
19          double face_sld,
20          double rim_sld,
21          double solvent_sld,
22          double theta,
23          double phi);
24
25
26double form_volume(double radius, double rim_thickness, double face_thickness, double length)
27{
28    return M_PI*(radius+rim_thickness)*(radius+rim_thickness)*(length+2*face_thickness);
29}
30
31static double
32bicelle_kernel(double qq,
33              double rad,
34              double radthick,
35              double facthick,
36              double length,
37              double rhoc,
38              double rhoh,
39              double rhor,
40              double rhosolv,
41              double dum)
42{
43    double si1,si2,be1,be2;
44
45    const double dr1 = rhoc-rhoh;
46    const double dr2 = rhor-rhosolv;
47    const double dr3 = rhoh-rhor;
48    const double vol1 = M_PI*rad*rad*(2.0*length);
49    const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick);
50    const double vol3 = M_PI*rad*rad*2.0*(length+facthick);
51    double sn,cn;
52    SINCOS(dum, sn, cn);
53    double besarg1 = qq*rad*sn;
54    double besarg2 = qq*(rad+radthick)*sn;
55    double sinarg1 = qq*length*cn;
56    double sinarg2 = qq*(length+facthick)*cn;
57
58    be1 = sas_J1c(besarg1);
59    be2 = sas_J1c(besarg2);
60    si1 = sinc(sinarg1);
61    si2 = sinc(sinarg2);
62
63    const double t = vol1*dr1*si1*be1 +
64                     vol2*dr2*si2*be2 +
65                     vol3*dr3*si2*be1;
66
67    const double retval = t*t*sn;
68
69    return(retval);
70
71}
72
73static double
74bicelle_integration(double qq,
75                   double rad,
76                   double radthick,
77                   double facthick,
78                   double length,
79                   double rhoc,
80                   double rhoh,
81                   double rhor,
82                   double rhosolv)
83{
84    // set up the integration end points
85    const double uplim = M_PI/4;
86    const double halfheight = length/2.0;
87
88    double summ = 0.0;
89    for(int i=0;i<N_POINTS_76;i++) {
90        double zi = (Gauss76Z[i] + 1.0)*uplim;
91        double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick,
92                             halfheight, rhoc, rhoh, rhor,rhosolv, zi);
93        summ += yyy;
94    }
95
96    // calculate value of integral to return
97    double answer = uplim*summ;
98    return(answer);
99}
100
101static double
102bicelle_kernel_2d(double q, double q_x, double q_y,
103          double radius,
104          double rim_thickness,
105          double face_thickness,
106          double length,
107          double core_sld,
108          double face_sld,
109          double rim_sld,
110          double solvent_sld,
111          double theta,
112          double phi)
113{
114    //convert angle degree to radian
115    theta *= M_PI_180;
116    phi *= M_PI_180;
117
118    // Cylinder orientation
119    const double cyl_x = cos(theta) * cos(phi);
120    const double cyl_y = sin(theta);
121
122    // Compute the angle btw vector q and the axis of the cylinder
123    const double cos_val = cyl_x*q_x + cyl_y*q_y;
124    const double alpha = acos( cos_val );
125
126    // Get the kernel
127    double answer = bicelle_kernel(q, radius, rim_thickness, face_thickness,
128                           length/2.0, core_sld, face_sld, rim_sld,
129                           solvent_sld, alpha) / fabs(sin(alpha));
130
131    answer *= 1.0e-4;
132
133    return answer;
134}
135
136double Iq(double q,
137          double radius,
138          double rim_thickness,
139          double face_thickness,
140          double length,
141          double core_sld,
142          double face_sld,
143          double rim_sld,
144          double solvent_sld)
145{
146    double intensity = bicelle_integration(q, radius, rim_thickness, face_thickness,
147                       length, core_sld, face_sld, rim_sld, solvent_sld);
148    return intensity*1.0e-4;
149}
150
151
152double Iqxy(double qx, double qy,
153          double radius,
154          double rim_thickness,
155          double face_thickness,
156          double length,
157          double core_sld,
158          double face_sld,
159          double rim_sld,
160          double solvent_sld,
161          double theta,
162          double phi)
163{
164    double q;
165    q = sqrt(qx*qx+qy*qy);
166    double intensity = bicelle_kernel_2d(q, qx/q, qy/q,
167                      radius,
168                      rim_thickness,
169                      face_thickness,
170                      length,
171                      core_sld,
172                      face_sld,
173                      rim_sld,
174                      solvent_sld,
175                      theta,
176                      phi);
177
178    return intensity;
179}
Note: See TracBrowser for help on using the repository browser.