[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 | { |
---|
[2a0b2b1] | 5 | return length_a * length_b * length_c + |
---|
| 6 | 2.0 * thick_rim_a * length_b * length_c + |
---|
[2222134] | 7 | 2.0 * thick_rim_b * length_a * length_c + |
---|
| 8 | 2.0 * thick_rim_c * length_a * length_b; |
---|
[44bd2be] | 9 | } |
---|
| 10 | |
---|
[becded3] | 11 | static double |
---|
| 12 | Iq(double q, |
---|
[44bd2be] | 13 | double core_sld, |
---|
| 14 | double arim_sld, |
---|
| 15 | double brim_sld, |
---|
| 16 | double crim_sld, |
---|
| 17 | double solvent_sld, |
---|
[2222134] | 18 | double length_a, |
---|
| 19 | double length_b, |
---|
| 20 | double length_c, |
---|
| 21 | double thick_rim_a, |
---|
| 22 | double thick_rim_b, |
---|
| 23 | double thick_rim_c) |
---|
[44bd2be] | 24 | { |
---|
[3a1fc7d] | 25 | // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c |
---|
[44bd2be] | 26 | // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) |
---|
[c69d6d6] | 27 | //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) |
---|
[2a0b2b1] | 28 | |
---|
[14838a3] | 29 | const double mu = 0.5 * q * length_b; |
---|
[2a0b2b1] | 30 | |
---|
[44bd2be] | 31 | // Scale sides by B |
---|
[3a1fc7d] | 32 | const double a_over_b = length_a / length_b; |
---|
| 33 | const double c_over_b = length_c / length_b; |
---|
[14838a3] | 34 | |
---|
[3a1fc7d] | 35 | double tA_over_b = a_over_b + 2.0*thick_rim_a/length_b; |
---|
| 36 | double tB_over_b = 1+ 2.0*thick_rim_b/length_b; |
---|
| 37 | double tC_over_b = c_over_b + 2.0*thick_rim_c/length_b; |
---|
[14838a3] | 38 | |
---|
| 39 | double Vin = length_a * length_b * length_c; |
---|
[3a1fc7d] | 40 | double VtA = (2.0 * thick_rim_a * length_b * length_c); |
---|
| 41 | double VtB = (2.0 * length_a * thick_rim_b * length_c); |
---|
| 42 | double VtC = (2.0 * length_a * length_b * thick_rim_c); |
---|
[14838a3] | 43 | |
---|
| 44 | // Scale factors (note that drC is not used later) |
---|
[3a1fc7d] | 45 | const double dr0 = (core_sld-solvent_sld); |
---|
| 46 | const double drA = (arim_sld-solvent_sld); |
---|
| 47 | const double drB = (brim_sld-solvent_sld); |
---|
| 48 | const double drC = (crim_sld-solvent_sld); |
---|
[14838a3] | 49 | |
---|
| 50 | // Precompute scale factors for combining cross terms from the shape |
---|
[3a1fc7d] | 51 | const double dr0_Vin = dr0*Vin; |
---|
| 52 | const double drA_VtA = drA*VtA; |
---|
| 53 | const double drB_VtB = drB*VtB; |
---|
| 54 | const double drC_VtC = drC*VtC; |
---|
| 55 | const double drV_delta = dr0_Vin - drA_VtA - drB_VtB - drC_VtC; |
---|
| 56 | |
---|
| 57 | /* *************** algorithm description ****************** |
---|
| 58 | |
---|
| 59 | // Rewrite f as x*siC + y*siCt to move the siC/siCt calculation out |
---|
| 60 | // of the inner loop. That is: |
---|
| 61 | |
---|
| 62 | f = (di-da-db-dc) sa sb sc + da sa' sb sc + db sa sb' sc + dc sa sb sc' |
---|
| 63 | = [ (di-da-db-dc) sa sb + da sa' sb + db sa sb' ] sc + [dc sa sb] sc' |
---|
| 64 | = x sc + y sc' |
---|
| 65 | |
---|
| 66 | // where: |
---|
| 67 | di = delta rho_core V_core |
---|
| 68 | da = delta rho_rimA V_rimA |
---|
| 69 | db = delta rho_rimB V_rimB |
---|
| 70 | dc = delta rho_rimC V_rimC |
---|
| 71 | sa = j_0 (q_a a/2) // siA the code |
---|
| 72 | sb = j_0 (q_b b/2) |
---|
| 73 | sc = j_0 (q_c c/2) |
---|
| 74 | sa' = j_0(q_a a_rim/2) // siAt the code |
---|
| 75 | sb' = j_0(q_b b_rim/2) |
---|
| 76 | sc' = j_0(q_c c_rim/2) |
---|
| 77 | |
---|
| 78 | // qa, qb, and qc are generated using polar coordinates, with the |
---|
| 79 | // outer loop integrating over [0,1] after the u-substitution |
---|
| 80 | // sigma = cos(theta), sqrt(1-sigma^2) = sin(theta) |
---|
| 81 | // and inner loop integrating over [0, pi/2] as |
---|
| 82 | // uu = phi |
---|
| 83 | |
---|
| 84 | ************************************************************ */ |
---|
[14838a3] | 85 | |
---|
| 86 | // outer integral (with gauss points), integration limits = 0, 1 |
---|
[3a1fc7d] | 87 | double outer_sum = 0; //initialize integral |
---|
[14838a3] | 88 | for( int i=0; i<76; i++) { |
---|
| 89 | double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); |
---|
| 90 | double mu_proj = mu * sqrt(1.0-sigma*sigma); |
---|
| 91 | |
---|
[3a1fc7d] | 92 | // inner integral (with gauss points), integration limits = 0, pi/2 |
---|
| 93 | const double siC = sas_sinx_x(mu * sigma * c_over_b); |
---|
| 94 | const double siCt = sas_sinx_x(mu * sigma * tC_over_b); |
---|
| 95 | double inner_sum = 0.0; |
---|
[14838a3] | 96 | for(int j=0; j<76; j++) { |
---|
| 97 | const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); |
---|
| 98 | double sin_uu, cos_uu; |
---|
| 99 | SINCOS(M_PI_2*uu, sin_uu, cos_uu); |
---|
[3a1fc7d] | 100 | const double siA = sas_sinx_x(mu_proj * sin_uu * a_over_b); |
---|
| 101 | const double siB = sas_sinx_x(mu_proj * cos_uu ); |
---|
| 102 | const double siAt = sas_sinx_x(mu_proj * sin_uu * tA_over_b); |
---|
| 103 | const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b); |
---|
| 104 | |
---|
| 105 | const double x = drV_delta*siA*siB + drA_VtA*siB*siAt + drB_VtB*siA*siBt; |
---|
| 106 | const double form = x*siC + drC_VtC*siA*siB*siCt; |
---|
| 107 | |
---|
| 108 | inner_sum += Gauss76Wt[j] * form * form; |
---|
[44bd2be] | 109 | } |
---|
[3a1fc7d] | 110 | inner_sum *= 0.5; |
---|
[14838a3] | 111 | // now sum up the outer integral |
---|
[3a1fc7d] | 112 | outer_sum += Gauss76Wt[i] * inner_sum; |
---|
[44bd2be] | 113 | } |
---|
[3a1fc7d] | 114 | outer_sum *= 0.5; |
---|
[44bd2be] | 115 | |
---|
[14838a3] | 116 | //convert from [1e-12 A-1] to [cm-1] |
---|
[3a1fc7d] | 117 | return 1.0e-4 * outer_sum; |
---|
[44bd2be] | 118 | } |
---|
| 119 | |
---|
[becded3] | 120 | static double |
---|
| 121 | Iqxy(double qa, double qb, double qc, |
---|
[44bd2be] | 122 | double core_sld, |
---|
| 123 | double arim_sld, |
---|
| 124 | double brim_sld, |
---|
| 125 | double crim_sld, |
---|
| 126 | double solvent_sld, |
---|
[2222134] | 127 | double length_a, |
---|
| 128 | double length_b, |
---|
| 129 | double length_c, |
---|
| 130 | double thick_rim_a, |
---|
| 131 | double thick_rim_b, |
---|
[becded3] | 132 | double thick_rim_c) |
---|
[44bd2be] | 133 | { |
---|
[14838a3] | 134 | // cspkernel in csparallelepiped recoded here |
---|
| 135 | const double dr0 = core_sld-solvent_sld; |
---|
| 136 | const double drA = arim_sld-solvent_sld; |
---|
| 137 | const double drB = brim_sld-solvent_sld; |
---|
| 138 | const double drC = crim_sld-solvent_sld; |
---|
| 139 | |
---|
| 140 | double Vin = length_a * length_b * length_c; |
---|
[3a1fc7d] | 141 | double VtA = 2.0 * thick_rim_a * length_b * length_c; |
---|
| 142 | double VtB = 2.0 * length_a * thick_rim_b * length_c; |
---|
| 143 | double VtC = 2.0 * length_a * length_b * thick_rim_c; |
---|
[14838a3] | 144 | |
---|
[44bd2be] | 145 | // The definitions of ta, tb, tc are not the same as in the 1D case because there is no |
---|
[c69d6d6] | 146 | // the scaling by B. |
---|
[3a1fc7d] | 147 | double tA = length_a + 2.0*thick_rim_a; |
---|
| 148 | double tB = length_b + 2.0*thick_rim_b; |
---|
| 149 | double tC = length_c + 2.0*thick_rim_c; |
---|
[44bd2be] | 150 | //handle arg=0 separately, as sin(t)/t -> 1 as t->0 |
---|
[2a0b2b1] | 151 | double siA = sas_sinx_x(0.5*length_a*qa); |
---|
| 152 | double siB = sas_sinx_x(0.5*length_b*qb); |
---|
| 153 | double siC = sas_sinx_x(0.5*length_c*qc); |
---|
[3a1fc7d] | 154 | double siAt = sas_sinx_x(0.5*tA*qa); |
---|
| 155 | double siBt = sas_sinx_x(0.5*tB*qb); |
---|
| 156 | double siCt = sas_sinx_x(0.5*tC*qc); |
---|
[2a0b2b1] | 157 | |
---|
[44bd2be] | 158 | |
---|
| 159 | // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed |
---|
| 160 | // in the 1D code, but should be checked! |
---|
[3a1fc7d] | 161 | double f = ( dr0*Vin*siA*siB*siC |
---|
| 162 | + drA*VtA*(siAt-siA)*siB*siC |
---|
| 163 | + drB*VtB*siA*(siBt-siB)*siC |
---|
| 164 | + drC*VtC*siA*siB*(siCt-siC)); |
---|
[2a0b2b1] | 165 | |
---|
[44bd2be] | 166 | return 1.0e-4 * f * f; |
---|
| 167 | } |
---|