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

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

core shell bicelle: code cleanup
  • Property mode set to 100644
File size: 4.5 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 q,
33              double rad,
34              double radthick,
35              double facthick,
36              double halflength,
37              double rhoc,
38              double rhoh,
39              double rhor,
40              double rhosolv,
41              double sin_alpha,
42              double cos_alpha)
43{
44    const double dr1 = rhoc-rhoh;
45    const double dr2 = rhor-rhosolv;
46    const double dr3 = rhoh-rhor;
47    const double vol1 = M_PI*square(rad)*2.0*(halflength);
48    const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick);
49    const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick);
50
51    const double be1 = sas_2J1x_x(q*(rad)*sin_alpha);
52    const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha);
53    const double si1 = sas_sinx_x(q*(halflength)*cos_alpha);
54    const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha);
55
56    const double t = vol1*dr1*si1*be1 +
57                     vol2*dr2*si2*be2 +
58                     vol3*dr3*si2*be1;
59
60    const double retval = t*t;
61
62    return retval;
63
64}
65
66static double
67bicelle_integration(double q,
68                   double rad,
69                   double radthick,
70                   double facthick,
71                   double length,
72                   double rhoc,
73                   double rhoh,
74                   double rhor,
75                   double rhosolv)
76{
77    // set up the integration end points
78    const double uplim = M_PI_4;
79    const double halflength = 0.5*length;
80
81    double summ = 0.0;
82    for(int i=0;i<N_POINTS_76;i++) {
83        double alpha = (Gauss76Z[i] + 1.0)*uplim;
84        double sin_alpha, cos_alpha; // slots to hold sincos function output
85        SINCOS(alpha, sin_alpha, cos_alpha);
86        double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick,
87                             halflength, rhoc, rhoh, rhor, rhosolv,
88                             sin_alpha, cos_alpha);
89        summ += yyy*sin_alpha;
90    }
91
92    // calculate value of integral to return
93    double answer = uplim*summ;
94    return answer;
95}
96
97static double
98bicelle_kernel_2d(double qx, double qy,
99          double radius,
100          double thick_rim,
101          double thick_face,
102          double length,
103          double core_sld,
104          double face_sld,
105          double rim_sld,
106          double solvent_sld,
107          double theta,
108          double phi)
109{
110    double q, sin_alpha, cos_alpha;
111    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
112
113    double answer = bicelle_kernel(q, radius, thick_rim, thick_face,
114                           0.5*length, core_sld, face_sld, rim_sld,
115                           solvent_sld, sin_alpha, cos_alpha);
116    return 1.0e-4*answer;
117}
118
119double Iq(double q,
120          double radius,
121          double thick_rim,
122          double thick_face,
123          double length,
124          double core_sld,
125          double face_sld,
126          double rim_sld,
127          double solvent_sld)
128{
129    double intensity = bicelle_integration(q, radius, thick_rim, thick_face,
130                       length, core_sld, face_sld, rim_sld, solvent_sld);
131    return intensity*1.0e-4;
132}
133
134
135double Iqxy(double qx, double qy,
136          double radius,
137          double thick_rim,
138          double thick_face,
139          double length,
140          double core_sld,
141          double face_sld,
142          double rim_sld,
143          double solvent_sld,
144          double theta,
145          double phi)
146{
147    double intensity = bicelle_kernel_2d(qx, qy,
148                      radius,
149                      thick_rim,
150                      thick_face,
151                      length,
152                      core_sld,
153                      face_sld,
154                      rim_sld,
155                      solvent_sld,
156                      theta,
157                      phi);
158
159    return intensity;
160}
Note: See TracBrowser for help on using the repository browser.