source: sasmodels/sasmodels/models/core_shell_ellipsoid.c @ 2222134

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 2222134 was 2222134, checked in by butler, 8 years ago

Updating parameter names regarding #649
  • Property mode set to 100644
File size: 4.7 KB
RevLine 
[2222134]1double form_volume(double radius_equat_core,
2                   double radius_polar_core,
3                   double radius_equat_shell,
4                   double radius_polar_shell);
[81dd619]5double Iq(double q,
[2222134]6          double radius_equat_core,
7          double radius_polar_core,
8          double radius_equat_shell,
9          double radius_polar_shell,
[29172aa]10          double sld_core,
11          double sld_shell,
12          double sld_solvent);
[81dd619]13
14
15double Iqxy(double qx, double qy,
[2222134]16          double radius_equat_core,
17          double radius_polar_core,
18          double radius_equat_shell,
19          double radius_polar_shell,
[29172aa]20          double sld_core,
21          double sld_shell,
22          double sld_solvent,
[81dd619]23          double theta,
24          double phi);
25
26
[2222134]27double form_volume(double radius_equat_core,
28                   double radius_polar_core,
29                   double radius_equat_shell,
30                   double radius_polar_shell)
[81dd619]31{
[2222134]32    double vol = 4.0*M_PI/3.0*radius_equat_shell*radius_equat_shell*radius_polar_shell;
[81dd619]33    return vol;
34}
35
36static double
37core_shell_ellipsoid_kernel(double q,
[2222134]38          double radius_equat_core,
39          double radius_polar_core,
40          double radius_equat_shell,
41          double radius_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,
[2222134]59                                  radius_equat_core,
60                                  radius_polar_core,
61                                  radius_equat_shell,
62                                  radius_polar_shell,
[e7678b2]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,
[2222134]79          double radius_equat_core,
80          double radius_polar_core,
81          double radius_equat_shell,
82          double radius_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,
[2222134]107                  radius_equat_core,
108                  radius_polar_core,
109                  radius_equat_shell,
110                  radius_polar_shell,
[81dd619]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,
[2222134]122          double radius_equat_core,
123          double radius_polar_core,
124          double radius_equat_shell,
125          double radius_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,
[2222134]131           radius_equat_core,
132           radius_polar_core,
133           radius_equat_shell,
134           radius_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,
[2222134]144          double radius_equat_core,
145          double radius_polar_core,
146          double radius_equat_shell,
147          double radius_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,
[2222134]157                       radius_equat_core,
158                       radius_polar_core,
159                       radius_equat_shell,
160                       radius_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.