source: sasmodels/sasmodels/models/fcc.c @ d507c3a

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

Converted SCCrystalModel

  • Property mode set to 100644
File size: 5.4 KB
Line 
1double form_volume(double radius);
2double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld);
3double Iqxy(double qx, double qy, double dnn,
4    double d_factor, double radius,double sld, double solvent_sld,
5    double theta, double phi, double psi);
6
7double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi);
8double _FCCeval(double Theta, double Phi, double temp1, double temp3);
9
10
11double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) {
12
13        const double Da = d_factor*dnn;
14        const double temp1 = q*q*Da*Da;
15        const double temp3 = q*dnn;
16
17        double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI);
18        return(retVal);
19}
20
21double _FCCeval(double Theta, double Phi, double temp1, double temp3) {
22
23        double result;
24        double sin_theta,cos_theta,sin_phi,cos_phi;
25        SINCOS(Theta, sin_theta, cos_theta);
26        SINCOS(Phi, sin_phi, cos_phi);
27
28        const double temp6 =  sin_theta;
29        const double temp7 =  sin_theta*sin_phi + cos_theta;
30        const double temp8 = -sin_theta*cos_phi + cos_theta;
31        const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi;
32
33        const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9)));
34        result = pow((1.0-(temp10*temp10)),3)*temp6
35            / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10)
36              * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10)
37              * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10));
38
39        return (result);
40}
41
42double form_volume(double radius){
43    return 1.333333333333333*M_PI*radius*radius*radius;
44}
45
46
47double Iq(double q, double dnn,
48  double d_factor, double radius,
49  double sld, double solvent_sld){
50
51        //Volume fraction calculated from lattice symmetry and sphere radius
52        const double s1 = dnn*sqrt(2.0);
53        const double latticescale = 4.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
54
55    const double va = 0.0;
56    const double vb = 2.0*M_PI;
57    const double vaj = 0.0;
58    const double vbj = M_PI;
59
60    double summ = 0.0;
61    double answer = 0.0;
62        for(int i=0; i<150; i++) {
63                //setup inner integral over the ellipsoidal cross-section
64                double summj=0.0;
65                const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi
66                for(int j=0;j<150;j++) {
67                        //20 gauss points for the inner integral
68                        double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta
69                        double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi);
70                        summj += yyy;
71                }
72                //now calculate the value of the inner integral
73                double answer = (vbj-vaj)/2.0*summj;
74
75                //now calculate outer integral
76                summ = summ+(Gauss150Wt[i] * answer);
77        }               //final scaling is done at the end of the function, after the NT_FP64 case
78
79        answer = (vb-va)/2.0*summ;
80        answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale;
81
82    return answer;
83
84
85}
86
87
88double Iqxy(double qx, double qy, double dnn,
89    double d_factor, double radius,double sld, double solvent_sld,
90    double theta, double phi, double psi){
91
92  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
93  // double q_z;
94  double cos_val_b3, cos_val_b2, cos_val_b1;
95  double a1_dot_q, a2_dot_q,a3_dot_q;
96  double answer;
97  double Zq, Fkq, Fkq_2;
98
99  //convert to q and make scaled values
100  double q = sqrt(qx*qx+qy*qy);
101  double q_x = qx/q;
102  double q_y = qy/q;
103
104  //convert angle degree to radian
105  theta = theta * M_PI_180;
106  phi = phi * M_PI_180;
107  psi = psi * M_PI_180;
108
109  const double Da = d_factor*dnn;
110  const double s1 = dnn/sqrt(0.75);
111
112
113  //the occupied volume of the lattice
114  const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1);
115  // q vector
116  // q_z = 0.0; // for SANS; assuming qz is negligible
117  /// Angles here are respect to detector coordinate
118  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
119    // b3 axis orientation
120    b3_x = cos(theta) * cos(phi);
121    b3_y = sin(theta);
122    //b3_z = -cos(theta) * sin(phi);
123    cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
124
125    //alpha = acos(cos_val_b3);
126    // b1 axis orientation
127    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
128    b1_y = sin(psi)*cos(theta);
129    cos_val_b1 = b1_x*q_x + b1_y*q_y;
130    // b2 axis orientation
131    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
132        b2_y = cos(theta)*cos(psi);
133    cos_val_b2 = b2_x*q_x + b2_y*q_y;
134
135    // The following test should always pass
136    if (fabs(cos_val_b3)>1.0) {
137      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
138      cos_val_b3 = 1.0;
139    }
140    if (fabs(cos_val_b2)>1.0) {
141      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
142      cos_val_b2 = 1.0;
143    }
144    if (fabs(cos_val_b1)>1.0) {
145      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
146      cos_val_b1 = 1.0;
147    }
148    // Compute the angle btw vector q and the a3 axis
149    a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
150
151    // a1 axis
152    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
153
154    // a2 axis
155    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
156
157
158    // Get Fkq and Fkq_2
159    Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
160    Fkq_2 = Fkq*Fkq;
161    // Call Zq=Z1*Z2*Z3
162    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
163    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
164    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
165
166  // Use SphereForm directly from libigor
167  answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
168
169  return answer;
170 }
Note: See TracBrowser for help on using the repository browser.