Changes in sasmodels/models/core_shell_parallelepiped.c [c69d6d6:e077231] in sasmodels
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/core_shell_parallelepiped.c
rc69d6d6 re077231 1 double form_volume(double length_a, double length_b, double length_c, 2 double thick_rim_a, double thick_rim_b, double thick_rim_c); 3 double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 4 double solvent_sld, double length_a, double length_b, double length_c, 5 double thick_rim_a, double thick_rim_b, double thick_rim_c); 6 double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld, 7 double crim_sld, double solvent_sld, double length_a, double length_b, 8 double length_c, double thick_rim_a, double thick_rim_b, 9 double thick_rim_c, double theta, double phi, double psi); 10 11 double form_volume(double length_a, double length_b, double length_c, 12 double thick_rim_a, double thick_rim_b, double thick_rim_c) 1 // Set OVERLAPPING to 1 in order to fill in the edges of the box, with 2 // c endcaps and b overlapping a. With the proper choice of parameters, 3 // (setting rim slds to sld, core sld to solvent, rim thickness to thickness 4 // and subtracting 2*thickness from length, this should match the hollow 5 // rectangular prism.) Set it to 0 for the documented behaviour. 6 #define OVERLAPPING 0 7 static double 8 form_volume(double length_a, double length_b, double length_c, 9 double thick_rim_a, double thick_rim_b, double thick_rim_c) 13 10 { 14 //return length_a * length_b * length_c; 15 return length_a * length_b * length_c + 16 2.0 * thick_rim_a * length_b * length_c + 17 2.0 * thick_rim_b * length_a * length_c + 18 2.0 * thick_rim_c * length_a * length_b; 11 return 12 #if OVERLAPPING 13 // Hollow rectangular prism only includes the volume of the shell 14 // so uncomment the next line when comparing. Solid rectangular 15 // prism, or parallelepiped want filled cores, so comment when 16 // comparing. 17 //-length_a * length_b * length_c + 18 (length_a + 2.0*thick_rim_a) * 19 (length_b + 2.0*thick_rim_b) * 20 (length_c + 2.0*thick_rim_c); 21 #else 22 length_a * length_b * length_c + 23 2.0 * thick_rim_a * length_b * length_c + 24 2.0 * length_a * thick_rim_b * length_c + 25 2.0 * length_a * length_b * thick_rim_c; 26 #endif 19 27 } 20 28 21 double Iq(double q, 29 static double 30 Iq(double q, 22 31 double core_sld, 23 32 double arim_sld, … … 32 41 double thick_rim_c) 33 42 { 34 // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c _scaled43 // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 35 44 // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 36 //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 45 // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker) 46 // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK) 37 47 38 const double mu = 0.5 * q * length_b;48 const double half_q = 0.5*q; 39 49 40 //calculate volume before rescaling (in original code, but not used) 41 //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 42 //double vol = length_a * length_b * length_c + 43 // 2.0 * thick_rim_a * length_b * length_c + 44 // 2.0 * thick_rim_b * length_a * length_c + 45 // 2.0 * thick_rim_c * length_a * length_b; 50 const double tA = length_a + 2.0*thick_rim_a; 51 const double tB = length_b + 2.0*thick_rim_b; 52 const double tC = length_c + 2.0*thick_rim_c; 46 53 47 // Scale sides by B 48 const double a_scaled = length_a / length_b; 49 const double c_scaled = length_c / length_b; 50 51 double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b; 52 double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b; 53 double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present 54 55 double Vin = length_a * length_b * length_c; 56 //double Vot = (length_a * length_b * length_c + 57 // 2.0 * thick_rim_a * length_b * length_c + 58 // 2.0 * length_a * thick_rim_b * length_c + 59 // 2.0 * length_a * length_b * thick_rim_c); 60 double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 61 double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 62 double V3 = (2.0 * length_a * length_b * thick_rim_c); //not present 63 double Vot = Vin + V1 + V2 + V3; 64 65 // Scale factors (note that drC is not used later) 66 const double drho0 = (core_sld-solvent_sld); 67 const double drhoA = (arim_sld-solvent_sld); 68 const double drhoB = (brim_sld-solvent_sld); 69 const double drhoC = (crim_sld-solvent_sld); // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; 70 71 72 // Precompute scale factors for combining cross terms from the shape 73 const double scale23 = drhoA*V1; 74 const double scale14 = drhoB*V2; 75 const double scale24 = drhoC*V3; 76 const double scale11 = drho0*Vin; 77 const double scale12 = drho0*Vin - scale23 - scale14 - scale24; 54 // Scale factors 55 const double dr0 = (core_sld-solvent_sld); 56 const double drA = (arim_sld-solvent_sld); 57 const double drB = (brim_sld-solvent_sld); 58 const double drC = (crim_sld-solvent_sld); 78 59 79 60 // outer integral (with gauss points), integration limits = 0, 1 80 double outer_total = 0; //initialize integral 61 double outer_sum = 0; //initialize integral 62 for( int i=0; i<GAUSS_N; i++) { 63 const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 64 const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 81 65 82 for( int i=0; i<76; i++) { 83 double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 84 double mu_proj = mu * sqrt(1.0-sigma*sigma); 66 // inner integral (with gauss points), integration limits = 0, pi/2 67 const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 68 const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 69 double inner_sum = 0.0; 70 for(int j=0; j<GAUSS_N; j++) { 71 const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 72 double sin_beta, cos_beta; 73 SINCOS(M_PI_2*beta, sin_beta, cos_beta); 74 const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 75 const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 76 const double siAt = tA * sas_sinx_x(tA * mu * sin_beta); 77 const double siBt = tB * sas_sinx_x(tB * mu * cos_beta); 85 78 86 // inner integral (with gauss points), integration limits = 0, 1 87 double inner_total = 0.0;88 double inner_total_crim = 0.0;89 for(int j=0; j<76; j++) {90 const double uu = 0.5 * ( Gauss76Z[j] + 1.0);91 double sin_uu, cos_uu; 92 SINCOS(M_PI_2*uu, sin_uu, cos_uu);93 const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);94 const double si2 = sas_sinx_x(mu_proj * cos_uu );95 const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);96 const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 79 #if OVERLAPPING 80 const double f = dr0*siA*siB*siC 81 + drA*(siAt-siA)*siB*siC 82 + drB*siAt*(siBt-siB)*siC 83 + drC*siAt*siBt*(siCt-siC); 84 #else 85 const double f = dr0*siA*siB*siC 86 + drA*(siAt-siA)*siB*siC 87 + drB*siA*(siBt-siB)*siC 88 + drC*siA*siB*(siCt-siC); 89 #endif 97 90 98 // Expression in libCylinder.c (neither drC nor Vot are used) 99 const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 100 const double form_crim = scale11*si1*si2; 101 102 103 // correct FF : sum of square of phase factors 104 inner_total += Gauss76Wt[j] * form * form; 105 inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; 91 inner_sum += GAUSS_W[j] * f * f; 106 92 } 107 inner_total *= 0.5; 108 inner_total_crim *= 0.5; 93 inner_sum *= 0.5; 109 94 // now sum up the outer integral 110 const double si = sas_sinx_x(mu * c_scaled * sigma); 111 const double si_crim = sas_sinx_x(mu * tc * sigma); 112 outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim); 95 outer_sum += GAUSS_W[i] * inner_sum; 113 96 } 114 outer_ total*= 0.5;97 outer_sum *= 0.5; 115 98 116 99 //convert from [1e-12 A-1] to [cm-1] 117 return 1.0e-4 * outer_ total;100 return 1.0e-4 * outer_sum; 118 101 } 119 102 120 double Iqxy(double qx, double qy, 103 static double 104 Iqabc(double qa, double qb, double qc, 121 105 double core_sld, 122 106 double arim_sld, … … 129 113 double thick_rim_a, 130 114 double thick_rim_b, 131 double thick_rim_c, 132 double theta, 133 double phi, 134 double psi) 115 double thick_rim_c) 135 116 { 136 double q, zhat, yhat, xhat;137 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);138 139 117 // cspkernel in csparallelepiped recoded here 140 118 const double dr0 = core_sld-solvent_sld; … … 143 121 const double drC = crim_sld-solvent_sld; 144 122 145 double Vin = length_a * length_b * length_c;146 double V1 = 2.0 * thick_rim_a * length_b * length_c; // incorrect V1(aa*bb*cc+2*ta*bb*cc)147 double V2 = 2.0 * length_a * thick_rim_b * length_c; // incorrect V2(aa*bb*cc+2*aa*tb*cc)148 double V3 = 2.0 * length_a * length_b * thick_rim_c;149 // As for the 1D case, Vot is not used150 //double Vot = (length_a * length_b * length_c +151 // 2.0 * thick_rim_a * length_b * length_c +152 // 2.0 * length_a * thick_rim_b * length_c +153 // 2.0 * length_a * length_b * thick_rim_c);123 const double tA = length_a + 2.0*thick_rim_a; 124 const double tB = length_b + 2.0*thick_rim_b; 125 const double tC = length_c + 2.0*thick_rim_c; 126 const double siA = length_a*sas_sinx_x(0.5*length_a*qa); 127 const double siB = length_b*sas_sinx_x(0.5*length_b*qb); 128 const double siC = length_c*sas_sinx_x(0.5*length_c*qc); 129 const double siAt = tA*sas_sinx_x(0.5*tA*qa); 130 const double siBt = tB*sas_sinx_x(0.5*tB*qb); 131 const double siCt = tC*sas_sinx_x(0.5*tC*qc); 154 132 155 // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 156 // the scaling by B. 157 double ta = length_a + 2.0*thick_rim_a; 158 double tb = length_b + 2.0*thick_rim_b; 159 double tc = length_c + 2.0*thick_rim_c; 160 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 161 double siA = sas_sinx_x(0.5*q*length_a*xhat); 162 double siB = sas_sinx_x(0.5*q*length_b*yhat); 163 double siC = sas_sinx_x(0.5*q*length_c*zhat); 164 double siAt = sas_sinx_x(0.5*q*ta*xhat); 165 double siBt = sas_sinx_x(0.5*q*tb*yhat); 166 double siCt = sas_sinx_x(0.5*q*tc*zhat); 167 168 169 // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 170 // in the 1D code, but should be checked! 171 double f = ( dr0*siA*siB*siC*Vin 172 + drA*(siAt-siA)*siB*siC*V1 173 + drB*siA*(siBt-siB)*siC*V2 174 + drC*siA*siB*(siCt-siC)*V3); 133 #if OVERLAPPING 134 const double f = dr0*siA*siB*siC 135 + drA*(siAt-siA)*siB*siC 136 + drB*siAt*(siBt-siB)*siC 137 + drC*siAt*siBt*(siCt-siC); 138 #else 139 const double f = dr0*siA*siB*siC 140 + drA*(siAt-siA)*siB*siC 141 + drB*siA*(siBt-siB)*siC 142 + drC*siA*siB*(siCt-siC); 143 #endif 175 144 176 145 return 1.0e-4 * f * f;
Note: See TracChangeset
for help on using the changeset viewer.