[becded3] | 1 | static double |
---|
| 2 | form_volume(double length_a, double length_b, double length_c, |
---|
| 3 | double thick_rim_a, double thick_rim_b, double thick_rim_c) |
---|
[44bd2be] | 4 | { |
---|
[2222134] | 5 | //return length_a * length_b * length_c; |
---|
[2a0b2b1] | 6 | return length_a * length_b * length_c + |
---|
| 7 | 2.0 * thick_rim_a * length_b * length_c + |
---|
[2222134] | 8 | 2.0 * thick_rim_b * length_a * length_c + |
---|
| 9 | 2.0 * thick_rim_c * length_a * length_b; |
---|
[44bd2be] | 10 | } |
---|
| 11 | |
---|
[becded3] | 12 | static double |
---|
| 13 | Iq(double q, |
---|
[44bd2be] | 14 | double core_sld, |
---|
| 15 | double arim_sld, |
---|
| 16 | double brim_sld, |
---|
| 17 | double crim_sld, |
---|
| 18 | double solvent_sld, |
---|
[2222134] | 19 | double length_a, |
---|
| 20 | double length_b, |
---|
| 21 | double length_c, |
---|
| 22 | double thick_rim_a, |
---|
| 23 | double thick_rim_b, |
---|
| 24 | double thick_rim_c) |
---|
[44bd2be] | 25 | { |
---|
| 26 | // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled |
---|
| 27 | // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) |
---|
[c69d6d6] | 28 | //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) |
---|
[2a0b2b1] | 29 | |
---|
[14838a3] | 30 | const double mu = 0.5 * q * length_b; |
---|
[2a0b2b1] | 31 | |
---|
[44bd2be] | 32 | //calculate volume before rescaling (in original code, but not used) |
---|
[2a0b2b1] | 33 | //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); |
---|
| 34 | //double vol = length_a * length_b * length_c + |
---|
| 35 | // 2.0 * thick_rim_a * length_b * length_c + |
---|
[2222134] | 36 | // 2.0 * thick_rim_b * length_a * length_c + |
---|
| 37 | // 2.0 * thick_rim_c * length_a * length_b; |
---|
[2a0b2b1] | 38 | |
---|
[44bd2be] | 39 | // Scale sides by B |
---|
[14838a3] | 40 | const double a_scaled = length_a / length_b; |
---|
| 41 | const double c_scaled = length_c / length_b; |
---|
| 42 | |
---|
[c69d6d6] | 43 | double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b; |
---|
| 44 | double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b; |
---|
| 45 | double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present |
---|
[14838a3] | 46 | |
---|
| 47 | double Vin = length_a * length_b * length_c; |
---|
| 48 | //double Vot = (length_a * length_b * length_c + |
---|
| 49 | // 2.0 * thick_rim_a * length_b * length_c + |
---|
| 50 | // 2.0 * length_a * thick_rim_b * length_c + |
---|
| 51 | // 2.0 * length_a * length_b * thick_rim_c); |
---|
| 52 | double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) |
---|
| 53 | double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) |
---|
[c69d6d6] | 54 | double V3 = (2.0 * length_a * length_b * thick_rim_c); //not present |
---|
[14838a3] | 55 | |
---|
| 56 | // Scale factors (note that drC is not used later) |
---|
| 57 | const double drho0 = (core_sld-solvent_sld); |
---|
| 58 | const double drhoA = (arim_sld-solvent_sld); |
---|
| 59 | const double drhoB = (brim_sld-solvent_sld); |
---|
[c69d6d6] | 60 | const double drhoC = (crim_sld-solvent_sld); // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; |
---|
| 61 | |
---|
[14838a3] | 62 | |
---|
| 63 | // Precompute scale factors for combining cross terms from the shape |
---|
| 64 | const double scale23 = drhoA*V1; |
---|
| 65 | const double scale14 = drhoB*V2; |
---|
[c69d6d6] | 66 | const double scale24 = drhoC*V3; |
---|
| 67 | const double scale11 = drho0*Vin; |
---|
| 68 | const double scale12 = drho0*Vin - scale23 - scale14 - scale24; |
---|
[14838a3] | 69 | |
---|
| 70 | // outer integral (with gauss points), integration limits = 0, 1 |
---|
| 71 | double outer_total = 0; //initialize integral |
---|
| 72 | |
---|
| 73 | for( int i=0; i<76; i++) { |
---|
| 74 | double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); |
---|
| 75 | double mu_proj = mu * sqrt(1.0-sigma*sigma); |
---|
| 76 | |
---|
| 77 | // inner integral (with gauss points), integration limits = 0, 1 |
---|
| 78 | double inner_total = 0.0; |
---|
[c69d6d6] | 79 | double inner_total_crim = 0.0; |
---|
[14838a3] | 80 | for(int j=0; j<76; j++) { |
---|
| 81 | const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); |
---|
| 82 | double sin_uu, cos_uu; |
---|
| 83 | SINCOS(M_PI_2*uu, sin_uu, cos_uu); |
---|
[1e7b0db0] | 84 | const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); |
---|
[c69d6d6] | 85 | const double si2 = sas_sinx_x(mu_proj * cos_uu ); |
---|
[1e7b0db0] | 86 | const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); |
---|
| 87 | const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); |
---|
[abdd01c] | 88 | |
---|
[44bd2be] | 89 | // Expression in libCylinder.c (neither drC nor Vot are used) |
---|
[14838a3] | 90 | const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; |
---|
[c69d6d6] | 91 | const double form_crim = scale11*si1*si2; |
---|
[2a0b2b1] | 92 | |
---|
[14838a3] | 93 | // correct FF : sum of square of phase factors |
---|
| 94 | inner_total += Gauss76Wt[j] * form * form; |
---|
[c69d6d6] | 95 | inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; |
---|
[44bd2be] | 96 | } |
---|
[14838a3] | 97 | inner_total *= 0.5; |
---|
[c69d6d6] | 98 | inner_total_crim *= 0.5; |
---|
[14838a3] | 99 | // now sum up the outer integral |
---|
[1e7b0db0] | 100 | const double si = sas_sinx_x(mu * c_scaled * sigma); |
---|
[c69d6d6] | 101 | const double si_crim = sas_sinx_x(mu * tc * sigma); |
---|
| 102 | outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim); |
---|
[44bd2be] | 103 | } |
---|
[14838a3] | 104 | outer_total *= 0.5; |
---|
[44bd2be] | 105 | |
---|
[14838a3] | 106 | //convert from [1e-12 A-1] to [cm-1] |
---|
| 107 | return 1.0e-4 * outer_total; |
---|
[44bd2be] | 108 | } |
---|
| 109 | |
---|
[becded3] | 110 | static double |
---|
| 111 | Iqxy(double qa, double qb, double qc, |
---|
[44bd2be] | 112 | double core_sld, |
---|
| 113 | double arim_sld, |
---|
| 114 | double brim_sld, |
---|
| 115 | double crim_sld, |
---|
| 116 | double solvent_sld, |
---|
[2222134] | 117 | double length_a, |
---|
| 118 | double length_b, |
---|
| 119 | double length_c, |
---|
| 120 | double thick_rim_a, |
---|
| 121 | double thick_rim_b, |
---|
[becded3] | 122 | double thick_rim_c) |
---|
[44bd2be] | 123 | { |
---|
[14838a3] | 124 | // cspkernel in csparallelepiped recoded here |
---|
| 125 | const double dr0 = core_sld-solvent_sld; |
---|
| 126 | const double drA = arim_sld-solvent_sld; |
---|
| 127 | const double drB = brim_sld-solvent_sld; |
---|
| 128 | const double drC = crim_sld-solvent_sld; |
---|
| 129 | |
---|
| 130 | double Vin = length_a * length_b * length_c; |
---|
| 131 | double V1 = 2.0 * thick_rim_a * length_b * length_c; // incorrect V1(aa*bb*cc+2*ta*bb*cc) |
---|
| 132 | double V2 = 2.0 * length_a * thick_rim_b * length_c; // incorrect V2(aa*bb*cc+2*aa*tb*cc) |
---|
| 133 | double V3 = 2.0 * length_a * length_b * thick_rim_c; |
---|
[44bd2be] | 134 | // As for the 1D case, Vot is not used |
---|
[14838a3] | 135 | //double Vot = (length_a * length_b * length_c + |
---|
[2222134] | 136 | // 2.0 * thick_rim_a * length_b * length_c + |
---|
| 137 | // 2.0 * length_a * thick_rim_b * length_c + |
---|
| 138 | // 2.0 * length_a * length_b * thick_rim_c); |
---|
[14838a3] | 139 | |
---|
[44bd2be] | 140 | // The definitions of ta, tb, tc are not the same as in the 1D case because there is no |
---|
[c69d6d6] | 141 | // the scaling by B. |
---|
[14838a3] | 142 | double ta = length_a + 2.0*thick_rim_a; |
---|
[c69d6d6] | 143 | double tb = length_b + 2.0*thick_rim_b; |
---|
| 144 | double tc = length_c + 2.0*thick_rim_c; |
---|
[44bd2be] | 145 | //handle arg=0 separately, as sin(t)/t -> 1 as t->0 |
---|
[2a0b2b1] | 146 | double siA = sas_sinx_x(0.5*length_a*qa); |
---|
| 147 | double siB = sas_sinx_x(0.5*length_b*qb); |
---|
| 148 | double siC = sas_sinx_x(0.5*length_c*qc); |
---|
| 149 | double siAt = sas_sinx_x(0.5*ta*qa); |
---|
| 150 | double siBt = sas_sinx_x(0.5*tb*qb); |
---|
| 151 | double siCt = sas_sinx_x(0.5*tc*qc); |
---|
| 152 | |
---|
[44bd2be] | 153 | |
---|
| 154 | // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed |
---|
| 155 | // in the 1D code, but should be checked! |
---|
[14838a3] | 156 | double f = ( dr0*siA*siB*siC*Vin |
---|
| 157 | + drA*(siAt-siA)*siB*siC*V1 |
---|
| 158 | + drB*siA*(siBt-siB)*siC*V2 |
---|
[c69d6d6] | 159 | + drC*siA*siB*(siCt-siC)*V3); |
---|
[2a0b2b1] | 160 | |
---|
[44bd2be] | 161 | return 1.0e-4 * f * f; |
---|
| 162 | } |
---|