1 | double form_volume(double radius); |
---|
2 | double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); |
---|
3 | double 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 | |
---|
7 | double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); |
---|
8 | double _BCCeval(double Theta, double Phi, double temp1, double temp3); |
---|
9 | double _sphereform(double q, double radius, double sld, double solvent_sld); |
---|
10 | |
---|
11 | |
---|
12 | double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { |
---|
13 | |
---|
14 | const double Da = d_factor*dnn; |
---|
15 | const double temp1 = q*q*Da*Da; |
---|
16 | const double temp3 = q*dnn; |
---|
17 | |
---|
18 | double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); |
---|
19 | return(retVal); |
---|
20 | } |
---|
21 | |
---|
22 | double _BCCeval(double Theta, double Phi, double temp1, double temp3) { |
---|
23 | |
---|
24 | double result; |
---|
25 | double sin_theta,cos_theta,sin_phi,cos_phi; |
---|
26 | SINCOS(Theta, sin_theta, cos_theta); |
---|
27 | SINCOS(Phi, sin_phi, cos_phi); |
---|
28 | |
---|
29 | const double temp6 = sin_theta; |
---|
30 | const double temp7 = sin_theta*cos_phi + sin_theta*sin_phi + cos_theta; |
---|
31 | const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; |
---|
32 | const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; |
---|
33 | |
---|
34 | const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); |
---|
35 | result = cube(1.0 - (temp10*temp10))*temp6 |
---|
36 | / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) |
---|
37 | * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) |
---|
38 | * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); |
---|
39 | |
---|
40 | return (result); |
---|
41 | } |
---|
42 | |
---|
43 | double form_volume(double radius){ |
---|
44 | return sphere_volume(radius); |
---|
45 | } |
---|
46 | |
---|
47 | |
---|
48 | double Iq(double q, double dnn, |
---|
49 | double d_factor, double radius, |
---|
50 | double sld, double solvent_sld){ |
---|
51 | |
---|
52 | //Volume fraction calculated from lattice symmetry and sphere radius |
---|
53 | const double s1 = dnn/sqrt(0.75); |
---|
54 | const double latticescale = 2.0*sphere_volume(radius/s1); |
---|
55 | |
---|
56 | const double va = 0.0; |
---|
57 | const double vb = 2.0*M_PI; |
---|
58 | const double vaj = 0.0; |
---|
59 | const double vbj = M_PI; |
---|
60 | |
---|
61 | double summ = 0.0; |
---|
62 | double answer = 0.0; |
---|
63 | for(int i=0; i<150; i++) { |
---|
64 | //setup inner integral over the ellipsoidal cross-section |
---|
65 | double summj=0.0; |
---|
66 | const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; //the outer dummy is phi |
---|
67 | for(int j=0;j<150;j++) { |
---|
68 | //20 gauss points for the inner integral |
---|
69 | double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the inner dummy is theta |
---|
70 | double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi); |
---|
71 | summj += yyy; |
---|
72 | } |
---|
73 | //now calculate the value of the inner integral |
---|
74 | double answer = (vbj-vaj)/2.0*summj; |
---|
75 | |
---|
76 | //now calculate outer integral |
---|
77 | summ = summ+(Gauss150Wt[i] * answer); |
---|
78 | } //final scaling is done at the end of the function, after the NT_FP64 case |
---|
79 | |
---|
80 | answer = (vb-va)/2.0*summ; |
---|
81 | answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; |
---|
82 | |
---|
83 | return answer; |
---|
84 | } |
---|
85 | |
---|
86 | |
---|
87 | double Iqxy(double qx, double qy, |
---|
88 | double dnn, double d_factor, double radius, |
---|
89 | double sld, double solvent_sld, |
---|
90 | double theta, double phi, double psi) |
---|
91 | { |
---|
92 | double q, zhat, yhat, xhat; |
---|
93 | ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); |
---|
94 | |
---|
95 | const double a1 = +xhat - zhat + yhat; |
---|
96 | const double a2 = +xhat + zhat - yhat; |
---|
97 | const double a3 = -xhat + zhat + yhat; |
---|
98 | |
---|
99 | const double qd = 0.5*q*dnn; |
---|
100 | const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); |
---|
101 | const double tanh_qd = tanh(arg); |
---|
102 | const double cosh_qd = cosh(arg); |
---|
103 | const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) |
---|
104 | * tanh_qd/(1. - cos(qd*a2)/cosh_qd) |
---|
105 | * tanh_qd/(1. - cos(qd*a3)/cosh_qd); |
---|
106 | |
---|
107 | const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; |
---|
108 | //the occupied volume of the lattice |
---|
109 | const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn); |
---|
110 | return lattice_scale * Fq; |
---|
111 | } |
---|