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