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