source: sasmodels/sasmodels/models/core_shell_ellipsoid.c @ 58c3367

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 58c3367 was 29172aa, checked in by richardh, 8 years ago

rename params and sort equations in docs for core_shell_ellipsoid
  • Property mode set to 100644
File size: 4.4 KB
RevLine 
[81dd619]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,
[29172aa]10          double sld_core,
11          double sld_shell,
12          double sld_solvent);
[81dd619]13
14
15double Iqxy(double qx, double qy,
16          double equat_core,
17          double polar_core,
18          double equat_shell,
19          double polar_shell,
[29172aa]20          double sld_core,
21          double sld_shell,
22          double sld_solvent,
[81dd619]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,
[29172aa]42          double sld_core,
43          double sld_shell,
44          double sld_solvent)
[81dd619]45{
46
[e7678b2]47    //upper and lower integration limits
48    const double lolim = 0.0;
49    const double uplim = 1.0;
[81dd619]50
[e7678b2]51    double summ = 0.0;   //initialize intergral
[81dd619]52
[29172aa]53    const double delpc = sld_core - sld_shell;    //core - shell
54    const double delps = sld_shell - sld_solvent; //shell - solvent
[81dd619]55
[e7678b2]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    }
[81dd619]68
[e7678b2]69    double answer = (uplim-lolim)/2.0*summ;
[81dd619]70
[e7678b2]71    //convert to [cm-1]
72    answer *= 1.0e-4;
[81dd619]73
[e7678b2]74    return answer;
[81dd619]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,
[29172aa]83          double sld_core,
84          double sld_shell,
85          double sld_solvent,
[81dd619]86          double theta,
87          double phi)
88{
89    //convert angle degree to radian
[e7678b2]90    theta = theta * M_PI_180;
91    phi = phi * M_PI_180;
[81dd619]92
93
94    // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis.
[e7678b2]95    const double cyl_x = cos(theta) * cos(phi);
96    const double cyl_y = sin(theta);
[81dd619]97
[29172aa]98    const double sldcs = sld_core - sld_shell;
99    const double sldss = sld_shell- sld_solvent;
[81dd619]100
101    // Compute the angle btw vector q and the
102    // axis of the cylinder
[e7678b2]103    const double cos_val = cyl_x*q_x + cyl_y*q_y;
[81dd619]104
105    // Call the IGOR library function to get the kernel: MUST use gfn4 not gf2 because of the def of params.
[e7678b2]106    double answer = gfn4(cos_val,
[81dd619]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,
[29172aa]126          double sld_core,
127          double sld_shell,
128          double sld_solvent)
[81dd619]129{
130    double intensity = core_shell_ellipsoid_kernel(q,
131           equat_core,
132           polar_core,
133           equat_shell,
134           polar_shell,
[29172aa]135           sld_core,
136           sld_shell,
137           sld_solvent);
[81dd619]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,
[29172aa]148          double sld_core,
149          double sld_shell,
150          double sld_solvent,
[81dd619]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,
[29172aa]161                       sld_core,
162                       sld_shell,
163                       sld_solvent,
[81dd619]164                       theta,
165                       phi);
166
167    return intensity;
168}
Note: See TracBrowser for help on using the repository browser.