[2222134] | 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); |
---|
[44bd2be] | 3 | double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, |
---|
[2222134] | 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); |
---|
[44bd2be] | 6 | double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld, |
---|
[2222134] | 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); |
---|
[44bd2be] | 10 | |
---|
[2222134] | 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) |
---|
[44bd2be] | 13 | { |
---|
[2222134] | 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; |
---|
[44bd2be] | 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, |
---|
[2222134] | 27 | double length_a, |
---|
| 28 | double length_b, |
---|
| 29 | double length_c, |
---|
| 30 | double thick_rim_a, |
---|
| 31 | double thick_rim_b, |
---|
| 32 | double thick_rim_c) |
---|
[44bd2be] | 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 | |
---|
[14838a3] | 37 | const double mu = 0.5 * q * length_b; |
---|
[44bd2be] | 38 | |
---|
| 39 | //calculate volume before rescaling (in original code, but not used) |
---|
[14838a3] | 40 | //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); |
---|
[2222134] | 41 | //double vol = length_a * length_b * length_c + |
---|
| 42 | // 2.0 * thick_rim_a * length_b * length_c + |
---|
| 43 | // 2.0 * thick_rim_b * length_a * length_c + |
---|
| 44 | // 2.0 * thick_rim_c * length_a * length_b; |
---|
[44bd2be] | 45 | |
---|
| 46 | // Scale sides by B |
---|
[14838a3] | 47 | const double a_scaled = length_a / length_b; |
---|
| 48 | const double c_scaled = length_c / length_b; |
---|
| 49 | |
---|
| 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 | |
---|
| 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) |
---|
| 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 drC_Vot = (crim_sld-solvent_sld)*Vot; |
---|
| 70 | |
---|
| 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 | |
---|
| 76 | // outer integral (with gauss points), integration limits = 0, 1 |
---|
| 77 | double outer_total = 0; //initialize integral |
---|
| 78 | |
---|
| 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); |
---|
[1e7b0db0] | 89 | const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); |
---|
| 90 | const double si2 = sas_sinx_x(mu_proj * cos_uu); |
---|
| 91 | const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); |
---|
| 92 | const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); |
---|
[abdd01c] | 93 | |
---|
[44bd2be] | 94 | // Expression in libCylinder.c (neither drC nor Vot are used) |
---|
[14838a3] | 95 | const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; |
---|
| 96 | |
---|
[44bd2be] | 97 | // To note also that in csparallelepiped.cpp, there is a function called |
---|
| 98 | // cspkernel, which seems to make more sense and has the following comment: |
---|
| 99 | // This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010. |
---|
| 100 | // tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)* |
---|
| 101 | // ( 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 |
---|
| 102 | // This is the function called by csparallelepiped_analytical_2D_scaled, |
---|
| 103 | // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c |
---|
| 104 | |
---|
[14838a3] | 105 | // correct FF : sum of square of phase factors |
---|
| 106 | inner_total += Gauss76Wt[j] * form * form; |
---|
[44bd2be] | 107 | } |
---|
[14838a3] | 108 | inner_total *= 0.5; |
---|
[44bd2be] | 109 | |
---|
[14838a3] | 110 | // now sum up the outer integral |
---|
[1e7b0db0] | 111 | const double si = sas_sinx_x(mu * c_scaled * sigma); |
---|
[14838a3] | 112 | outer_total += Gauss76Wt[i] * inner_total * si * si; |
---|
[44bd2be] | 113 | } |
---|
[14838a3] | 114 | outer_total *= 0.5; |
---|
[44bd2be] | 115 | |
---|
[14838a3] | 116 | //convert from [1e-12 A-1] to [cm-1] |
---|
| 117 | return 1.0e-4 * outer_total; |
---|
[44bd2be] | 118 | } |
---|
| 119 | |
---|
| 120 | double Iqxy(double qx, double qy, |
---|
| 121 | double core_sld, |
---|
| 122 | double arim_sld, |
---|
| 123 | double brim_sld, |
---|
| 124 | double crim_sld, |
---|
| 125 | double solvent_sld, |
---|
[2222134] | 126 | double length_a, |
---|
| 127 | double length_b, |
---|
| 128 | double length_c, |
---|
| 129 | double thick_rim_a, |
---|
| 130 | double thick_rim_b, |
---|
| 131 | double thick_rim_c, |
---|
[44bd2be] | 132 | double theta, |
---|
| 133 | double phi, |
---|
| 134 | double psi) |
---|
| 135 | { |
---|
[92dfe0c] | 136 | double q, zhat, yhat, xhat; |
---|
| 137 | ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); |
---|
[14838a3] | 138 | |
---|
| 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; |
---|
| 144 | |
---|
| 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; |
---|
[44bd2be] | 149 | // As for the 1D case, Vot is not used |
---|
[14838a3] | 150 | //double Vot = (length_a * length_b * length_c + |
---|
[2222134] | 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); |
---|
[14838a3] | 154 | |
---|
[44bd2be] | 155 | // The definitions of ta, tb, tc are not the same as in the 1D case because there is no |
---|
[2222134] | 156 | // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, |
---|
[44bd2be] | 157 | // but for the moment I let it like this until understanding better the code. |
---|
[14838a3] | 158 | double ta = length_a + 2.0*thick_rim_a; |
---|
[2222134] | 159 | double tb = length_a + 2.0*thick_rim_b; |
---|
| 160 | double tc = length_a + 2.0*thick_rim_c; |
---|
[44bd2be] | 161 | //handle arg=0 separately, as sin(t)/t -> 1 as t->0 |
---|
[92dfe0c] | 162 | double siA = sas_sinx_x(0.5*q*length_a*xhat); |
---|
| 163 | double siB = sas_sinx_x(0.5*q*length_b*yhat); |
---|
| 164 | double siC = sas_sinx_x(0.5*q*length_c*zhat); |
---|
| 165 | double siAt = sas_sinx_x(0.5*q*ta*xhat); |
---|
| 166 | double siBt = sas_sinx_x(0.5*q*tb*yhat); |
---|
| 167 | double siCt = sas_sinx_x(0.5*q*tc*zhat); |
---|
[44bd2be] | 168 | |
---|
| 169 | |
---|
| 170 | // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed |
---|
| 171 | // in the 1D code, but should be checked! |
---|
[14838a3] | 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); |
---|
[44bd2be] | 176 | |
---|
| 177 | return 1.0e-4 * f * f; |
---|
| 178 | } |
---|