[f7930be] | 1 | FourShell(double dp[], double q) |
---|
| 2 | { |
---|
| 3 | // variables are: |
---|
| 4 | //[0] scale factor |
---|
| 5 | //[1] radius of core [ᅵ] |
---|
| 6 | //[2] SLD of the core [ᅵ-2] |
---|
| 7 | //[3] thickness of shell 1 [ᅵ] |
---|
| 8 | //[4] SLD of shell 1 |
---|
| 9 | //[5] thickness of shell 2 [ᅵ] |
---|
| 10 | //[6] SLD of shell 2 |
---|
| 11 | //[7] thickness of shell 3 |
---|
| 12 | //[8] SLD of shell 3 |
---|
| 13 | //[9] thickness of shell 3 |
---|
| 14 | //[10] SLD of shell 3 |
---|
| 15 | //[11] SLD of solvent |
---|
| 16 | //[12] background [cm-1] |
---|
| 17 | |
---|
| 18 | double x,pi; |
---|
| 19 | double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg; //my local names |
---|
| 20 | double bes,f,vol,qr,contr,f2; |
---|
| 21 | double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4; |
---|
| 22 | |
---|
| 23 | pi = 4.0*atan(1.0); |
---|
| 24 | x=q; |
---|
| 25 | |
---|
| 26 | scale = dp[0]; |
---|
| 27 | rcore = dp[1]; |
---|
| 28 | rhocore = dp[2]; |
---|
| 29 | thick1 = dp[3]; |
---|
| 30 | rhoshel1 = dp[4]; |
---|
| 31 | thick2 = dp[5]; |
---|
| 32 | rhoshel2 = dp[6]; |
---|
| 33 | thick3 = dp[7]; |
---|
| 34 | rhoshel3 = dp[8]; |
---|
| 35 | thick4 = dp[9]; |
---|
| 36 | rhoshel4 = dp[10]; |
---|
| 37 | rhosolv = dp[11]; |
---|
| 38 | bkg = dp[12]; |
---|
| 39 | |
---|
| 40 | // core first, then add in shells |
---|
| 41 | qr=x*rcore; |
---|
| 42 | contr = rhocore-rhoshel1; |
---|
| 43 | if(qr == 0){ |
---|
| 44 | bes = 1.0; |
---|
| 45 | }else{ |
---|
| 46 | bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); |
---|
| 47 | } |
---|
| 48 | vol = 4.0*pi/3.0*rcore*rcore*rcore; |
---|
| 49 | f = vol*bes*contr; |
---|
| 50 | //now the shell (1) |
---|
| 51 | qr=x*(rcore+thick1); |
---|
| 52 | contr = rhoshel1-rhoshel2; |
---|
| 53 | if(qr == 0){ |
---|
| 54 | bes = 1.0; |
---|
| 55 | }else{ |
---|
| 56 | bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); |
---|
| 57 | } |
---|
| 58 | vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); |
---|
| 59 | f += vol*bes*contr; |
---|
| 60 | //now the shell (2) |
---|
| 61 | qr=x*(rcore+thick1+thick2); |
---|
| 62 | contr = rhoshel2-rhoshel3; |
---|
| 63 | if(qr == 0){ |
---|
| 64 | bes = 1.0; |
---|
| 65 | }else{ |
---|
| 66 | bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); |
---|
| 67 | } |
---|
| 68 | vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); |
---|
| 69 | f += vol*bes*contr; |
---|
| 70 | //now the shell (3) |
---|
| 71 | qr=x*(rcore+thick1+thick2+thick3); |
---|
| 72 | contr = rhoshel3-rhoshel4; |
---|
| 73 | if(qr == 0){ |
---|
| 74 | bes = 1.0; |
---|
| 75 | }else{ |
---|
| 76 | bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); |
---|
| 77 | } |
---|
| 78 | vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); |
---|
| 79 | f += vol*bes*contr; |
---|
| 80 | //now the shell (4) |
---|
| 81 | qr=x*(rcore+thick1+thick2+thick3+thick4); |
---|
| 82 | contr = rhoshel4-rhosolv; |
---|
| 83 | if(qr == 0){ |
---|
| 84 | bes = 1.0; |
---|
| 85 | }else{ |
---|
| 86 | bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); |
---|
| 87 | } |
---|
| 88 | vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4); |
---|
| 89 | f += vol*bes*contr; |
---|
| 90 | |
---|
| 91 | |
---|
| 92 | // normalize to particle volume and rescale from [ᅵ-1] to [cm-1] |
---|
| 93 | f2 = f*f/vol*1.0e8; |
---|
| 94 | |
---|
| 95 | //scale if desired |
---|
| 96 | f2 *= scale; |
---|
| 97 | // then add in the background |
---|
| 98 | f2 += bkg; |
---|
| 99 | |
---|
| 100 | return(f2); |
---|
| 101 | } |
---|
| 102 | |
---|
| 103 | // above is cut and past with no changes from libigor four shell |
---|
| 104 | static double |
---|
| 105 | f_exp(double q, double r, double sld_in, double sld_out, |
---|
| 106 | double thickness, double A) |
---|
| 107 | { |
---|
| 108 | const double vol = 4.0/3.0 * M_PI * r * r * r; |
---|
| 109 | const double qr = q * r; |
---|
| 110 | const double alpha = A * r/thickness; |
---|
| 111 | const double bes = sph_j1c(qr); |
---|
| 112 | const double B = (sld_out - sld_in)/expm1(A); |
---|
| 113 | const double C = sld_in - B; |
---|
| 114 | double fun; |
---|
| 115 | if (qr == 0.0) { |
---|
| 116 | fun = 1.0; |
---|
| 117 | } else if (fabs(A) > 0.0) { |
---|
| 118 | const double qrsq = qr*qr; |
---|
| 119 | const double alphasq = alpha*alpha; |
---|
| 120 | const double sumsq = alphasq + qrsq; |
---|
| 121 | double sinqr, cosqr; |
---|
| 122 | SINCOS(qr, sinqr, cosqr); |
---|
| 123 | fun = -3.0( |
---|
| 124 | ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq |
---|
| 125 | - (alpha*sinqr/qr - cosqr) |
---|
| 126 | ) / sumsq; |
---|
| 127 | } else { |
---|
| 128 | fun = bes; |
---|
| 129 | } |
---|
| 130 | return vol * (B*fun + C*bes); |
---|
| 131 | } |
---|
| 132 | |
---|
| 133 | static double |
---|
| 134 | f_linear(double q, double r, double sld, double slope) |
---|
| 135 | { |
---|
| 136 | const double vol = 4.0/3.0 * M_PI * r * r * r; |
---|
| 137 | const double qr = q * r; |
---|
| 138 | const double bes = sph_j1c(qr); |
---|
| 139 | double fun = 0.0; |
---|
| 140 | if (qr > 0.0) { |
---|
| 141 | const double qrsq = qr*qr; |
---|
| 142 | double sinqr, cosqr; |
---|
| 143 | SINCOS(qr, sinqr, cosqr); |
---|
| 144 | // Jae-He's code seems to simplify to this |
---|
| 145 | // fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq); |
---|
| 146 | // Rederiving the math, we get the following instead: |
---|
| 147 | fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq); |
---|
| 148 | } |
---|
| 149 | return vol * (sld*bes + fun); |
---|
| 150 | } |
---|
| 151 | |
---|
| 152 | static double |
---|
| 153 | f_constant(double q, double r, double sld) |
---|
| 154 | { |
---|
| 155 | const double bes = sph_j1c(q * r); |
---|
| 156 | const double vol = 4.0/3.0 * M_PI * r * r * r; |
---|
| 157 | return sld * vol * bes; |
---|
| 158 | } |
---|
| 159 | |
---|
| 160 | static double |
---|
| 161 | form_volume(double core_radius, double n, double thickness[]) |
---|
| 162 | { |
---|
| 163 | int i; |
---|
| 164 | double r = core_radius; |
---|
| 165 | for (i=0; i < n; i++) { |
---|
| 166 | r += thickness[i]; |
---|
| 167 | } |
---|
| 168 | return 4.0/3.0 * M_PI * r * r * r; |
---|
| 169 | } |
---|
| 170 | |
---|
| 171 | static double |
---|
| 172 | Iq(double q, double core_sld, double core_radius, double solvent_sld, |
---|
| 173 | double n, double in_sld[], double out_sld[], double thickness[], |
---|
| 174 | double A[]) |
---|
| 175 | { |
---|
| 176 | int i; |
---|
| 177 | r = core_radius; |
---|
| 178 | f = f_constant(q, r, core_sld); |
---|
| 179 | for (i=0; i<n; i++){ |
---|
| 180 | const double r0 = r; |
---|
| 181 | r += thickness[i]; |
---|
| 182 | if (r == r0) { |
---|
| 183 | // no thickness, so nothing to add |
---|
[abdd01c] | 184 | } else if (fabs(A[i]) < 1.0e-16 || sld_out[i] == sld_in[i]) { |
---|
[f7930be] | 185 | f -= f_constant(q, r0, sld_in[i]); |
---|
| 186 | f += f_constant(q, r, sld_in[i]); |
---|
[abdd01c] | 187 | } else if (fabs(A[i]) < 1.0e-4) { |
---|
[f7930be] | 188 | const double slope = (sld_out[i] - sld_in[i])/thickness[i]; |
---|
| 189 | f -= f_linear(q, r0, sld_in[i], slope); |
---|
| 190 | f += f_linear(q, r, sld_out[i], slope); |
---|
| 191 | } else { |
---|
| 192 | f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]); |
---|
| 193 | f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]); |
---|
| 194 | } |
---|
| 195 | } |
---|
| 196 | f -= f_constant(q, r, solvent_sld); |
---|
| 197 | f2 = f * f * 1.0e-4; |
---|
| 198 | |
---|
| 199 | return f2; |
---|
| 200 | } |
---|