source: sasmodels/sasmodels/models/core_shell_bicelle.c @ 592343f

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 592343f was 592343f, checked in by wojciech, 7 years ago

sas_J1c translated to sas_2J1x_x
  • Property mode set to 100644
File size: 4.6 KB
Line 
1double form_volume(double radius, double thick_rim, double thick_face, double length);
2double Iq(double q,
3          double radius,
4          double thick_rim,
5          double thick_face,
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 thick_rim,
16          double thick_face,
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 thick_rim, double thick_face, double length)
27{
28    return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face);
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 sin_alpha,
42              double cos_alpha)
43{
44    double si1,si2,be1,be2;
45
46    const double dr1 = rhoc-rhoh;
47    const double dr2 = rhor-rhosolv;
48    const double dr3 = rhoh-rhor;
49    const double vol1 = M_PI*rad*rad*(2.0*length);
50    const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick);
51    const double vol3 = M_PI*rad*rad*2.0*(length+facthick);
52    double besarg1 = qq*rad*sin_alpha;
53    double besarg2 = qq*(rad+radthick)*sin_alpha;
54    double sinarg1 = qq*length*cos_alpha;
55    double sinarg2 = qq*(length+facthick)*cos_alpha;
56
57    be1 = sas_2J1x_x(besarg1);
58    be2 = sas_2J1x_x(besarg2);
59    si1 = sas_sinx_x(sinarg1);
60    si2 = sas_sinx_x(sinarg2);
61
62    const double t = vol1*dr1*si1*be1 +
63                     vol2*dr2*si2*be2 +
64                     vol3*dr3*si2*be1;
65
66    const double retval = t*t*sin_alpha;
67
68    return retval;
69
70}
71
72static double
73bicelle_integration(double qq,
74                   double rad,
75                   double radthick,
76                   double facthick,
77                   double length,
78                   double rhoc,
79                   double rhoh,
80                   double rhor,
81                   double rhosolv)
82{
83    // set up the integration end points
84    const double uplim = M_PI_4;
85    const double halfheight = 0.5*length;
86
87    double summ = 0.0;
88    for(int i=0;i<N_POINTS_76;i++) {
89        double alpha = (Gauss76Z[i] + 1.0)*uplim;
90        double sin_alpha, cos_alpha; // slots to hold sincos function output
91        SINCOS(alpha, sin_alpha, cos_alpha);
92        double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick,
93                             halfheight, rhoc, rhoh, rhor, rhosolv,
94                             sin_alpha, cos_alpha);
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 qx, double qy,
105          double radius,
106          double thick_rim,
107          double thick_face,
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    double q, sin_alpha, cos_alpha;
117    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
118
119    double answer = bicelle_kernel(q, radius, thick_rim, thick_face,
120                           0.5*length, core_sld, face_sld, rim_sld,
121                           solvent_sld, sin_alpha, cos_alpha) / fabs(sin_alpha);
122
123    answer *= 1.0e-4;
124
125    return answer;
126}
127
128double Iq(double q,
129          double radius,
130          double thick_rim,
131          double thick_face,
132          double length,
133          double core_sld,
134          double face_sld,
135          double rim_sld,
136          double solvent_sld)
137{
138    double intensity = bicelle_integration(q, radius, thick_rim, thick_face,
139                       length, core_sld, face_sld, rim_sld, solvent_sld);
140    return intensity*1.0e-4;
141}
142
143
144double Iqxy(double qx, double qy,
145          double radius,
146          double thick_rim,
147          double thick_face,
148          double length,
149          double core_sld,
150          double face_sld,
151          double rim_sld,
152          double solvent_sld,
153          double theta,
154          double phi)
155{
156    double intensity = bicelle_kernel_2d(qx, qy,
157                      radius,
158                      thick_rim,
159                      thick_face,
160                      length,
161                      core_sld,
162                      face_sld,
163                      rim_sld,
164                      solvent_sld,
165                      theta,
166                      phi);
167
168    return intensity;
169}
Note: See TracBrowser for help on using the repository browser.