source: sasmodels/sasmodels/models/core_shell_ellipsoid.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.4 KB
Line 
1double form_volume(double equat_core,
2                   double polar_core,
3                   double equat_shell,
4                   double polar_shell);
5double Iq(double q,
6          double equat_core,
7          double polar_core,
8          double equat_shell,
9          double polar_shell,
10          double core_sld,
11          double shell_sld,
12          double solvent_sld);
13
14
15double Iqxy(double qx, double qy,
16          double equat_core,
17          double polar_core,
18          double equat_shell,
19          double polar_shell,
20          double core_sld,
21          double shell_sld,
22          double solvent_sld,
23          double theta,
24          double phi);
25
26
27double form_volume(double equat_core,
28                   double polar_core,
29                   double equat_shell,
30                   double polar_shell)
31{
32    double vol = 4.0*M_PI/3.0*equat_shell*equat_shell*polar_shell;
33    return vol;
34}
35
36static double
37core_shell_ellipsoid_kernel(double q,
38          double equat_core,
39          double polar_core,
40          double equat_shell,
41          double polar_shell,
42          double core_sld,
43          double shell_sld,
44          double solvent_sld)
45{
46
47    //upper and lower integration limits
48    const double lolim = 0.0;
49    const double uplim = 1.0;
50
51    double summ = 0.0;   //initialize intergral
52
53    const double delpc = core_sld - shell_sld;    //core - shell
54    const double delps = shell_sld - solvent_sld; //shell - solvent
55
56    for(int i=0;i<N_POINTS_76;i++) {
57        double zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
58        double yyy = Gauss76Wt[i] * gfn4(zi,
59                                  equat_core,
60                                  polar_core,
61                                  equat_shell,
62                                  polar_shell,
63                                  delpc,
64                                  delps,
65                                  q);
66        summ += yyy;
67    }
68
69    double answer = (uplim-lolim)/2.0*summ;
70
71    //convert to [cm-1]
72    answer *= 1.0e-4;
73
74    return answer;
75}
76
77static double
78core_shell_ellipsoid_kernel_2d(double q, double q_x, double q_y,
79          double equat_core,
80          double polar_core,
81          double equat_shell,
82          double polar_shell,
83          double core_sld,
84          double shell_sld,
85          double solvent_sld,
86          double theta,
87          double phi)
88{
89    //convert angle degree to radian
90    theta = theta * M_PI_180;
91    phi = phi * M_PI_180;
92
93
94    // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis.
95    const double cyl_x = cos(theta) * cos(phi);
96    const double cyl_y = sin(theta);
97
98    const double sldcs = core_sld - shell_sld;
99    const double sldss = shell_sld- solvent_sld;
100
101    // Compute the angle btw vector q and the
102    // axis of the cylinder
103    const double cos_val = cyl_x*q_x + cyl_y*q_y;
104
105    // Call the IGOR library function to get the kernel: MUST use gfn4 not gf2 because of the def of params.
106    double answer = gfn4(cos_val,
107                  equat_core,
108                  polar_core,
109                  equat_shell,
110                  polar_shell,
111                  sldcs,
112                  sldss,
113                  q);
114
115    //convert to [cm-1]
116    answer *= 1.0e-4;
117
118    return answer;
119}
120
121double Iq(double q,
122          double equat_core,
123          double polar_core,
124          double equat_shell,
125          double polar_shell,
126          double core_sld,
127          double shell_sld,
128          double solvent_sld)
129{
130    double intensity = core_shell_ellipsoid_kernel(q,
131           equat_core,
132           polar_core,
133           equat_shell,
134           polar_shell,
135           core_sld,
136           shell_sld,
137           solvent_sld);
138
139    return intensity;
140}
141
142
143double Iqxy(double qx, double qy,
144          double equat_core,
145          double polar_core,
146          double equat_shell,
147          double polar_shell,
148          double core_sld,
149          double shell_sld,
150          double solvent_sld,
151          double theta,
152          double phi)
153{
154    double q;
155    q = sqrt(qx*qx+qy*qy);
156    double intensity = core_shell_ellipsoid_kernel_2d(q, qx/q, qy/q,
157                       equat_core,
158                       polar_core,
159                       equat_shell,
160                       polar_shell,
161                       core_sld,
162                       shell_sld,
163                       solvent_sld,
164                       theta,
165                       phi);
166
167    return intensity;
168}
Note: See TracBrowser for help on using the repository browser.