Changeset 2a0b2b1 in sasmodels for sasmodels/models/sc_paracrystal.c
- Timestamp:
- Apr 14, 2017 6:30:29 AM (7 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- fdd56a1
- Parents:
- 9901384
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/sc_paracrystal.c
r50beefe r2a0b2b1 1 double form_volume(double radius); 1 static double 2 _sq_sc(double qa, double qb, double qc, double dnn, double d_factor) 3 { 4 // Rewriting equations for efficiency, accuracy and readability, and so 5 // code is reusable between 1D and 2D models. 6 const double a1 = qa; 7 const double a2 = qb; 8 const double a3 = qc; 2 9 3 double Iq(double q, 4 double dnn, 5 double d_factor, 6 double radius, 7 double sphere_sld, 8 double solvent_sld); 10 const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 9 11 10 double Iqxy(double qx, double qy, 11 double dnn, 12 double d_factor, 13 double radius, 14 double sphere_sld, 15 double solvent_sld, 16 double theta, 17 double phi, 18 double psi); 12 // Numerator: (1 - exp(a)^2)^3 13 // => (-(exp(2a) - 1))^3 14 // => -expm1(2a)^3 15 // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2) 16 // => exp(a)^2 - 2 cos(xk) exp(a) + 1 17 // => (exp(a) - 2 cos(xk)) * exp(a) + 1 18 const double exp_arg = exp(-arg); 19 const double Sq = -cube(expm1(-2.0*arg)) 20 / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 21 * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 22 * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 19 23 20 double form_volume(double radius) 24 return Sq; 25 } 26 27 // occupied volume fraction calculated from lattice symmetry and sphere radius 28 static double 29 _sc_volume_fraction(double radius, double dnn) 30 { 31 return sphere_volume(radius/dnn); 32 } 33 34 static double 35 form_volume(double radius) 21 36 { 22 37 return sphere_volume(radius); 23 38 } 24 39 25 static double 26 sc_eval(double theta, double phi, double temp3, double temp4, double temp5) 40 41 static double Iq(double q, double dnn, 42 double d_factor, double radius, 43 double sld, double solvent_sld) 27 44 { 28 double cnt, snt; 29 SINCOS(theta, cnt, snt); 45 // translate a point in [-1,1] to a point in [0, 2 pi] 46 const double phi_m = M_PI_4; 47 const double phi_b = M_PI_4; 48 // translate a point in [-1,1] to a point in [0, pi] 49 const double theta_m = M_PI_4; 50 const double theta_b = M_PI_4; 30 51 31 double cnp, snp;32 SINCOS(phi, cnp, snp);33 52 34 double temp6 = snt; 35 double temp7 = -1.0*temp3*snt*cnp; 36 double temp8 = temp3*snt*snp; 37 double temp9 = temp3*cnt; 38 double result = temp6/((1.0-temp4*cos((temp7))+temp5)* 39 (1.0-temp4*cos((temp8))+temp5)* 40 (1.0-temp4*cos((temp9))+temp5)); 41 return (result); 53 double outer_sum = 0.0; 54 for(int i=0; i<150; i++) { 55 double inner_sum = 0.0; 56 const double theta = Gauss150Z[i]*theta_m + theta_b; 57 double sin_theta, cos_theta; 58 SINCOS(theta, sin_theta, cos_theta); 59 const double qc = q*cos_theta; 60 const double qab = q*sin_theta; 61 for(int j=0;j<150;j++) { 62 const double phi = Gauss150Z[j]*phi_m + phi_b; 63 double sin_phi, cos_phi; 64 SINCOS(phi, sin_phi, cos_phi); 65 const double qa = qab*cos_phi; 66 const double qb = qab*sin_phi; 67 const double fq = _sq_sc(qa, qb, qc, dnn, d_factor); 68 inner_sum += Gauss150Wt[j] * fq; 69 } 70 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 71 outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 72 } 73 outer_sum *= theta_m; 74 const double Sq = outer_sum/M_PI_2; 75 const double Pq = sphere_form(q, radius, sld, solvent_sld); 76 77 return _sc_volume_fraction(radius, dnn) * Pq * Sq; 42 78 } 43 79 44 static double45 sc_integrand(double dnn, double d_factor, double qq, double xx, double yy)46 {47 //Function to calculate integrand values for simple cubic structure48 80 49 double da = d_factor*dnn; 50 double temp1 = qq*qq*da*da; 51 double temp2 = cube(-expm1(-temp1)); 52 double temp3 = qq*dnn; 53 double temp4 = 2.0*exp(-0.5*temp1); 54 double temp5 = exp(-1.0*temp1); 55 56 double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2; 57 58 return(integrand); 59 } 60 61 double Iq(double q, 62 double dnn, 63 double d_factor, 64 double radius, 65 double sphere_sld, 66 double solvent_sld) 67 { 68 const double va = 0.0; 69 const double vb = M_PI_2; //orientation average, outer integral 70 71 double summ=0.0; 72 double answer=0.0; 73 74 for(int i=0;i<150;i++) { 75 //setup inner integral over the ellipsoidal cross-section 76 double summj=0.0; 77 double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; 78 for(int j=0;j<150;j++) { 79 //150 gauss points for the inner integral 80 double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0; 81 double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij); 82 summj += tmp; 83 } 84 //now calculate the value of the inner integral 85 answer = (vb-va)/2.0*summj; 86 87 //now calculate outer integral 88 double tmp = Gauss150Wt[i] * answer; 89 summ += tmp; 90 } //final scaling is done at the end of the function, after the NT_FP64 case 91 92 answer = (vb-va)/2.0*summ; 93 94 //Volume fraction calculated from lattice symmetry and sphere radius 95 // NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3 96 const double latticeScale = sphere_volume(radius/dnn); 97 98 answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale; 99 100 return answer; 101 } 102 103 double Iqxy(double qx, double qy, 104 double dnn, 105 double d_factor, 106 double radius, 107 double sphere_sld, 108 double solvent_sld, 109 double theta, 110 double phi, 111 double psi) 81 static double Iqxy(double qx, double qy, 82 double dnn, double d_factor, double radius, 83 double sld, double solvent_sld, 84 double theta, double phi, double psi) 112 85 { 113 86 double q, zhat, yhat, xhat; 114 87 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 88 const double qa = q*xhat; 89 const double qb = q*yhat; 90 const double qc = q*zhat; 115 91 116 const double qd = q*dnn; 117 const double arg = 0.5*square(qd*d_factor); 118 const double tanh_qd = tanh(arg); 119 const double cosh_qd = cosh(arg); 120 const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd) 121 * tanh_qd/(1. - cos(qd*yhat)/cosh_qd) 122 * tanh_qd/(1. - cos(qd*xhat)/cosh_qd); 123 124 const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 125 //the occupied volume of the lattice 126 const double lattice_scale = sphere_volume(radius/dnn); 127 return lattice_scale * Fq; 92 q = sqrt(qa*qa + qb*qb + qc*qc); 93 const double Pq = sphere_form(q, radius, sld, solvent_sld); 94 const double Sq = _sq_sc(qa, qb, qc, dnn, d_factor); 95 return _sc_volume_fraction(radius, dnn) * Pq * Sq; 128 96 }
Note: See TracChangeset
for help on using the changeset viewer.