[3271e20] | 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 | |
---|
[e7b3d7b] | 7 | double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); |
---|
| 8 | double _FCCeval(double Theta, double Phi, double temp1, double temp3); |
---|
[3271e20] | 9 | |
---|
| 10 | |
---|
[e7b3d7b] | 11 | double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { |
---|
[3271e20] | 12 | |
---|
| 13 | const double Da = d_factor*dnn; |
---|
| 14 | const double temp1 = q*q*Da*Da; |
---|
| 15 | const double temp3 = q*dnn; |
---|
| 16 | |
---|
[e7b3d7b] | 17 | double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); |
---|
[3271e20] | 18 | return(retVal); |
---|
| 19 | } |
---|
| 20 | |
---|
[e7b3d7b] | 21 | double _FCCeval(double Theta, double Phi, double temp1, double temp3) { |
---|
[3271e20] | 22 | |
---|
| 23 | double result; |
---|
[82d239a] | 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))); |
---|
[4962519] | 34 | result = cube(1.0-(temp10*temp10))*temp6 |
---|
[82d239a] | 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)); |
---|
[3271e20] | 38 | |
---|
| 39 | return (result); |
---|
| 40 | } |
---|
| 41 | |
---|
| 42 | double form_volume(double radius){ |
---|
[ad90df9] | 43 | return sphere_volume(radius); |
---|
[3271e20] | 44 | } |
---|
| 45 | |
---|
| 46 | |
---|
| 47 | double 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 |
---|
[e7b3d7b] | 52 | const double s1 = dnn*sqrt(2.0); |
---|
[ad90df9] | 53 | const double latticescale = 4.0*sphere_volume(radius/s1); |
---|
[3271e20] | 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 |
---|
[e7b3d7b] | 69 | double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi); |
---|
[3271e20] | 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; |
---|
[9aac25d] | 80 | answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; |
---|
[3271e20] | 81 | |
---|
| 82 | return answer; |
---|
| 83 | |
---|
| 84 | |
---|
| 85 | } |
---|
| 86 | |
---|
[11ca2ab] | 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 | { |
---|
[50beefe] | 92 | double q, zhat, yhat, xhat; |
---|
| 93 | ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); |
---|
[11ca2ab] | 94 | |
---|
[50beefe] | 95 | const double a1 = yhat + xhat; |
---|
| 96 | const double a2 = xhat + zhat; |
---|
| 97 | const double a3 = yhat + zhat; |
---|
[11ca2ab] | 98 | const double qd = 0.5*q*dnn; |
---|
[0b717c5] | 99 | const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); |
---|
| 100 | const double tanh_qd = tanh(arg); |
---|
| 101 | const double cosh_qd = cosh(arg); |
---|
| 102 | const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) |
---|
| 103 | * tanh_qd/(1. - cos(qd*a2)/cosh_qd) |
---|
| 104 | * tanh_qd/(1. - cos(qd*a3)/cosh_qd); |
---|
| 105 | |
---|
| 106 | //if (isnan(Zq)) printf("q:(%g,%g) qd: %g a1: %g a2: %g a3: %g arg: %g\n", qx, qy, qd, a1, a2, a3, arg); |
---|
[11ca2ab] | 107 | |
---|
| 108 | const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; |
---|
| 109 | //the occupied volume of the lattice |
---|
[ec9d329] | 110 | const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn); |
---|
[11ca2ab] | 111 | return lattice_scale * Fq; |
---|
| 112 | } |
---|