Changeset 14838a3 in sasmodels for sasmodels/models/core_shell_parallelepiped.c
- Timestamp:
- Oct 15, 2016 10:37:09 AM (7 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 2941abf
- Parents:
- 4962519
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/core_shell_parallelepiped.c
r3a48772 r14838a3 35 35 // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 36 36 37 double t1, t2, t3, t4, tmp, answer; 38 double mu = q * length_b; 37 const double mu = 0.5 * q * length_b; 39 38 40 39 //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); 40 //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 42 41 //double vol = length_a * length_b * length_c + 43 42 // 2.0 * thick_rim_a * length_b * length_c + … … 46 45 47 46 // Scale sides by B 48 double a_scaled = length_a / length_b;49 double c_scaled = length_c / length_b;47 const double a_scaled = length_a / length_b; 48 const double c_scaled = length_c / length_b; 50 49 51 // DelRho values (note that drC is not used later) 52 double dr0 = core_sld-solvent_sld; 53 double drA = arim_sld-solvent_sld; 54 double drB = brim_sld-solvent_sld; 55 //double drC = crim_sld-solvent_sld; 50 // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) 51 // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c, 52 // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 53 // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!) 54 double ta = (a_scaled + 2.0*thick_rim_a)/length_b; 55 double tb = (a_scaled + 2.0*thick_rim_b)/length_b; 56 56 57 //Order of integration 58 int nordi=76; 59 int nordj=76; 60 61 // outer integral (with nordi gauss points), integration limits = 0, 1 62 double summ = 0; //initialize integral 57 double Vin = length_a * length_b * length_c; 58 //double Vot = (length_a * length_b * length_c + 59 // 2.0 * thick_rim_a * length_b * length_c + 60 // 2.0 * length_a * thick_rim_b * length_c + 61 // 2.0 * length_a * length_b * thick_rim_c); 62 double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 63 double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 63 64 64 for( int i=0; i<nordi; i++) { 65 66 // inner integral (with nordj gauss points), integration limits = 0, 1 67 68 double summj = 0.0; 69 double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 70 71 for(int j=0; j<nordj; j++) { 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 drC_Vot = (crim_sld-solvent_sld)*Vot; 72 70 73 double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 74 double mudum = mu * sqrt(1.0-sigma*sigma); 71 // Precompute scale factors for combining cross terms from the shape 72 const double scale23 = drhoA*V1; 73 const double scale14 = drhoB*V2; 74 const double scale12 = drho0*Vin - scale23 - scale14; 75 75 76 double Vin = length_a * length_b * length_c; 77 //double Vot = (length_a * length_b * length_c + 78 // 2.0 * thick_rim_a * length_b * length_c + 79 // 2.0 * length_a * thick_rim_b * length_c + 80 // 2.0 * length_a * length_b * thick_rim_c); 81 double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 82 double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 83 84 // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) 85 // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c, 86 // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 87 // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!) 88 double ta = (a_scaled+2.0*thick_rim_a)/length_b; 89 double tb = (a_scaled+2.0*thick_rim_b)/length_b; 90 91 double arg1 = (0.5*mudum*a_scaled) * sin(M_PI_2*uu); 92 double arg2 = (0.5*mudum) * cos(M_PI_2*uu); 93 double arg3= (0.5*mudum*ta) * sin(M_PI_2*uu); 94 double arg4= (0.5*mudum*tb) * cos(M_PI_2*uu); 76 // outer integral (with gauss points), integration limits = 0, 1 77 double outer_total = 0; //initialize integral 95 78 96 if(arg1==0.0){ 97 t1 = 1.0; 98 } else { 99 t1 = (sin(arg1)/arg1); //defn for CSPP model sin(arg1)/arg1 test: (sin(arg1)/arg1)*(sin(arg1)/arg1) 100 } 101 if(arg2==0.0){ 102 t2 = 1.0; 103 } else { 104 t2 = (sin(arg2)/arg2); //defn for CSPP model sin(arg2)/arg2 test: (sin(arg2)/arg2)*(sin(arg2)/arg2) 105 } 106 if(arg3==0.0){ 107 t3 = 1.0; 108 } else { 109 t3 = sin(arg3)/arg3; 110 } 111 if(arg4==0.0){ 112 t4 = 1.0; 113 } else { 114 t4 = sin(arg4)/arg4; 115 } 116 79 for( int i=0; i<76; i++) { 80 double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 81 double mu_proj = mu * sqrt(1.0-sigma*sigma); 82 83 // inner integral (with gauss points), integration limits = 0, 1 84 double inner_total = 0.0; 85 for(int j=0; j<76; j++) { 86 const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 87 double sin_uu, cos_uu; 88 SINCOS(M_PI_2*uu, sin_uu, cos_uu); 89 const double si1 = sinc(mu_proj * sin_uu * a_scaled); 90 const double si2 = sinc(mu_proj * cos_uu); 91 const double si3 = sinc(mu_proj * sin_uu * ta); 92 const double si4 = sinc(mu_proj * cos_uu * tb); 93 117 94 // Expression in libCylinder.c (neither drC nor Vot are used) 118 tmp =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 ); // correct FF : square of sum of phase factors 119 95 const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 96 120 97 // To note also that in csparallelepiped.cpp, there is a function called 121 98 // cspkernel, which seems to make more sense and has the following comment: … … 126 103 // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 127 104 128 summj += Gauss76Wt[j] * tmp; 105 // correct FF : sum of square of phase factors 106 inner_total += Gauss76Wt[j] * form * form; 107 } 108 inner_total *= 0.5; 129 109 130 } 131 132 // value of the inner integral 133 answer = 0.5 * summj; 110 // now sum up the outer integral 111 const double si = sinc(mu * c_scaled * sigma); 112 outer_total += Gauss76Wt[i] * inner_total * si * si; 113 } 114 outer_total *= 0.5; 134 115 135 // finish the outer integral 136 double arg = 0.5 * mu* c_scaled *sigma; 137 if ( arg == 0.0 ) { 138 answer *= 1.0; 139 } else { 140 answer *= sin(arg)*sin(arg)/arg/arg; 141 } 142 143 // now sum up the outer integral 144 summ += Gauss76Wt[i] * answer; 145 146 } 147 148 answer = 0.5 * summ; 149 150 //convert from [1e-12 A-1] to [cm-1] 151 answer *= 1.0e-4; 152 153 return answer; 116 //convert from [1e-12 A-1] to [cm-1] 117 return 1.0e-4 * outer_total; 154 118 } 155 119 … … 170 134 double psi) 171 135 { 172 double tmp1, tmp2, tmp3, tmpt1, tmpt2, tmpt3; 136 double q, cos_val_a, cos_val_b, cos_val_c; 137 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 173 138 174 double q = sqrt(qx*qx+qy*qy); 175 double qx_scaled = qx/q; 176 double qy_scaled = qy/q; 139 // cspkernel in csparallelepiped recoded here 140 const double dr0 = core_sld-solvent_sld; 141 const double drA = arim_sld-solvent_sld; 142 const double drB = brim_sld-solvent_sld; 143 const double drC = crim_sld-solvent_sld; 177 144 178 // Convert angles given in degrees to radians 179 theta *= M_PI_180; 180 phi *= M_PI_180; 181 psi *= M_PI_180; 182 183 // Parallelepiped c axis orientation 184 double cparallel_x = cos(theta) * cos(phi); 185 double cparallel_y = sin(theta); 186 187 // Compute angle between q and parallelepiped axis 188 double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz; 189 190 // Parallelepiped a axis orientation 191 double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 192 double parallel_y = sin(psi)*cos(theta); 193 double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled; 194 195 // Parallelepiped b axis orientation 196 double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 197 double bparallel_y = cos(theta)*cos(psi); 198 double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled; 199 200 // The following tests should always pass 201 if (fabs(cos_val_c)>1.0) { 202 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 203 cos_val_c = 1.0; 204 } 205 if (fabs(cos_val_a)>1.0) { 206 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 207 cos_val_a = 1.0; 208 } 209 if (fabs(cos_val_b)>1.0) { 210 //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 211 cos_val_b = 1.0; 212 } 213 214 // cspkernel in csparallelepiped recoded here 215 double dr0 = core_sld-solvent_sld; 216 double drA = arim_sld-solvent_sld; 217 double drB = brim_sld-solvent_sld; 218 double drC = crim_sld-solvent_sld; 219 double Vin = length_a * length_b * length_c; 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; 220 149 // As for the 1D case, Vot is not used 221 150 //double Vot = (length_a * length_b * length_c + 222 151 // 2.0 * thick_rim_a * length_b * length_c + 223 152 // 2.0 * length_a * thick_rim_b * length_c + 224 153 // 2.0 * length_a * length_b * thick_rim_c); 225 double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 226 double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 227 double V3 = (2.0 * length_a * length_b * thick_rim_c); 154 228 155 // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 229 156 // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, 230 157 // but for the moment I let it like this until understanding better the code. 231 158 double ta = length_a + 2.0*thick_rim_a; 232 159 double tb = length_a + 2.0*thick_rim_b; 233 160 double tc = length_a + 2.0*thick_rim_c; 234 161 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 235 double argA = 0.5*q*length_a*cos_val_a;236 double argB = 0.5*q*length_b*cos_val_b;237 double argC = 0.5*q*length_c*cos_val_c;238 double argtA = 0.5*q*ta*cos_val_a;239 double argtB = 0.5*q*tb*cos_val_b;240 double argtC = 0.5*q*tc*cos_val_c;162 double siA = sinc(0.5*q*length_a*cos_val_a); 163 double siB = sinc(0.5*q*length_b*cos_val_b); 164 double siC = sinc(0.5*q*length_c*cos_val_c); 165 double siAt = sinc(0.5*q*ta*cos_val_a); 166 double siBt = sinc(0.5*q*tb*cos_val_b); 167 double siCt = sinc(0.5*q*tc*cos_val_c); 241 168 242 if(argA==0.0) {243 tmp1 = 1.0;244 } else {245 tmp1 = sin(argA)/argA;246 }247 if (argB==0.0) {248 tmp2 = 1.0;249 } else {250 tmp2 = sin(argB)/argB;251 }252 if (argC==0.0) {253 tmp3 = 1.0;254 } else {255 tmp3 = sin(argC)/argC;256 }257 if(argtA==0.0) {258 tmpt1 = 1.0;259 } else {260 tmpt1 = sin(argtA)/argtA;261 }262 if (argtB==0.0) {263 tmpt2 = 1.0;264 } else {265 tmpt2 = sin(argtB)/argtB;266 }267 if (argtC==0.0) {268 tmpt3 = 1.0;269 } else {270 tmpt3 = sin(argtC)*sin(argtC)/argtC/argtC;271 }272 169 273 170 // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 274 171 // in the 1D code, but should be checked! 275 double f = ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1 + 276 drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3); 172 double f = ( dr0*siA*siB*siC*Vin 173 + drA*(siAt-siA)*siB*siC*V1 174 + drB*siA*(siBt-siB)*siC*V2 175 + drC*siA*siB*(siCt*siCt-siC)*V3); 277 176 278 177 return 1.0e-4 * f * f;
Note: See TracChangeset
for help on using the changeset viewer.