[fa4a994] | 1 | double Iq(double q, double fp_case_num, |
---|
[d2bb604] | 2 | double N[], double Phi[], double v[], double L[], double b[], |
---|
[82c299f] | 3 | double Kab, double Kac, double Kad, |
---|
| 4 | double Kbc, double Kbd, double Kcd |
---|
| 5 | ); |
---|
| 6 | |
---|
[fa4a994] | 7 | double Iq(double q, double fp_case_num, |
---|
[51ec7e8] | 8 | double N[], // DEGREE OF POLYMERIZATION |
---|
| 9 | double Phi[], // VOL FRACTION |
---|
| 10 | double v[], // SPECIFIC VOLUME |
---|
| 11 | double L[], // SCATT. LENGTH |
---|
| 12 | double b[], // SEGMENT LENGTH |
---|
| 13 | double Kab, double Kac, double Kad, // CHI PARAM |
---|
[82c299f] | 14 | double Kbc, double Kbd, double Kcd |
---|
[51ec7e8] | 15 | ) |
---|
| 16 | { |
---|
[fa4a994] | 17 | int icase = (int)(fp_case_num+0.5); |
---|
[82c299f] | 18 | |
---|
[51ec7e8] | 19 | double Nab,Nac,Nad,Nbc,Nbd,Ncd; |
---|
| 20 | double Phiab,Phiac,Phiad,Phibc,Phibd,Phicd; |
---|
| 21 | double vab,vac,vad,vbc,vbd,vcd; |
---|
| 22 | double m; |
---|
| 23 | double Xa,Xb,Xc,Xd; |
---|
| 24 | double Paa,S0aa,Pab,S0ab,Pac,S0ac,Pad,S0ad; |
---|
| 25 | double S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd; |
---|
| 26 | double S0ca,S0cb,Pcc,S0cc,Pcd,S0cd; |
---|
[6351bfa] | 27 | double S0da,S0db,S0dc; |
---|
| 28 | double Pdd,S0dd; |
---|
| 29 | double Kaa,Kbb,Kcc; |
---|
| 30 | double Kba,Kca,Kcb; |
---|
| 31 | double Kda,Kdb,Kdc,Kdd; |
---|
[51ec7e8] | 32 | double Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc; |
---|
| 33 | double DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33; |
---|
| 34 | double Y1,Y2,Y3,X11,X12,X13,X21,X22,X23,X31,X32,X33; |
---|
| 35 | double ZZ,DenQ1,DenQ2,DenQ3,DenQ,Q11,Q12,Q13,Q21,Q22,Q23,Q31,Q32,Q33; |
---|
| 36 | double N11,N12,N13,N21,N22,N23,N31,N32,N33; |
---|
| 37 | double M11,M12,M13,M21,M22,M23,M31,M32,M33; |
---|
| 38 | double S11,S12,S13,S14,S21,S22,S23,S24; |
---|
| 39 | double S31,S32,S33,S34,S41,S42,S43,S44; |
---|
| 40 | double Lad,Lbd,Lcd,Nav,Intg; |
---|
[4329539] | 41 | |
---|
| 42 | // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc) |
---|
[51ec7e8] | 43 | //icase was shifted to N-1 from the original code |
---|
| 44 | if (icase <= 1){ |
---|
| 45 | Phi[0] = Phi[1] = 0.0000001; |
---|
| 46 | N[0] = N[1] = 1000.0; |
---|
| 47 | L[0] = L[1] = 1.e-12; |
---|
| 48 | v[0] = v[1] = 100.0; |
---|
| 49 | b[0] = b[1] = 5.0; |
---|
| 50 | Kab = Kac = Kad = Kbc = Kbd = -0.0004; |
---|
| 51 | } |
---|
| 52 | else if ((icase > 1) && (icase <= 4)){ |
---|
| 53 | Phi[0] = 0.0000001; |
---|
| 54 | N[0] = 1000.0; |
---|
| 55 | L[0] = 1.e-12; |
---|
| 56 | v[0] = 100.0; |
---|
| 57 | b[0] = 5.0; |
---|
| 58 | Kab = Kac = Kad = -0.0004; |
---|
| 59 | } |
---|
[4329539] | 60 | |
---|
[bb73096] | 61 | // Set volume fraction of component D based on constraint that sum of vol frac =1 |
---|
| 62 | Phi[3]=1.0-Phi[0]-Phi[1]-Phi[2]; |
---|
[51ec7e8] | 63 | |
---|
[4329539] | 64 | //set up values for cross terms in case of block copolymers (1,3,4,6,7,8,9) |
---|
[51ec7e8] | 65 | Nab=sqrt(N[0]*N[1]); |
---|
| 66 | Nac=sqrt(N[0]*N[2]); |
---|
| 67 | Nad=sqrt(N[0]*N[3]); |
---|
| 68 | Nbc=sqrt(N[1]*N[2]); |
---|
| 69 | Nbd=sqrt(N[1]*N[3]); |
---|
| 70 | Ncd=sqrt(N[2]*N[3]); |
---|
| 71 | |
---|
| 72 | vab=sqrt(v[0]*v[1]); |
---|
| 73 | vac=sqrt(v[0]*v[2]); |
---|
| 74 | vad=sqrt(v[0]*v[3]); |
---|
| 75 | vbc=sqrt(v[1]*v[2]); |
---|
| 76 | vbd=sqrt(v[1]*v[3]); |
---|
| 77 | vcd=sqrt(v[2]*v[3]); |
---|
| 78 | |
---|
| 79 | Phiab=sqrt(Phi[0]*Phi[1]); |
---|
| 80 | Phiac=sqrt(Phi[0]*Phi[2]); |
---|
| 81 | Phiad=sqrt(Phi[0]*Phi[3]); |
---|
| 82 | Phibc=sqrt(Phi[1]*Phi[2]); |
---|
| 83 | Phibd=sqrt(Phi[1]*Phi[3]); |
---|
| 84 | Phicd=sqrt(Phi[2]*Phi[3]); |
---|
| 85 | |
---|
[4329539] | 86 | // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk |
---|
[51ec7e8] | 87 | Xa=q*q*b[0]*b[0]*N[0]/6.0; |
---|
| 88 | Xb=q*q*b[1]*b[1]*N[1]/6.0; |
---|
[d680a2b] | 89 | Xc=q*q*b[2]*b[2]*N[2]/6.0; |
---|
| 90 | Xd=q*q*b[3]*b[3]*N[3]/6.0; |
---|
[51ec7e8] | 91 | |
---|
[4329539] | 92 | //calculate all partial structure factors Pij and normalize n^2 |
---|
| 93 | Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); // free A chain form factor |
---|
| 94 | S0aa=N[0]*Phi[0]*v[0]*Paa; // Phi * Vp * P(Q)= I(Q0)/delRho^2 |
---|
| 95 | Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); //AB diblock (anchored Paa * anchored Pbb) partial form factor |
---|
[51ec7e8] | 96 | S0ab=(Phiab*vab*Nab)*Pab; |
---|
[4329539] | 97 | Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor |
---|
[51ec7e8] | 98 | S0ac=(Phiac*vac*Nac)*Pac; |
---|
[4329539] | 99 | Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block |
---|
[51ec7e8] | 100 | S0ad=(Phiad*vad*Nad)*Pad; |
---|
| 101 | |
---|
| 102 | S0ba=S0ab; |
---|
[4329539] | 103 | Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain |
---|
[51ec7e8] | 104 | S0bb=N[1]*Phi[1]*v[1]*Pbb; |
---|
[4329539] | 105 | Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock |
---|
[51ec7e8] | 106 | S0bc=(Phibc*vbc*Nbc)*Pbc; |
---|
[4329539] | 107 | Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock |
---|
[51ec7e8] | 108 | S0bd=(Phibd*vbd*Nbd)*Pbd; |
---|
| 109 | |
---|
| 110 | S0ca=S0ac; |
---|
| 111 | S0cb=S0bc; |
---|
[4329539] | 112 | Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain |
---|
[51ec7e8] | 113 | S0cc=N[2]*Phi[2]*v[2]*Pcc; |
---|
[4329539] | 114 | Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock |
---|
[51ec7e8] | 115 | S0cd=(Phicd*vcd*Ncd)*Pcd; |
---|
| 116 | |
---|
| 117 | S0da=S0ad; |
---|
| 118 | S0db=S0bd; |
---|
| 119 | S0dc=S0cd; |
---|
[4329539] | 120 | Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain |
---|
[51ec7e8] | 121 | S0dd=N[3]*Phi[3]*v[3]*Pdd; |
---|
[4329539] | 122 | |
---|
| 123 | // Reset all unused partial structure factors to 0 (depends on case) |
---|
[51ec7e8] | 124 | //icase was shifted to N-1 from the original code |
---|
| 125 | switch(icase){ |
---|
| 126 | case 0: |
---|
| 127 | S0aa=0.000001; |
---|
| 128 | S0ab=0.000002; |
---|
| 129 | S0ac=0.000003; |
---|
| 130 | S0ad=0.000004; |
---|
| 131 | S0bb=0.000005; |
---|
| 132 | S0bc=0.000006; |
---|
| 133 | S0bd=0.000007; |
---|
| 134 | S0cd=0.000008; |
---|
| 135 | break; |
---|
| 136 | case 1: |
---|
| 137 | S0aa=0.000001; |
---|
| 138 | S0ab=0.000002; |
---|
| 139 | S0ac=0.000003; |
---|
| 140 | S0ad=0.000004; |
---|
| 141 | S0bb=0.000005; |
---|
| 142 | S0bc=0.000006; |
---|
| 143 | S0bd=0.000007; |
---|
| 144 | break; |
---|
| 145 | case 2: |
---|
| 146 | S0aa=0.000001; |
---|
| 147 | S0ab=0.000002; |
---|
| 148 | S0ac=0.000003; |
---|
| 149 | S0ad=0.000004; |
---|
| 150 | S0bc=0.000005; |
---|
| 151 | S0bd=0.000006; |
---|
| 152 | S0cd=0.000007; |
---|
| 153 | break; |
---|
| 154 | case 3: |
---|
| 155 | S0aa=0.000001; |
---|
| 156 | S0ab=0.000002; |
---|
| 157 | S0ac=0.000003; |
---|
| 158 | S0ad=0.000004; |
---|
| 159 | S0bc=0.000005; |
---|
| 160 | S0bd=0.000006; |
---|
| 161 | break; |
---|
| 162 | case 4: |
---|
| 163 | S0aa=0.000001; |
---|
| 164 | S0ab=0.000002; |
---|
| 165 | S0ac=0.000003; |
---|
| 166 | S0ad=0.000004; |
---|
| 167 | break; |
---|
| 168 | case 5: |
---|
| 169 | S0ab=0.000001; |
---|
| 170 | S0ac=0.000002; |
---|
| 171 | S0ad=0.000003; |
---|
| 172 | S0bc=0.000004; |
---|
| 173 | S0bd=0.000005; |
---|
| 174 | S0cd=0.000006; |
---|
| 175 | break; |
---|
| 176 | case 6: |
---|
| 177 | S0ab=0.000001; |
---|
| 178 | S0ac=0.000002; |
---|
| 179 | S0ad=0.000003; |
---|
| 180 | S0bc=0.000004; |
---|
| 181 | S0bd=0.000005; |
---|
| 182 | break; |
---|
| 183 | case 7: |
---|
| 184 | S0ab=0.000001; |
---|
| 185 | S0ac=0.000002; |
---|
| 186 | S0ad=0.000003; |
---|
| 187 | break; |
---|
| 188 | case 8: |
---|
| 189 | S0ac=0.000001; |
---|
| 190 | S0ad=0.000002; |
---|
| 191 | S0bc=0.000003; |
---|
| 192 | S0bd=0.000004; |
---|
| 193 | break; |
---|
| 194 | default : //case 9: |
---|
| 195 | break; |
---|
| 196 | } |
---|
| 197 | S0ba=S0ab; |
---|
| 198 | S0ca=S0ac; |
---|
| 199 | S0cb=S0bc; |
---|
| 200 | S0da=S0ad; |
---|
| 201 | S0db=S0bd; |
---|
| 202 | S0dc=S0cd; |
---|
| 203 | |
---|
[4329539] | 204 | // self chi parameter is 0 ... of course |
---|
[51ec7e8] | 205 | Kaa=0.0; |
---|
| 206 | Kbb=0.0; |
---|
| 207 | Kcc=0.0; |
---|
| 208 | Kdd=0.0; |
---|
| 209 | |
---|
| 210 | Kba=Kab; |
---|
| 211 | Kca=Kac; |
---|
| 212 | Kcb=Kbc; |
---|
| 213 | Kda=Kad; |
---|
| 214 | Kdb=Kbd; |
---|
| 215 | Kdc=Kcd; |
---|
| 216 | |
---|
| 217 | Zaa=Kaa-Kad-Kad; |
---|
| 218 | Zab=Kab-Kad-Kbd; |
---|
| 219 | Zac=Kac-Kad-Kcd; |
---|
| 220 | Zba=Kba-Kbd-Kad; |
---|
| 221 | Zbb=Kbb-Kbd-Kbd; |
---|
| 222 | Zbc=Kbc-Kbd-Kcd; |
---|
| 223 | Zca=Kca-Kcd-Kad; |
---|
| 224 | Zcb=Kcb-Kcd-Kbd; |
---|
| 225 | Zcc=Kcc-Kcd-Kcd; |
---|
| 226 | |
---|
| 227 | DenT=(-(S0ac*S0bb*S0ca) + S0ab*S0bc*S0ca + S0ac*S0ba*S0cb - S0aa*S0bc*S0cb - S0ab*S0ba*S0cc + S0aa*S0bb*S0cc); |
---|
| 228 | |
---|
| 229 | T11= (-(S0bc*S0cb) + S0bb*S0cc)/DenT; |
---|
| 230 | T12= (S0ac*S0cb - S0ab*S0cc)/DenT; |
---|
| 231 | T13= (-(S0ac*S0bb) + S0ab*S0bc)/DenT; |
---|
| 232 | T21= (S0bc*S0ca - S0ba*S0cc)/DenT; |
---|
| 233 | T22= (-(S0ac*S0ca) + S0aa*S0cc)/DenT; |
---|
| 234 | T23= (S0ac*S0ba - S0aa*S0bc)/DenT; |
---|
| 235 | T31= (-(S0bb*S0ca) + S0ba*S0cb)/DenT; |
---|
| 236 | T32= (S0ab*S0ca - S0aa*S0cb)/DenT; |
---|
| 237 | T33= (-(S0ab*S0ba) + S0aa*S0bb)/DenT; |
---|
| 238 | |
---|
| 239 | Y1=T11*S0ad+T12*S0bd+T13*S0cd+1.0; |
---|
| 240 | Y2=T21*S0ad+T22*S0bd+T23*S0cd+1.0; |
---|
| 241 | Y3=T31*S0ad+T32*S0bd+T33*S0cd+1.0; |
---|
| 242 | |
---|
| 243 | X11=Y1*Y1; |
---|
| 244 | X12=Y1*Y2; |
---|
| 245 | X13=Y1*Y3; |
---|
| 246 | X21=Y2*Y1; |
---|
| 247 | X22=Y2*Y2; |
---|
| 248 | X23=Y2*Y3; |
---|
| 249 | X31=Y3*Y1; |
---|
| 250 | X32=Y3*Y2; |
---|
| 251 | X33=Y3*Y3; |
---|
| 252 | |
---|
| 253 | ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); |
---|
| 254 | |
---|
[4329539] | 255 | // D is considered the matrix or background component so enters here |
---|
[51ec7e8] | 256 | m=1.0/(S0dd-ZZ); |
---|
| 257 | |
---|
| 258 | N11=m*X11+Zaa; |
---|
| 259 | N12=m*X12+Zab; |
---|
| 260 | N13=m*X13+Zac; |
---|
| 261 | N21=m*X21+Zba; |
---|
| 262 | N22=m*X22+Zbb; |
---|
| 263 | N23=m*X23+Zbc; |
---|
| 264 | N31=m*X31+Zca; |
---|
| 265 | N32=m*X32+Zcb; |
---|
| 266 | N33=m*X33+Zcc; |
---|
| 267 | |
---|
| 268 | M11= N11*S0aa + N12*S0ab + N13*S0ac; |
---|
| 269 | M12= N11*S0ab + N12*S0bb + N13*S0bc; |
---|
| 270 | M13= N11*S0ac + N12*S0bc + N13*S0cc; |
---|
| 271 | M21= N21*S0aa + N22*S0ab + N23*S0ac; |
---|
| 272 | M22= N21*S0ab + N22*S0bb + N23*S0bc; |
---|
| 273 | M23= N21*S0ac + N22*S0bc + N23*S0cc; |
---|
| 274 | M31= N31*S0aa + N32*S0ab + N33*S0ac; |
---|
| 275 | M32= N31*S0ab + N32*S0bb + N33*S0bc; |
---|
| 276 | M33= N31*S0ac + N32*S0bc + N33*S0cc; |
---|
| 277 | |
---|
| 278 | DenQ1=1.0+M11-M12*M21+M22+M11*M22-M13*M31-M13*M22*M31; |
---|
| 279 | DenQ2= M12*M23*M31+M13*M21*M32-M23*M32-M11*M23*M32+M33+M11*M33; |
---|
| 280 | DenQ3= -M12*M21*M33+M22*M33+M11*M22*M33; |
---|
| 281 | DenQ=DenQ1+DenQ2+DenQ3; |
---|
| 282 | |
---|
| 283 | Q11= (1.0 + M22-M23*M32 + M33 + M22*M33)/DenQ; |
---|
| 284 | Q12= (-M12 + M13*M32 - M12*M33)/DenQ; |
---|
| 285 | Q13= (-M13 - M13*M22 + M12*M23)/DenQ; |
---|
| 286 | Q21= (-M21 + M23*M31 - M21*M33)/DenQ; |
---|
| 287 | Q22= (1.0 + M11 - M13*M31 + M33 + M11*M33)/DenQ; |
---|
| 288 | Q23= (M13*M21 - M23 - M11*M23)/DenQ; |
---|
| 289 | Q31= (-M31 - M22*M31 + M21*M32)/DenQ; |
---|
| 290 | Q32= (M12*M31 - M32 - M11*M32)/DenQ; |
---|
| 291 | Q33= (1.0 + M11 - M12*M21 + M22 + M11*M22)/DenQ; |
---|
| 292 | |
---|
| 293 | S11= Q11*S0aa + Q21*S0ab + Q31*S0ac; |
---|
| 294 | S12= Q12*S0aa + Q22*S0ab + Q32*S0ac; |
---|
| 295 | S13= Q13*S0aa + Q23*S0ab + Q33*S0ac; |
---|
| 296 | S14=-S11-S12-S13; |
---|
| 297 | S21= Q11*S0ba + Q21*S0bb + Q31*S0bc; |
---|
| 298 | S22= Q12*S0ba + Q22*S0bb + Q32*S0bc; |
---|
| 299 | S23= Q13*S0ba + Q23*S0bb + Q33*S0bc; |
---|
| 300 | S24=-S21-S22-S23; |
---|
| 301 | S31= Q11*S0ca + Q21*S0cb + Q31*S0cc; |
---|
| 302 | S32= Q12*S0ca + Q22*S0cb + Q32*S0cc; |
---|
| 303 | S33= Q13*S0ca + Q23*S0cb + Q33*S0cc; |
---|
| 304 | S34=-S31-S32-S33; |
---|
| 305 | S41=S14; |
---|
| 306 | S42=S24; |
---|
| 307 | S43=S34; |
---|
| 308 | S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; |
---|
| 309 | |
---|
[4329539] | 310 | //calculate contrast where L[i] is the scattering length of i and D is the matrix |
---|
[20c856a] | 311 | //Note that should multiply by Nav to get units of SLD which will become |
---|
| 312 | // Nav*2 in the next line (SLD^2) but then normalization in that line would |
---|
| 313 | //need to divide by Nav leaving only Nav or sqrt(Nav) before squaring. |
---|
[51ec7e8] | 314 | Nav=6.022045e+23; |
---|
| 315 | Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); |
---|
| 316 | Lbd=(L[1]/v[1]-L[3]/v[3])*sqrt(Nav); |
---|
| 317 | Lcd=(L[2]/v[2]-L[3]/v[3])*sqrt(Nav); |
---|
| 318 | |
---|
| 319 | Intg=Lad*Lad*S11+Lbd*Lbd*S22+Lcd*Lcd*S33+2.0*Lad*Lbd*S12+2.0*Lbd*Lcd*S23+2.0*Lad*Lcd*S13; |
---|
| 320 | |
---|
[20c856a] | 321 | //rescale for units of Lij^2 (fm^2 to cm^2) |
---|
| 322 | Intg *= 1.0e-26; |
---|
[acfb094] | 323 | |
---|
[51ec7e8] | 324 | return Intg; |
---|
| 325 | |
---|
| 326 | |
---|
| 327 | /* Attempts at a new implementation --- supressed for now |
---|
| 328 | #if 1 // Sasview defaults |
---|
[82c299f] | 329 | if (icase <= 1) { |
---|
[d2bb604] | 330 | N[0]=N[1]=1000.0; |
---|
| 331 | Phi[0]=Phi[1]=0.0000001; |
---|
[82c299f] | 332 | Kab=Kac=Kad=Kbc=Kbd=-0.0004; |
---|
[51ec7e8] | 333 | L[0]=L[1]=1.0e-12; |
---|
| 334 | v[0]=v[1]=100.0; |
---|
| 335 | b[0]=b[1]=5.0; |
---|
[82c299f] | 336 | } else if (icase <= 4) { |
---|
[d2bb604] | 337 | Phi[0]=0.0000001; |
---|
[82c299f] | 338 | Kab=Kac=Kad=-0.0004; |
---|
[51ec7e8] | 339 | L[0]=1.0e-12; |
---|
| 340 | v[0]=100.0; |
---|
| 341 | b[0]=5.0; |
---|
[82c299f] | 342 | } |
---|
| 343 | #else |
---|
| 344 | if (icase <= 1) { |
---|
[d2bb604] | 345 | N[0]=N[1]=0.0; |
---|
| 346 | Phi[0]=Phi[1]=0.0; |
---|
[82c299f] | 347 | Kab=Kac=Kad=Kbc=Kbd=0.0; |
---|
[d2bb604] | 348 | L[0]=L[1]=L[3]; |
---|
| 349 | v[0]=v[1]=v[3]; |
---|
| 350 | b[0]=b[1]=0.0; |
---|
[82c299f] | 351 | } else if (icase <= 4) { |
---|
[d2bb604] | 352 | N[0] = 0.0; |
---|
| 353 | Phi[0]=0.0; |
---|
[82c299f] | 354 | Kab=Kac=Kad=0.0; |
---|
[d2bb604] | 355 | L[0]=L[3]; |
---|
| 356 | v[0]=v[3]; |
---|
| 357 | b[0]=0.0; |
---|
[82c299f] | 358 | } |
---|
| 359 | #endif |
---|
| 360 | |
---|
[d2bb604] | 361 | const double Xa = q*q*b[0]*b[0]*N[0]/6.0; |
---|
| 362 | const double Xb = q*q*b[1]*b[1]*N[1]/6.0; |
---|
| 363 | const double Xc = q*q*b[2]*b[2]*N[2]/6.0; |
---|
| 364 | const double Xd = q*q*b[3]*b[3]*N[3]/6.0; |
---|
[82c299f] | 365 | |
---|
| 366 | // limit as Xa goes to 0 is 1 |
---|
| 367 | const double Pa = Xa==0 ? 1.0 : -expm1(-Xa)/Xa; |
---|
| 368 | const double Pb = Xb==0 ? 1.0 : -expm1(-Xb)/Xb; |
---|
| 369 | const double Pc = Xc==0 ? 1.0 : -expm1(-Xc)/Xc; |
---|
| 370 | const double Pd = Xd==0 ? 1.0 : -expm1(-Xd)/Xd; |
---|
| 371 | |
---|
| 372 | // limit as Xa goes to 0 is 1 |
---|
| 373 | const double Paa = Xa==0 ? 1.0 : 2.0*(1.0-Pa)/Xa; |
---|
| 374 | const double Pbb = Xb==0 ? 1.0 : 2.0*(1.0-Pb)/Xb; |
---|
| 375 | const double Pcc = Xc==0 ? 1.0 : 2.0*(1.0-Pc)/Xc; |
---|
| 376 | const double Pdd = Xd==0 ? 1.0 : 2.0*(1.0-Pd)/Xd; |
---|
| 377 | |
---|
| 378 | |
---|
| 379 | // Note: S0ij only defined for copolymers; otherwise set to zero |
---|
| 380 | // 0: C/D binary mixture |
---|
| 381 | // 1: C-D diblock copolymer |
---|
| 382 | // 2: B/C/D ternery mixture |
---|
| 383 | // 3: B/C-D binary mixture,1 homopolymer, 1 diblock copolymer |
---|
| 384 | // 4: B-C-D triblock copolymer |
---|
| 385 | // 5: A/B/C/D quaternary mixture |
---|
| 386 | // 6: A/B/C-D ternery mixture, 2 homopolymer, 1 diblock copolymer |
---|
| 387 | // 7: A/B-C-D binary mixture, 1 homopolymer, 1 triblock copolymer |
---|
| 388 | // 8: A-B/C-D binary mixture, 2 diblock copolymer |
---|
| 389 | // 9: A-B-C-D tetra-block copolymer |
---|
| 390 | #if 0 |
---|
| 391 | const double S0aa = icase<5 |
---|
[d2bb604] | 392 | ? 1.0 : N[0]*Phi[0]*v[0]*Paa; |
---|
[82c299f] | 393 | const double S0bb = icase<2 |
---|
[d2bb604] | 394 | ? 1.0 : N[1]*Phi[1]*v[1]*Pbb; |
---|
| 395 | const double S0cc = N[2]*Phi[2]*v[2]*Pcc; |
---|
| 396 | const double S0dd = N[3]*Phi[3]*v[3]*Pdd; |
---|
[82c299f] | 397 | const double S0ab = icase<8 |
---|
[d2bb604] | 398 | ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; |
---|
[82c299f] | 399 | const double S0ac = icase<9 |
---|
[d2bb604] | 400 | ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); |
---|
[82c299f] | 401 | const double S0ad = icase<9 |
---|
[d2bb604] | 402 | ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); |
---|
[82c299f] | 403 | const double S0bc = (icase!=4 && icase!=7 && icase!= 9) |
---|
[d2bb604] | 404 | ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; |
---|
[82c299f] | 405 | const double S0bd = (icase!=4 && icase!=7 && icase!= 9) |
---|
[d2bb604] | 406 | ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); |
---|
[82c299f] | 407 | const double S0cd = (icase==0 || icase==2 || icase==5) |
---|
[d2bb604] | 408 | ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; |
---|
[82c299f] | 409 | #else // sasview equivalent |
---|
[d2bb604] | 410 | //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc); |
---|
| 411 | double S0aa = N[0]*Phi[0]*v[0]*Paa; |
---|
| 412 | double S0bb = N[1]*Phi[1]*v[1]*Pbb; |
---|
| 413 | double S0cc = N[2]*Phi[2]*v[2]*Pcc; |
---|
| 414 | double S0dd = N[3]*Phi[3]*v[3]*Pdd; |
---|
| 415 | double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; |
---|
| 416 | double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); |
---|
| 417 | double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); |
---|
| 418 | double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; |
---|
| 419 | double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); |
---|
| 420 | double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; |
---|
[82c299f] | 421 | switch(icase){ |
---|
| 422 | case 0: |
---|
| 423 | S0aa=0.000001; |
---|
| 424 | S0ab=0.000002; |
---|
| 425 | S0ac=0.000003; |
---|
| 426 | S0ad=0.000004; |
---|
| 427 | S0bb=0.000005; |
---|
| 428 | S0bc=0.000006; |
---|
| 429 | S0bd=0.000007; |
---|
| 430 | S0cd=0.000008; |
---|
| 431 | break; |
---|
| 432 | case 1: |
---|
| 433 | S0aa=0.000001; |
---|
| 434 | S0ab=0.000002; |
---|
| 435 | S0ac=0.000003; |
---|
| 436 | S0ad=0.000004; |
---|
| 437 | S0bb=0.000005; |
---|
| 438 | S0bc=0.000006; |
---|
| 439 | S0bd=0.000007; |
---|
| 440 | break; |
---|
| 441 | case 2: |
---|
| 442 | S0aa=0.000001; |
---|
| 443 | S0ab=0.000002; |
---|
| 444 | S0ac=0.000003; |
---|
| 445 | S0ad=0.000004; |
---|
| 446 | S0bc=0.000005; |
---|
| 447 | S0bd=0.000006; |
---|
| 448 | S0cd=0.000007; |
---|
| 449 | break; |
---|
| 450 | case 3: |
---|
| 451 | S0aa=0.000001; |
---|
| 452 | S0ab=0.000002; |
---|
| 453 | S0ac=0.000003; |
---|
| 454 | S0ad=0.000004; |
---|
| 455 | S0bc=0.000005; |
---|
| 456 | S0bd=0.000006; |
---|
| 457 | break; |
---|
| 458 | case 4: |
---|
| 459 | S0aa=0.000001; |
---|
| 460 | S0ab=0.000002; |
---|
| 461 | S0ac=0.000003; |
---|
| 462 | S0ad=0.000004; |
---|
| 463 | break; |
---|
| 464 | case 5: |
---|
| 465 | S0ab=0.000001; |
---|
| 466 | S0ac=0.000002; |
---|
| 467 | S0ad=0.000003; |
---|
| 468 | S0bc=0.000004; |
---|
| 469 | S0bd=0.000005; |
---|
| 470 | S0cd=0.000006; |
---|
| 471 | break; |
---|
| 472 | case 6: |
---|
| 473 | S0ab=0.000001; |
---|
| 474 | S0ac=0.000002; |
---|
| 475 | S0ad=0.000003; |
---|
| 476 | S0bc=0.000004; |
---|
| 477 | S0bd=0.000005; |
---|
| 478 | break; |
---|
| 479 | case 7: |
---|
| 480 | S0ab=0.000001; |
---|
| 481 | S0ac=0.000002; |
---|
| 482 | S0ad=0.000003; |
---|
| 483 | break; |
---|
| 484 | case 8: |
---|
| 485 | S0ac=0.000001; |
---|
| 486 | S0ad=0.000002; |
---|
| 487 | S0bc=0.000003; |
---|
| 488 | S0bd=0.000004; |
---|
| 489 | break; |
---|
| 490 | default : //case 9: |
---|
| 491 | break; |
---|
| 492 | } |
---|
| 493 | #endif |
---|
| 494 | |
---|
| 495 | // eq 12a: \kappa_{ij}^F = \chi_{ij}^F - \chi_{i0}^F - \chi_{j0}^F |
---|
| 496 | const double Kaa = 0.0; |
---|
| 497 | const double Kbb = 0.0; |
---|
| 498 | const double Kcc = 0.0; |
---|
[13ed84c] | 499 | //const double Kdd = 0.0; |
---|
[82c299f] | 500 | const double Zaa = Kaa - Kad - Kad; |
---|
| 501 | const double Zab = Kab - Kad - Kbd; |
---|
| 502 | const double Zac = Kac - Kad - Kcd; |
---|
| 503 | const double Zbb = Kbb - Kbd - Kbd; |
---|
| 504 | const double Zbc = Kbc - Kbd - Kcd; |
---|
| 505 | const double Zcc = Kcc - Kcd - Kcd; |
---|
| 506 | //printf("Za: %10.5g %10.5g %10.5g\n", Zaa, Zab, Zac); |
---|
| 507 | //printf("Zb: %10.5g %10.5g %10.5g\n", Zab, Zbb, Zbc); |
---|
| 508 | //printf("Zc: %10.5g %10.5g %10.5g\n", Zac, Zbc, Zcc); |
---|
| 509 | |
---|
| 510 | // T = inv(S0) |
---|
| 511 | const double DenT = (- S0ac*S0bb*S0ac + S0ab*S0bc*S0ac + S0ac*S0ab*S0bc |
---|
| 512 | - S0aa*S0bc*S0bc - S0ab*S0ab*S0cc + S0aa*S0bb*S0cc); |
---|
| 513 | const double T11 = (-S0bc*S0bc + S0bb*S0cc)/DenT; |
---|
| 514 | const double T12 = ( S0ac*S0bc - S0ab*S0cc)/DenT; |
---|
| 515 | const double T13 = (-S0ac*S0bb + S0ab*S0bc)/DenT; |
---|
| 516 | const double T22 = (-S0ac*S0ac + S0aa*S0cc)/DenT; |
---|
| 517 | const double T23 = ( S0ac*S0ab - S0aa*S0bc)/DenT; |
---|
| 518 | const double T33 = (-S0ab*S0ab + S0aa*S0bb)/DenT; |
---|
| 519 | |
---|
| 520 | //printf("T1: %10.5g %10.5g %10.5g\n", T11, T12, T13); |
---|
| 521 | //printf("T2: %10.5g %10.5g %10.5g\n", T12, T22, T23); |
---|
| 522 | //printf("T3: %10.5g %10.5g %10.5g\n", T13, T23, T33); |
---|
| 523 | |
---|
| 524 | // eq 18e: m = 1/(S0_{dd} - s0^T inv(S0) s0) |
---|
| 525 | const double ZZ = S0ad*(T11*S0ad + T12*S0bd + T13*S0cd) |
---|
| 526 | + S0bd*(T12*S0ad + T22*S0bd + T23*S0cd) |
---|
| 527 | + S0cd*(T13*S0ad + T23*S0bd + T33*S0cd); |
---|
| 528 | |
---|
| 529 | const double m=1.0/(S0dd-ZZ); |
---|
| 530 | |
---|
| 531 | // eq 18d: Y = inv(S0)s0 + e |
---|
| 532 | const double Y1 = T11*S0ad + T12*S0bd + T13*S0cd + 1.0; |
---|
| 533 | const double Y2 = T12*S0ad + T22*S0bd + T23*S0cd + 1.0; |
---|
| 534 | const double Y3 = T13*S0ad + T23*S0bd + T33*S0cd + 1.0; |
---|
| 535 | |
---|
| 536 | // N = mYY^T + \kappa^F |
---|
| 537 | const double N11 = m*Y1*Y1 + Zaa; |
---|
| 538 | const double N12 = m*Y1*Y2 + Zab; |
---|
| 539 | const double N13 = m*Y1*Y3 + Zac; |
---|
| 540 | const double N22 = m*Y2*Y2 + Zbb; |
---|
| 541 | const double N23 = m*Y2*Y3 + Zbc; |
---|
| 542 | const double N33 = m*Y3*Y3 + Zcc; |
---|
| 543 | |
---|
| 544 | //printf("N1: %10.5g %10.5g %10.5g\n", N11, N12, N13); |
---|
| 545 | //printf("N2: %10.5g %10.5g %10.5g\n", N12, N22, N23); |
---|
| 546 | //printf("N3: %10.5g %10.5g %10.5g\n", N13, N23, N33); |
---|
| 547 | //printf("S0a: %10.5g %10.5g %10.5g\n", S0aa, S0ab, S0ac); |
---|
| 548 | //printf("S0b: %10.5g %10.5g %10.5g\n", S0ab, S0bb, S0bc); |
---|
| 549 | //printf("S0c: %10.5g %10.5g %10.5g\n", S0ac, S0bc, S0cc); |
---|
| 550 | |
---|
| 551 | // M = I + S0 N |
---|
| 552 | const double Maa = N11*S0aa + N12*S0ab + N13*S0ac + 1.0; |
---|
| 553 | const double Mab = N11*S0ab + N12*S0bb + N13*S0bc; |
---|
| 554 | const double Mac = N11*S0ac + N12*S0bc + N13*S0cc; |
---|
| 555 | const double Mba = N12*S0aa + N22*S0ab + N23*S0ac; |
---|
| 556 | const double Mbb = N12*S0ab + N22*S0bb + N23*S0bc + 1.0; |
---|
| 557 | const double Mbc = N12*S0ac + N22*S0bc + N23*S0cc; |
---|
| 558 | const double Mca = N13*S0aa + N23*S0ab + N33*S0ac; |
---|
| 559 | const double Mcb = N13*S0ab + N23*S0bb + N33*S0bc; |
---|
| 560 | const double Mcc = N13*S0ac + N23*S0bc + N33*S0cc + 1.0; |
---|
| 561 | //printf("M1: %10.5g %10.5g %10.5g\n", Maa, Mab, Mac); |
---|
| 562 | //printf("M2: %10.5g %10.5g %10.5g\n", Mba, Mbb, Mbc); |
---|
| 563 | //printf("M3: %10.5g %10.5g %10.5g\n", Mca, Mcb, Mcc); |
---|
| 564 | |
---|
| 565 | // Q = inv(M) = inv(I + S0 N) |
---|
| 566 | const double DenQ = (+ Maa*Mbb*Mcc - Maa*Mbc*Mcb - Mab*Mba*Mcc |
---|
| 567 | + Mab*Mbc*Mca + Mac*Mba*Mcb - Mac*Mbb*Mca); |
---|
| 568 | |
---|
| 569 | const double Q11 = ( Mbb*Mcc - Mbc*Mcb)/DenQ; |
---|
| 570 | const double Q12 = (-Mab*Mcc + Mac*Mcb)/DenQ; |
---|
| 571 | const double Q13 = ( Mab*Mbc - Mac*Mbb)/DenQ; |
---|
[13ed84c] | 572 | //const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ; |
---|
[82c299f] | 573 | const double Q22 = ( Maa*Mcc - Mac*Mca)/DenQ; |
---|
| 574 | const double Q23 = (-Maa*Mbc + Mac*Mba)/DenQ; |
---|
[13ed84c] | 575 | //const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ; |
---|
| 576 | //const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ; |
---|
[82c299f] | 577 | const double Q33 = ( Maa*Mbb - Mab*Mba)/DenQ; |
---|
| 578 | |
---|
| 579 | //printf("Q1: %10.5g %10.5g %10.5g\n", Q11, Q12, Q13); |
---|
| 580 | //printf("Q2: %10.5g %10.5g %10.5g\n", Q21, Q22, Q23); |
---|
| 581 | //printf("Q3: %10.5g %10.5g %10.5g\n", Q31, Q32, Q33); |
---|
| 582 | // eq 18c: inv(S) = inv(S0) + mYY^T + \kappa^F |
---|
| 583 | // eq A1 in the appendix |
---|
| 584 | // To solve for S, use: |
---|
| 585 | // S = inv(inv(S^0) + N) inv(S^0) S^0 |
---|
| 586 | // = inv(S^0 inv(S^0) + N) S^0 |
---|
| 587 | // = inv(I + S^0 N) S^0 |
---|
| 588 | // = Q S^0 |
---|
| 589 | const double S11 = Q11*S0aa + Q12*S0ab + Q13*S0ac; |
---|
| 590 | const double S12 = Q12*S0aa + Q22*S0ab + Q23*S0ac; |
---|
| 591 | const double S13 = Q13*S0aa + Q23*S0ab + Q33*S0ac; |
---|
| 592 | const double S22 = Q12*S0ab + Q22*S0bb + Q23*S0bc; |
---|
| 593 | const double S23 = Q13*S0ab + Q23*S0bb + Q33*S0bc; |
---|
| 594 | const double S33 = Q13*S0ac + Q23*S0bc + Q33*S0cc; |
---|
| 595 | // If the full S is needed...it isn't since Ldd = (rho_d - rho_d) = 0 below |
---|
| 596 | //const double S14=-S11-S12-S13; |
---|
| 597 | //const double S24=-S12-S22-S23; |
---|
| 598 | //const double S34=-S13-S23-S33; |
---|
| 599 | //const double S44=S11+S22+S33 + 2.0*(S12+S13+S23); |
---|
| 600 | |
---|
| 601 | // eq 12 of Akcasu, 1990: I(q) = L^T S L |
---|
| 602 | // Note: eliminate cases without A and B polymers by setting Lij to 0 |
---|
| 603 | // Note: 1e-13 to convert from fm to cm for scattering length |
---|
| 604 | const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; |
---|
[d2bb604] | 605 | const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav; |
---|
| 606 | const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav; |
---|
| 607 | const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav; |
---|
[82c299f] | 608 | |
---|
| 609 | const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 |
---|
| 610 | + 2.0*(Lad*Lbd*S12 + Lbd*Lcd*S23 + Lad*Lcd*S13); |
---|
| 611 | |
---|
| 612 | return result; |
---|
[51ec7e8] | 613 | */ |
---|
[82c299f] | 614 | } |
---|