[82c299f] | 1 | double Iq(double q, double case_num, |
---|
| 2 | double Na, double Phia, double va, double a_sld, double ba, |
---|
| 3 | double Nb, double Phib, double vb, double b_sld, double bb, |
---|
| 4 | double Nc, double Phic, double vc, double c_sld, double bc, |
---|
| 5 | double Nd, double Phid, double vd, double d_sld, double bd, |
---|
| 6 | double Kab, double Kac, double Kad, |
---|
| 7 | double Kbc, double Kbd, double Kcd |
---|
| 8 | ); |
---|
| 9 | |
---|
| 10 | double Iqxy(double qx, double qy, double case_num, |
---|
| 11 | double Na, double Phia, double va, double a_sld, double ba, |
---|
| 12 | double Nb, double Phib, double vb, double b_sld, double bb, |
---|
| 13 | double Nc, double Phic, double vc, double c_sld, double bc, |
---|
| 14 | double Nd, double Phid, double vd, double d_sld, double bd, |
---|
| 15 | double Kab, double Kac, double Kad, |
---|
| 16 | double Kbc, double Kbd, double Kcd |
---|
| 17 | ); |
---|
| 18 | |
---|
| 19 | double form_volume(void); |
---|
| 20 | |
---|
| 21 | double form_volume(void) |
---|
| 22 | { |
---|
| 23 | return 1.0; |
---|
| 24 | } |
---|
| 25 | |
---|
| 26 | double Iq(double q, double case_num, |
---|
| 27 | double Na, double Phia, double va, double La, double ba, |
---|
| 28 | double Nb, double Phib, double vb, double Lb, double bb, |
---|
| 29 | double Nc, double Phic, double vc, double Lc, double bc, |
---|
| 30 | double Nd, double Phid, double vd, double Ld, double bd, |
---|
| 31 | double Kab, double Kac, double Kad, |
---|
| 32 | double Kbc, double Kbd, double Kcd |
---|
| 33 | ) { |
---|
| 34 | int icase = (int)case_num; |
---|
| 35 | |
---|
| 36 | #if 0 // Sasview defaults |
---|
| 37 | if (icase <= 1) { |
---|
| 38 | Na=Nb=1000.0; |
---|
| 39 | Phia=Phib=0.0000001; |
---|
| 40 | Kab=Kac=Kad=Kbc=Kbd=-0.0004; |
---|
| 41 | La=Lb=1e-12; |
---|
| 42 | va=vb=100.0; |
---|
| 43 | ba=bb=5.0; |
---|
| 44 | } else if (icase <= 4) { |
---|
| 45 | Phia=0.0000001; |
---|
| 46 | Kab=Kac=Kad=-0.0004; |
---|
| 47 | La=1e-12; |
---|
| 48 | va=100.0; |
---|
| 49 | ba=5.0; |
---|
| 50 | } |
---|
| 51 | #else |
---|
| 52 | if (icase <= 1) { |
---|
| 53 | Na=Nb=0.0; |
---|
| 54 | Phia=Phib=0.0; |
---|
| 55 | Kab=Kac=Kad=Kbc=Kbd=0.0; |
---|
| 56 | La=Lb=Ld; |
---|
| 57 | va=vb=vd; |
---|
| 58 | ba=bb=0.0; |
---|
| 59 | } else if (icase <= 4) { |
---|
| 60 | Na = 0.0; |
---|
| 61 | Phia=0.0; |
---|
| 62 | Kab=Kac=Kad=0.0; |
---|
| 63 | La=Ld; |
---|
| 64 | va=vd; |
---|
| 65 | ba=0.0; |
---|
| 66 | } |
---|
| 67 | #endif |
---|
| 68 | |
---|
| 69 | const double Xa = q*q*ba*ba*Na/6.0; |
---|
| 70 | const double Xb = q*q*bb*bb*Nb/6.0; |
---|
| 71 | const double Xc = q*q*bc*bc*Nc/6.0; |
---|
| 72 | const double Xd = q*q*bd*bd*Nd/6.0; |
---|
| 73 | |
---|
| 74 | // limit as Xa goes to 0 is 1 |
---|
| 75 | const double Pa = Xa==0 ? 1.0 : -expm1(-Xa)/Xa; |
---|
| 76 | const double Pb = Xb==0 ? 1.0 : -expm1(-Xb)/Xb; |
---|
| 77 | const double Pc = Xc==0 ? 1.0 : -expm1(-Xc)/Xc; |
---|
| 78 | const double Pd = Xd==0 ? 1.0 : -expm1(-Xd)/Xd; |
---|
| 79 | |
---|
| 80 | // limit as Xa goes to 0 is 1 |
---|
| 81 | const double Paa = Xa==0 ? 1.0 : 2.0*(1.0-Pa)/Xa; |
---|
| 82 | const double Pbb = Xb==0 ? 1.0 : 2.0*(1.0-Pb)/Xb; |
---|
| 83 | const double Pcc = Xc==0 ? 1.0 : 2.0*(1.0-Pc)/Xc; |
---|
| 84 | const double Pdd = Xd==0 ? 1.0 : 2.0*(1.0-Pd)/Xd; |
---|
| 85 | |
---|
| 86 | |
---|
| 87 | // Note: S0ij only defined for copolymers; otherwise set to zero |
---|
| 88 | // 0: C/D binary mixture |
---|
| 89 | // 1: C-D diblock copolymer |
---|
| 90 | // 2: B/C/D ternery mixture |
---|
| 91 | // 3: B/C-D binary mixture,1 homopolymer, 1 diblock copolymer |
---|
| 92 | // 4: B-C-D triblock copolymer |
---|
| 93 | // 5: A/B/C/D quaternary mixture |
---|
| 94 | // 6: A/B/C-D ternery mixture, 2 homopolymer, 1 diblock copolymer |
---|
| 95 | // 7: A/B-C-D binary mixture, 1 homopolymer, 1 triblock copolymer |
---|
| 96 | // 8: A-B/C-D binary mixture, 2 diblock copolymer |
---|
| 97 | // 9: A-B-C-D tetra-block copolymer |
---|
| 98 | #if 0 |
---|
| 99 | const double S0aa = icase<5 |
---|
| 100 | ? 1.0 : Na*Phia*va*Paa; |
---|
| 101 | const double S0bb = icase<2 |
---|
| 102 | ? 1.0 : Nb*Phib*vb*Pbb; |
---|
| 103 | const double S0cc = Nc*Phic*vc*Pcc; |
---|
| 104 | const double S0dd = Nd*Phid*vd*Pdd; |
---|
| 105 | const double S0ab = icase<8 |
---|
| 106 | ? 0.0 : sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb; |
---|
| 107 | const double S0ac = icase<9 |
---|
| 108 | ? 0.0 : sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb); |
---|
| 109 | const double S0ad = icase<9 |
---|
| 110 | ? 0.0 : sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc); |
---|
| 111 | const double S0bc = (icase!=4 && icase!=7 && icase!= 9) |
---|
| 112 | ? 0.0 : sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc; |
---|
| 113 | const double S0bd = (icase!=4 && icase!=7 && icase!= 9) |
---|
| 114 | ? 0.0 : sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc); |
---|
| 115 | const double S0cd = (icase==0 || icase==2 || icase==5) |
---|
| 116 | ? 0.0 : sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd; |
---|
| 117 | #else // sasview equivalent |
---|
| 118 | //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,Nc,Phic,vc,Pcc); |
---|
| 119 | double S0aa = Na*Phia*va*Paa; |
---|
| 120 | double S0bb = Nb*Phib*vb*Pbb; |
---|
| 121 | double S0cc = Nc*Phic*vc*Pcc; |
---|
| 122 | double S0dd = Nd*Phid*vd*Pdd; |
---|
| 123 | double S0ab = sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb; |
---|
| 124 | double S0ac = sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb); |
---|
| 125 | double S0ad = sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc); |
---|
| 126 | double S0bc = sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc; |
---|
| 127 | double S0bd = sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc); |
---|
| 128 | double S0cd = sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd; |
---|
| 129 | switch(icase){ |
---|
| 130 | case 0: |
---|
| 131 | S0aa=0.000001; |
---|
| 132 | S0ab=0.000002; |
---|
| 133 | S0ac=0.000003; |
---|
| 134 | S0ad=0.000004; |
---|
| 135 | S0bb=0.000005; |
---|
| 136 | S0bc=0.000006; |
---|
| 137 | S0bd=0.000007; |
---|
| 138 | S0cd=0.000008; |
---|
| 139 | break; |
---|
| 140 | case 1: |
---|
| 141 | S0aa=0.000001; |
---|
| 142 | S0ab=0.000002; |
---|
| 143 | S0ac=0.000003; |
---|
| 144 | S0ad=0.000004; |
---|
| 145 | S0bb=0.000005; |
---|
| 146 | S0bc=0.000006; |
---|
| 147 | S0bd=0.000007; |
---|
| 148 | break; |
---|
| 149 | case 2: |
---|
| 150 | S0aa=0.000001; |
---|
| 151 | S0ab=0.000002; |
---|
| 152 | S0ac=0.000003; |
---|
| 153 | S0ad=0.000004; |
---|
| 154 | S0bc=0.000005; |
---|
| 155 | S0bd=0.000006; |
---|
| 156 | S0cd=0.000007; |
---|
| 157 | break; |
---|
| 158 | case 3: |
---|
| 159 | S0aa=0.000001; |
---|
| 160 | S0ab=0.000002; |
---|
| 161 | S0ac=0.000003; |
---|
| 162 | S0ad=0.000004; |
---|
| 163 | S0bc=0.000005; |
---|
| 164 | S0bd=0.000006; |
---|
| 165 | break; |
---|
| 166 | case 4: |
---|
| 167 | S0aa=0.000001; |
---|
| 168 | S0ab=0.000002; |
---|
| 169 | S0ac=0.000003; |
---|
| 170 | S0ad=0.000004; |
---|
| 171 | break; |
---|
| 172 | case 5: |
---|
| 173 | S0ab=0.000001; |
---|
| 174 | S0ac=0.000002; |
---|
| 175 | S0ad=0.000003; |
---|
| 176 | S0bc=0.000004; |
---|
| 177 | S0bd=0.000005; |
---|
| 178 | S0cd=0.000006; |
---|
| 179 | break; |
---|
| 180 | case 6: |
---|
| 181 | S0ab=0.000001; |
---|
| 182 | S0ac=0.000002; |
---|
| 183 | S0ad=0.000003; |
---|
| 184 | S0bc=0.000004; |
---|
| 185 | S0bd=0.000005; |
---|
| 186 | break; |
---|
| 187 | case 7: |
---|
| 188 | S0ab=0.000001; |
---|
| 189 | S0ac=0.000002; |
---|
| 190 | S0ad=0.000003; |
---|
| 191 | break; |
---|
| 192 | case 8: |
---|
| 193 | S0ac=0.000001; |
---|
| 194 | S0ad=0.000002; |
---|
| 195 | S0bc=0.000003; |
---|
| 196 | S0bd=0.000004; |
---|
| 197 | break; |
---|
| 198 | default : //case 9: |
---|
| 199 | break; |
---|
| 200 | } |
---|
| 201 | #endif |
---|
| 202 | |
---|
| 203 | // eq 12a: \kappa_{ij}^F = \chi_{ij}^F - \chi_{i0}^F - \chi_{j0}^F |
---|
| 204 | const double Kaa = 0.0; |
---|
| 205 | const double Kbb = 0.0; |
---|
| 206 | const double Kcc = 0.0; |
---|
| 207 | const double Kdd = 0.0; |
---|
| 208 | const double Zaa = Kaa - Kad - Kad; |
---|
| 209 | const double Zab = Kab - Kad - Kbd; |
---|
| 210 | const double Zac = Kac - Kad - Kcd; |
---|
| 211 | const double Zbb = Kbb - Kbd - Kbd; |
---|
| 212 | const double Zbc = Kbc - Kbd - Kcd; |
---|
| 213 | const double Zcc = Kcc - Kcd - Kcd; |
---|
| 214 | //printf("Za: %10.5g %10.5g %10.5g\n", Zaa, Zab, Zac); |
---|
| 215 | //printf("Zb: %10.5g %10.5g %10.5g\n", Zab, Zbb, Zbc); |
---|
| 216 | //printf("Zc: %10.5g %10.5g %10.5g\n", Zac, Zbc, Zcc); |
---|
| 217 | |
---|
| 218 | // T = inv(S0) |
---|
| 219 | const double DenT = (- S0ac*S0bb*S0ac + S0ab*S0bc*S0ac + S0ac*S0ab*S0bc |
---|
| 220 | - S0aa*S0bc*S0bc - S0ab*S0ab*S0cc + S0aa*S0bb*S0cc); |
---|
| 221 | const double T11 = (-S0bc*S0bc + S0bb*S0cc)/DenT; |
---|
| 222 | const double T12 = ( S0ac*S0bc - S0ab*S0cc)/DenT; |
---|
| 223 | const double T13 = (-S0ac*S0bb + S0ab*S0bc)/DenT; |
---|
| 224 | const double T22 = (-S0ac*S0ac + S0aa*S0cc)/DenT; |
---|
| 225 | const double T23 = ( S0ac*S0ab - S0aa*S0bc)/DenT; |
---|
| 226 | const double T33 = (-S0ab*S0ab + S0aa*S0bb)/DenT; |
---|
| 227 | |
---|
| 228 | //printf("T1: %10.5g %10.5g %10.5g\n", T11, T12, T13); |
---|
| 229 | //printf("T2: %10.5g %10.5g %10.5g\n", T12, T22, T23); |
---|
| 230 | //printf("T3: %10.5g %10.5g %10.5g\n", T13, T23, T33); |
---|
| 231 | |
---|
| 232 | // eq 18e: m = 1/(S0_{dd} - s0^T inv(S0) s0) |
---|
| 233 | const double ZZ = S0ad*(T11*S0ad + T12*S0bd + T13*S0cd) |
---|
| 234 | + S0bd*(T12*S0ad + T22*S0bd + T23*S0cd) |
---|
| 235 | + S0cd*(T13*S0ad + T23*S0bd + T33*S0cd); |
---|
| 236 | |
---|
| 237 | const double m=1.0/(S0dd-ZZ); |
---|
| 238 | |
---|
| 239 | // eq 18d: Y = inv(S0)s0 + e |
---|
| 240 | const double Y1 = T11*S0ad + T12*S0bd + T13*S0cd + 1.0; |
---|
| 241 | const double Y2 = T12*S0ad + T22*S0bd + T23*S0cd + 1.0; |
---|
| 242 | const double Y3 = T13*S0ad + T23*S0bd + T33*S0cd + 1.0; |
---|
| 243 | |
---|
| 244 | // N = mYY^T + \kappa^F |
---|
| 245 | const double N11 = m*Y1*Y1 + Zaa; |
---|
| 246 | const double N12 = m*Y1*Y2 + Zab; |
---|
| 247 | const double N13 = m*Y1*Y3 + Zac; |
---|
| 248 | const double N22 = m*Y2*Y2 + Zbb; |
---|
| 249 | const double N23 = m*Y2*Y3 + Zbc; |
---|
| 250 | const double N33 = m*Y3*Y3 + Zcc; |
---|
| 251 | |
---|
| 252 | //printf("N1: %10.5g %10.5g %10.5g\n", N11, N12, N13); |
---|
| 253 | //printf("N2: %10.5g %10.5g %10.5g\n", N12, N22, N23); |
---|
| 254 | //printf("N3: %10.5g %10.5g %10.5g\n", N13, N23, N33); |
---|
| 255 | //printf("S0a: %10.5g %10.5g %10.5g\n", S0aa, S0ab, S0ac); |
---|
| 256 | //printf("S0b: %10.5g %10.5g %10.5g\n", S0ab, S0bb, S0bc); |
---|
| 257 | //printf("S0c: %10.5g %10.5g %10.5g\n", S0ac, S0bc, S0cc); |
---|
| 258 | |
---|
| 259 | // M = I + S0 N |
---|
| 260 | const double Maa = N11*S0aa + N12*S0ab + N13*S0ac + 1.0; |
---|
| 261 | const double Mab = N11*S0ab + N12*S0bb + N13*S0bc; |
---|
| 262 | const double Mac = N11*S0ac + N12*S0bc + N13*S0cc; |
---|
| 263 | const double Mba = N12*S0aa + N22*S0ab + N23*S0ac; |
---|
| 264 | const double Mbb = N12*S0ab + N22*S0bb + N23*S0bc + 1.0; |
---|
| 265 | const double Mbc = N12*S0ac + N22*S0bc + N23*S0cc; |
---|
| 266 | const double Mca = N13*S0aa + N23*S0ab + N33*S0ac; |
---|
| 267 | const double Mcb = N13*S0ab + N23*S0bb + N33*S0bc; |
---|
| 268 | const double Mcc = N13*S0ac + N23*S0bc + N33*S0cc + 1.0; |
---|
| 269 | //printf("M1: %10.5g %10.5g %10.5g\n", Maa, Mab, Mac); |
---|
| 270 | //printf("M2: %10.5g %10.5g %10.5g\n", Mba, Mbb, Mbc); |
---|
| 271 | //printf("M3: %10.5g %10.5g %10.5g\n", Mca, Mcb, Mcc); |
---|
| 272 | |
---|
| 273 | // Q = inv(M) = inv(I + S0 N) |
---|
| 274 | const double DenQ = (+ Maa*Mbb*Mcc - Maa*Mbc*Mcb - Mab*Mba*Mcc |
---|
| 275 | + Mab*Mbc*Mca + Mac*Mba*Mcb - Mac*Mbb*Mca); |
---|
| 276 | |
---|
| 277 | const double Q11 = ( Mbb*Mcc - Mbc*Mcb)/DenQ; |
---|
| 278 | const double Q12 = (-Mab*Mcc + Mac*Mcb)/DenQ; |
---|
| 279 | const double Q13 = ( Mab*Mbc - Mac*Mbb)/DenQ; |
---|
| 280 | const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ; |
---|
| 281 | const double Q22 = ( Maa*Mcc - Mac*Mca)/DenQ; |
---|
| 282 | const double Q23 = (-Maa*Mbc + Mac*Mba)/DenQ; |
---|
| 283 | const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ; |
---|
| 284 | const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ; |
---|
| 285 | const double Q33 = ( Maa*Mbb - Mab*Mba)/DenQ; |
---|
| 286 | |
---|
| 287 | //printf("Q1: %10.5g %10.5g %10.5g\n", Q11, Q12, Q13); |
---|
| 288 | //printf("Q2: %10.5g %10.5g %10.5g\n", Q21, Q22, Q23); |
---|
| 289 | //printf("Q3: %10.5g %10.5g %10.5g\n", Q31, Q32, Q33); |
---|
| 290 | // eq 18c: inv(S) = inv(S0) + mYY^T + \kappa^F |
---|
| 291 | // eq A1 in the appendix |
---|
| 292 | // To solve for S, use: |
---|
| 293 | // S = inv(inv(S^0) + N) inv(S^0) S^0 |
---|
| 294 | // = inv(S^0 inv(S^0) + N) S^0 |
---|
| 295 | // = inv(I + S^0 N) S^0 |
---|
| 296 | // = Q S^0 |
---|
| 297 | const double S11 = Q11*S0aa + Q12*S0ab + Q13*S0ac; |
---|
| 298 | const double S12 = Q12*S0aa + Q22*S0ab + Q23*S0ac; |
---|
| 299 | const double S13 = Q13*S0aa + Q23*S0ab + Q33*S0ac; |
---|
| 300 | const double S22 = Q12*S0ab + Q22*S0bb + Q23*S0bc; |
---|
| 301 | const double S23 = Q13*S0ab + Q23*S0bb + Q33*S0bc; |
---|
| 302 | const double S33 = Q13*S0ac + Q23*S0bc + Q33*S0cc; |
---|
| 303 | // If the full S is needed...it isn't since Ldd = (rho_d - rho_d) = 0 below |
---|
| 304 | //const double S14=-S11-S12-S13; |
---|
| 305 | //const double S24=-S12-S22-S23; |
---|
| 306 | //const double S34=-S13-S23-S33; |
---|
| 307 | //const double S44=S11+S22+S33 + 2.0*(S12+S13+S23); |
---|
| 308 | |
---|
| 309 | // eq 12 of Akcasu, 1990: I(q) = L^T S L |
---|
| 310 | // Note: eliminate cases without A and B polymers by setting Lij to 0 |
---|
| 311 | // Note: 1e-13 to convert from fm to cm for scattering length |
---|
| 312 | const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; |
---|
| 313 | const double Lad = icase<5 ? 0.0 : (La/va - Ld/vd)*sqrt_Nav; |
---|
| 314 | const double Lbd = icase<2 ? 0.0 : (Lb/vb - Ld/vd)*sqrt_Nav; |
---|
| 315 | const double Lcd = (Lc/vc - Ld/vd)*sqrt_Nav; |
---|
| 316 | |
---|
| 317 | const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 |
---|
| 318 | + 2.0*(Lad*Lbd*S12 + Lbd*Lcd*S23 + Lad*Lcd*S13); |
---|
| 319 | |
---|
| 320 | return result; |
---|
| 321 | |
---|
| 322 | } |
---|
| 323 | |
---|
| 324 | double Iqxy(double qx, double qy, |
---|
| 325 | double case_num, |
---|
| 326 | double Na, double Phia, double va, double a_sld, double ba, |
---|
| 327 | double Nb, double Phib, double vb, double b_sld, double bb, |
---|
| 328 | double Nc, double Phic, double vc, double c_sld, double bc, |
---|
| 329 | double Nd, double Phid, double vd, double d_sld, double bd, |
---|
| 330 | double Kab, double Kac, double Kad, |
---|
| 331 | double Kbc, double Kbd, double Kcd |
---|
| 332 | ) |
---|
| 333 | { |
---|
| 334 | double q = sqrt(qx*qx + qy*qy); |
---|
| 335 | return Iq(q, |
---|
| 336 | case_num, |
---|
| 337 | Na, Phia, va, a_sld, ba, |
---|
| 338 | Nb, Phib, vb, b_sld, bb, |
---|
| 339 | Nc, Phic, vc, c_sld, bc, |
---|
| 340 | Nd, Phid, vd, d_sld, bd, |
---|
| 341 | Kab, Kac, Kad, |
---|
| 342 | Kbc, Kbd, Kcd); |
---|
| 343 | } |
---|