source: sasmodels/sasmodels/models/core_shell_bicelle.c @ 2573fa1

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 2573fa1 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
RevLine 
[2222134]1double form_volume(double radius, double thick_rim, double thick_face, double length);
[8007311]2double Iq(double q,
3          double radius,
[2222134]4          double thick_rim,
5          double thick_face,
[8007311]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,
[2222134]15          double thick_rim,
16          double thick_face,
[8007311]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
[2222134]26double form_volume(double radius, double thick_rim, double thick_face, double length)
[8007311]27{
[5bddd89]28    return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face);
[8007311]29}
30
31static double
[b260926]32bicelle_kernel(double q,
[8007311]33              double rad,
34              double radthick,
35              double facthick,
[b260926]36              double halflength,
[8007311]37              double rhoc,
38              double rhoh,
39              double rhor,
40              double rhosolv,
[5bddd89]41              double sin_alpha,
42              double cos_alpha)
[8007311]43{
[e7678b2]44    const double dr1 = rhoc-rhoh;
45    const double dr2 = rhor-rhosolv;
46    const double dr3 = rhoh-rhor;
[b260926]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);
[e7678b2]55
56    const double t = vol1*dr1*si1*be1 +
57                     vol2*dr2*si2*be2 +
58                     vol3*dr3*si2*be1;
59
[b260926]60    const double retval = t*t;
[e7678b2]61
[5bddd89]62    return retval;
[8007311]63
64}
65
66static double
[b260926]67bicelle_integration(double q,
[8007311]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{
[e7678b2]77    // set up the integration end points
[5bddd89]78    const double uplim = M_PI_4;
[b260926]79    const double halflength = 0.5*length;
[e7678b2]80
81    double summ = 0.0;
82    for(int i=0;i<N_POINTS_76;i++) {
[5bddd89]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);
[b260926]86        double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick,
87                             halflength, rhoc, rhoh, rhor, rhosolv,
[5bddd89]88                             sin_alpha, cos_alpha);
[b260926]89        summ += yyy*sin_alpha;
[e7678b2]90    }
91
92    // calculate value of integral to return
93    double answer = uplim*summ;
[5bddd89]94    return answer;
[8007311]95}
96
97static double
[5bddd89]98bicelle_kernel_2d(double qx, double qy,
[8007311]99          double radius,
[2222134]100          double thick_rim,
101          double thick_face,
[8007311]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{
[5bddd89]110    double q, sin_alpha, cos_alpha;
111    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
[8007311]112
[2222134]113    double answer = bicelle_kernel(q, radius, thick_rim, thick_face,
[5bddd89]114                           0.5*length, core_sld, face_sld, rim_sld,
[b260926]115                           solvent_sld, sin_alpha, cos_alpha);
116    return 1.0e-4*answer;
[8007311]117}
118
119double Iq(double q,
120          double radius,
[2222134]121          double thick_rim,
122          double thick_face,
[8007311]123          double length,
124          double core_sld,
125          double face_sld,
126          double rim_sld,
127          double solvent_sld)
128{
[2222134]129    double intensity = bicelle_integration(q, radius, thick_rim, thick_face,
[8007311]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,
[2222134]137          double thick_rim,
138          double thick_face,
[8007311]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{
[5bddd89]147    double intensity = bicelle_kernel_2d(qx, qy,
[8007311]148                      radius,
[2222134]149                      thick_rim,
150                      thick_face,
[8007311]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.