source: sasmodels/sasmodels/models/elliptical_cylinder.c @ a8b3cdb

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since a8b3cdb was a8b3cdb, checked in by Doucet, Mathieu <doucetm@…>, 8 years ago

Converted elliptical cylinder model

  • Property mode set to 100644
File size: 5.4 KB
Line 
1double form_volume(double r_minor, double r_ratio, double length);
2double Iq(double q, double r_minor, double r_ratio, double length,
3          double sld, double solvent_sld);
4double Iqxy(double qx, double qy, double r_minor, double r_ratio, double length,
5            double sld, double solvent_sld, double theta, double phi, double psi);
6
7
8double _elliptical_cylinder_kernel(double q, double r_minor, double r_ratio, double theta);
9
10double _elliptical_cylinder_kernel(double q, double r_minor, double r_ratio, double theta)
11{
12    // This is the function LAMBDA1^2 in Feigin's notation
13    // q is the q-value for the calculation (1/A)
14    // r_minor is the transformed radius"a" in Feigin's notation
15    // r_ratio is the ratio (major radius)/(minor radius) of the Ellipsoid [=] ---
16    // theta is the dummy variable of the integration
17
18    double retval,arg;
19
20    arg = q*r_minor*sqrt((1.0+r_ratio*r_ratio)/2+(1.0-r_ratio*r_ratio)*cos(theta)/2);
21    if (arg == 0.0){
22        retval = 1.0;
23    }else{
24        retval = 2.0*NR_BessJ1(arg)/arg;
25    }
26    return retval*retval ;
27}
28
29
30double form_volume(double r_minor, double r_ratio, double length)
31{
32    return M_PI * r_minor * r_minor * r_ratio * length;
33}
34
35double Iq(double q, double r_minor, double r_ratio, double length,
36          double sld, double solvent_sld) {
37
38    const int nordi=76; //order of integration
39    const int nordj=20;
40    double va,vb;       //upper and lower integration limits
41    double summ,zi,yyy,answer;         //running tally of integration
42    double summj,vaj,vbj,zij,arg,si;            //for the inner integration
43
44    // orientational average limits
45    va = 0.0;
46    vb = 1.0;
47    // inner integral limits
48    vaj=0.0;
49    vbj=M_PI;
50
51    //initialize integral
52    summ = 0.0;
53
54    const double delrho = sld - solvent_sld;
55
56    for(int i=0;i<nordi;i++) {
57        //setup inner integral over the ellipsoidal cross-section
58        summj=0;
59        zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;     //the "x" dummy
60        arg = r_minor*sqrt(1.0-zi*zi);
61        for(int j=0;j<nordj;j++) {
62            //20 gauss points for the inner integral
63            zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0;        //the "y" dummy
64            yyy = Gauss20Wt[j] * _elliptical_cylinder_kernel(q, arg, r_ratio, zij);
65            summj += yyy;
66        }
67        //now calculate the value of the inner integral
68        answer = (vbj-vaj)/2.0*summj;
69        //divide integral by Pi
70        answer /= M_PI;
71
72        //now calculate outer integral
73        arg = q*length*zi/2.0;
74        if (arg == 0.0){
75            si = 1.0;
76        }else{
77            si = sin(arg) * sin(arg) / arg / arg;
78        }
79        yyy = Gauss76Wt[i] * answer * si;
80        summ += yyy;
81    }
82
83    answer = (vb-va)/2.0*summ;
84    // Multiply by contrast^2
85    answer *= delrho*delrho;
86
87    const double vol = form_volume(r_minor, r_ratio, length);
88    return answer*vol*vol*1e-4;
89}
90
91
92double Iqxy(double qx, double qy, double r_minor, double r_ratio, double length,
93            double sld, double solvent_sld, double theta, double phi, double psi) {
94    const double _theta = theta * M_PI / 180.0;
95    const double _phi = phi * M_PI / 180.0;
96    const double _psi = psi * M_PI / 180.0;
97    const double q = sqrt(qx*qx+qy*qy);
98    const double q_x = qx/q;
99    const double q_y = qy/q;
100
101    //Cylinder orientation
102    double cyl_x = cos(_theta) * cos(_phi);
103    double cyl_y = sin(_theta);
104
105    //cyl_z = -cos(_theta) * sin(_phi);
106
107    // q vector
108    //q_z = 0;
109
110    // Note: cos(alpha) = 0 and 1 will get an
111    // undefined value from CylKernel
112    //alpha = acos( cos_val );
113
114    //ellipse orientation:
115    // the elliptical corss section was transformed and projected
116    // into the detector plane already through sin(alpha)and furthermore psi remains as same
117    // on the detector plane.
118    // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt
119    // the wave vector q.
120
121    //x- y- component on the detector plane.
122    const double ella_x =  -cos(_phi)*sin(_psi) * sin(_theta)+sin(_phi)*cos(_psi);
123    const double ella_y =  sin(_psi)*cos(_theta);
124    const double ellb_x =  -sin(_theta)*cos(_psi)*cos(_phi)-sin(_psi)*sin(_phi);
125    const double ellb_y =  cos(_theta)*cos(_psi);
126
127    // Compute the angle btw vector q and the
128    // axis of the cylinder
129    double cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
130
131    // calculate the axis of the ellipse wrt q-coord.
132    double cos_nu = ella_x*q_x + ella_y*q_y;
133    double cos_mu = ellb_x*q_x + ellb_y*q_y;
134
135    // The following test should always pass
136    if (fabs(cos_val)>1.0) {
137      //printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
138      cos_val = 1.0;
139    }
140    if (fabs(cos_nu)>1.0) {
141      //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
142      cos_nu = 1.0;
143    }
144    if (fabs(cos_mu)>1.0) {
145      //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
146      cos_mu = 1.0;
147    }
148
149    const double r_major = r_ratio * r_minor;
150    const double qr = q*sqrt( r_major*r_major*cos_nu*cos_nu + r_minor*r_minor*cos_mu*cos_mu );
151    const double qL = q*length*cos_val/2.0;
152
153    double Be;
154    if (qr==0){
155      Be = 0.5;
156    }else{
157      Be = NR_BessJ1(qr)/qr;
158    }
159
160    double Si;
161    if (qL==0){
162      Si = 1.0;
163    }else{
164      Si = sin(qL)/qL;
165    }
166
167    const double k = 2.0 * Be * Si;
168    const double vol = form_volume(r_minor, r_ratio, length);
169    return (sld - solvent_sld) * (sld - solvent_sld) * k * k *vol*vol*1.0e-4;
170}
Note: See TracBrowser for help on using the repository browser.