source: sasmodels/sasmodels/models/elliptical_cylinder.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 40a87fa, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

lint and latex cleanup

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