[44bd2be] | 1 | double form_volume(double a_side, double b_side, double c_side, |
---|
| 2 | double arim_thickness, double brim_thickness, double crim_thickness); |
---|
| 3 | double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, |
---|
| 4 | double solvent_sld, double a_side, double b_side, double c_side, |
---|
| 5 | double arim_thickness, double brim_thickness, double crim_thickness); |
---|
| 6 | double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld, |
---|
| 7 | double crim_sld, double solvent_sld, double a_side, double b_side, |
---|
| 8 | double c_side, double arim_thickness, double brim_thickness, |
---|
| 9 | double crim_thickness, double theta, double phi, double psi); |
---|
| 10 | |
---|
| 11 | double form_volume(double a_side, double b_side, double c_side, |
---|
| 12 | double arim_thickness, double brim_thickness, double crim_thickness) |
---|
| 13 | { |
---|
| 14 | //return a_side * b_side * c_side; |
---|
| 15 | return a_side * b_side * c_side + |
---|
| 16 | 2.0 * arim_thickness * b_side * c_side + |
---|
| 17 | 2.0 * brim_thickness * a_side * c_side + |
---|
| 18 | 2.0 * crim_thickness * a_side * b_side; |
---|
| 19 | } |
---|
| 20 | |
---|
| 21 | double Iq(double q, |
---|
| 22 | double core_sld, |
---|
| 23 | double arim_sld, |
---|
| 24 | double brim_sld, |
---|
| 25 | double crim_sld, |
---|
| 26 | double solvent_sld, |
---|
| 27 | double a_side, |
---|
| 28 | double b_side, |
---|
| 29 | double c_side, |
---|
| 30 | double arim_thickness, |
---|
| 31 | double brim_thickness, |
---|
| 32 | double crim_thickness) |
---|
| 33 | { |
---|
| 34 | // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled |
---|
| 35 | // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) |
---|
| 36 | |
---|
| 37 | double t1, t2, t3, t4, tmp, answer; |
---|
| 38 | double mu = q * b_side; |
---|
| 39 | |
---|
| 40 | //calculate volume before rescaling (in original code, but not used) |
---|
| 41 | //double vol = form_volume(a_side, b_side, c_side, arim_thickness, brim_thickness, crim_thickness); |
---|
| 42 | //double vol = a_side * b_side * c_side + |
---|
| 43 | // 2.0 * arim_thickness * b_side * c_side + |
---|
| 44 | // 2.0 * brim_thickness * a_side * c_side + |
---|
| 45 | // 2.0 * crim_thickness * a_side * b_side; |
---|
| 46 | |
---|
| 47 | // Scale sides by B |
---|
| 48 | double a_scaled = a_side / b_side; |
---|
| 49 | double c_scaled = c_side / b_side; |
---|
[abdd01c] | 50 | |
---|
[44bd2be] | 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; |
---|
| 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 |
---|
| 63 | |
---|
| 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++) { |
---|
| 72 | |
---|
| 73 | double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); |
---|
| 74 | double mudum = mu * sqrt(1.0-sigma*sigma); |
---|
| 75 | |
---|
| 76 | double Vin = a_side * b_side * c_side; |
---|
| 77 | //double Vot = (a_side * b_side * c_side + |
---|
| 78 | // 2.0 * arim_thickness * b_side * c_side + |
---|
| 79 | // 2.0 * a_side * brim_thickness * c_side + |
---|
| 80 | // 2.0 * a_side * b_side * crim_thickness); |
---|
| 81 | double V1 = (2.0 * arim_thickness * b_side * c_side); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) |
---|
| 82 | double V2 = (2.0 * a_side * brim_thickness * c_side); // 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*arim_thickness)/b_side; |
---|
| 89 | double tb = (a_scaled+2.0*brim_thickness)/b_side; |
---|
| 90 | |
---|
| 91 | double arg1 = (0.5*mudum*a_scaled) * sin(0.5*M_PI*uu); |
---|
| 92 | double arg2 = (0.5*mudum) * cos(0.5*M_PI*uu); |
---|
| 93 | double arg3= (0.5*mudum*ta) * sin(0.5*M_PI*uu); |
---|
| 94 | double arg4= (0.5*mudum*tb) * cos(0.5*M_PI*uu); |
---|
| 95 | |
---|
| 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 | |
---|
| 117 | // 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 | |
---|
| 120 | // To note also that in csparallelepiped.cpp, there is a function called |
---|
| 121 | // cspkernel, which seems to make more sense and has the following comment: |
---|
| 122 | // This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010. |
---|
| 123 | // tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)* |
---|
| 124 | // ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3); // correct FF : square of sum of phase factors |
---|
| 125 | // This is the function called by csparallelepiped_analytical_2D_scaled, |
---|
| 126 | // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c |
---|
| 127 | |
---|
| 128 | summj += Gauss76Wt[j] * tmp; |
---|
| 129 | |
---|
| 130 | } |
---|
| 131 | |
---|
| 132 | // value of the inner integral |
---|
| 133 | answer = 0.5 * summj; |
---|
| 134 | |
---|
| 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; |
---|
| 154 | } |
---|
| 155 | |
---|
| 156 | double Iqxy(double qx, double qy, |
---|
| 157 | double core_sld, |
---|
| 158 | double arim_sld, |
---|
| 159 | double brim_sld, |
---|
| 160 | double crim_sld, |
---|
| 161 | double solvent_sld, |
---|
| 162 | double a_side, |
---|
| 163 | double b_side, |
---|
| 164 | double c_side, |
---|
| 165 | double arim_thickness, |
---|
| 166 | double brim_thickness, |
---|
| 167 | double crim_thickness, |
---|
| 168 | double theta, |
---|
| 169 | double phi, |
---|
| 170 | double psi) |
---|
| 171 | { |
---|
| 172 | double tmp1, tmp2, tmp3, tmpt1, tmpt2, tmpt3; |
---|
| 173 | |
---|
| 174 | double q = sqrt(qx*qx+qy*qy); |
---|
| 175 | double qx_scaled = qx/q; |
---|
| 176 | double qy_scaled = qy/q; |
---|
| 177 | |
---|
| 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 = a_side * b_side * c_side; |
---|
| 220 | // As for the 1D case, Vot is not used |
---|
| 221 | //double Vot = (a_side * b_side * c_side + |
---|
| 222 | // 2.0 * arim_thickness * b_side * c_side + |
---|
| 223 | // 2.0 * a_side * brim_thickness * c_side + |
---|
| 224 | // 2.0 * a_side * b_side * crim_thickness); |
---|
| 225 | double V1 = (2.0 * arim_thickness * b_side * c_side); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) |
---|
| 226 | double V2 = (2.0 * a_side * brim_thickness * c_side); // incorrect V2(aa*bb*cc+2*aa*tb*cc) |
---|
| 227 | double V3 = (2.0 * a_side * b_side * crim_thickness); |
---|
| 228 | // The definitions of ta, tb, tc are not the same as in the 1D case because there is no |
---|
| 229 | // the scaling by B. The use of a_side for the 3 of them seems clearly a mistake to me, |
---|
| 230 | // but for the moment I let it like this until understanding better the code. |
---|
| 231 | double ta = a_side + 2.0*arim_thickness; |
---|
| 232 | double tb = a_side + 2.0*brim_thickness; |
---|
| 233 | double tc = a_side + 2.0*crim_thickness; |
---|
| 234 | //handle arg=0 separately, as sin(t)/t -> 1 as t->0 |
---|
| 235 | double argA = 0.5*q*a_side*cos_val_a; |
---|
| 236 | double argB = 0.5*q*b_side*cos_val_b; |
---|
| 237 | double argC = 0.5*q*c_side*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; |
---|
| 241 | |
---|
| 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 | |
---|
| 273 | // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed |
---|
| 274 | // 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); |
---|
| 277 | |
---|
| 278 | return 1.0e-4 * f * f; |
---|
| 279 | } |
---|