source: sasmodels/sasmodels/models/core_shell_bicelle.c @ 8007311

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

Converted CoreShellBicelle?
  • Property mode set to 100644
File size: 5.0 KB
Line 
1double form_volume(double radius, double rim_thickness, double face_thickness, double length);
2double Iq(double q,
3          double radius,
4          double rim_thickness,
5          double face_thickness,
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,
15          double rim_thickness,
16          double face_thickness,
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
26double form_volume(double radius, double rim_thickness, double face_thickness, double length)
27{
28    return M_PI*(radius+rim_thickness)*(radius+rim_thickness)*(length+2*face_thickness);
29}
30
31static double
32bicelle_kernel(double qq,
33              double rad,
34              double radthick,
35              double facthick,
36              double length,
37              double rhoc,
38              double rhoh,
39              double rhor,
40              double rhosolv,
41              double dum)
42{
43        double dr1,dr2,dr3;
44        double besarg1,besarg2;
45        double vol1,vol2,vol3;
46        double sinarg1,sinarg2;
47        double t1,t2,t3;
48        double retval,si1,si2,be1,be2;
49
50        dr1 = rhoc-rhoh;
51        dr2 = rhor-rhosolv;
52        dr3=  rhoh-rhor;
53        vol1 = M_PI*rad*rad*(2.0*length);
54        vol2 = M_PI*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
55        vol3= M_PI*(rad)*(rad)*(2.0*length+2.0*facthick);
56        besarg1 = qq*rad*sin(dum);
57        besarg2 = qq*(rad+radthick)*sin(dum);
58        sinarg1 = qq*length*cos(dum);
59        sinarg2 = qq*(length+facthick)*cos(dum);
60
61        if(besarg1 == 0) {
62                be1 = 0.5;
63        } else {
64                be1 = J1(besarg1)/besarg1;
65        }
66        if(besarg2 == 0) {
67                be2 = 0.5;
68        } else {
69                be2 = J1(besarg2)/besarg2;
70        }
71        if(sinarg1 == 0) {
72                si1 = 1.0;
73        } else {
74                si1 = sin(sinarg1)/sinarg1;
75        }
76        if(sinarg2 == 0) {
77                si2 = 1.0;
78        } else {
79                si2 = sin(sinarg2)/sinarg2;
80        }
81        t1 = 2.0*vol1*dr1*si1*be1;
82        t2 = 2.0*vol2*dr2*si2*be2;
83        t3 = 2.0*vol3*dr3*si2*be1;
84
85        retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum);
86        return(retval);
87
88}
89
90static double
91bicelle_integration(double qq,
92                   double rad,
93                   double radthick,
94                   double facthick,
95                   double length,
96                   double rhoc,
97                   double rhoh,
98                   double rhor,
99                   double rhosolv)
100{
101
102
103        double answer,halfheight;
104        double lolim,uplim,summ,yyy,zi;
105        int nord,i;
106
107        // set up the integration end points
108        nord = 76;
109        lolim = 0.0;
110        uplim = M_PI/2;
111        halfheight = length/2.0;
112
113        summ = 0.0;
114        i=0;
115        for(i=0;i<nord;i++) {
116                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
117                yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick,
118                                     halfheight, rhoc, rhoh, rhor,rhosolv, zi);
119                summ += yyy;
120        }
121
122        // calculate value of integral to return
123        answer = (uplim-lolim)/2.0*summ;
124        return(answer);
125}
126
127static double
128bicelle_kernel_2d(double q, double q_x, double q_y,
129          double radius,
130          double rim_thickness,
131          double face_thickness,
132          double length,
133          double core_sld,
134          double face_sld,
135          double rim_sld,
136          double solvent_sld,
137          double theta,
138          double phi)
139{
140    double cyl_x, cyl_y;
141    double alpha, cos_val;
142    double answer;
143
144    //convert angle degree to radian
145    theta *= M_PI/180.0;
146    phi *= M_PI/180.0;
147
148    // Cylinder orientation
149    cyl_x = cos(theta) * cos(phi);
150    cyl_y = sin(theta);
151
152    // Compute the angle btw vector q and the axis of the cylinder
153    cos_val = cyl_x*q_x + cyl_y*q_y;
154    alpha = acos( cos_val );
155
156    // Get the kernel
157    answer = bicelle_kernel(q, radius, rim_thickness, face_thickness,
158                           length/2.0, core_sld, face_sld, rim_sld,
159                           solvent_sld, alpha) / fabs(sin(alpha));
160
161    answer *= 1.0e-4;
162
163    return answer;
164}
165
166double Iq(double q,
167          double radius,
168          double rim_thickness,
169          double face_thickness,
170          double length,
171          double core_sld,
172          double face_sld,
173          double rim_sld,
174          double solvent_sld)
175{
176    double intensity = bicelle_integration(q, radius, rim_thickness, face_thickness,
177                       length, core_sld, face_sld, rim_sld, solvent_sld);
178    return intensity*1.0e-4;
179}
180
181
182double Iqxy(double qx, double qy,
183          double radius,
184          double rim_thickness,
185          double face_thickness,
186          double length,
187          double core_sld,
188          double face_sld,
189          double rim_sld,
190          double solvent_sld,
191          double theta,
192          double phi)
193{
194    double q;
195    q = sqrt(qx*qx+qy*qy);
196    double intensity = bicelle_kernel_2d(q, qx/q, qy/q,
197                      radius,
198                      rim_thickness,
199                      face_thickness,
200                      length,
201                      core_sld,
202                      face_sld,
203                      rim_sld,
204                      solvent_sld,
205                      theta,
206                      phi);
207
208    return intensity;
209}
Note: See TracBrowser for help on using the repository browser.