source: sasmodels/sasmodels/models/core_shell_bicelle.c @ e7678b2

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

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