| 1 | #pragma rtGlobals=1 // Use modern global access method. |
|---|
| 2 | #pragma IgorVersion = 6.0 |
|---|
| 3 | |
|---|
| 4 | //CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
|---|
| 5 | //CCCCCCCC |
|---|
| 6 | //C SUBROUTINE FOR THE CALCULATION OF THE SCATTERING FUNCTIONS |
|---|
| 7 | //C OF RODLIKE MICELLES. METHODLOGY FOLLOWS THAT OF PEDERSEN AND |
|---|
| 8 | //C SCHURTENBERGER, MACORMOLECULES, VOL 29,PG 7602, 1996. |
|---|
| 9 | //C WITH EXCULDED VOLUME EFFECTS (METHOD 3) |
|---|
| 10 | // |
|---|
| 11 | // - copied directly from FORTRAN code supplied by Jan Pedersen |
|---|
| 12 | // SRK - 2002, but shows discontinuity at Qlb = 3.1 |
|---|
| 13 | // |
|---|
| 14 | // Jan 2006 - re-worked FORTRAN correcting typos in paper: now is smooth, but |
|---|
| 15 | // the splicing is actually at Qlb = 2, which is not what the paper |
|---|
| 16 | // says is to be done (and not from earlier models) |
|---|
| 17 | // |
|---|
| 18 | // July 2006 - now is CORRECT with Wei-Ren's changes to the code |
|---|
| 19 | // Matlab code was not too difficult to convert to Igor (only a few hours...) |
|---|
| 20 | // |
|---|
| 21 | Proc PlotFlexExclVolCyl(num,qmin,qmax) |
|---|
| 22 | Variable num=128,qmin=0.001,qmax=0.7 |
|---|
| 23 | Prompt num "Enter number of data points for model: " |
|---|
| 24 | Prompt qmin "Enter minimum q-value (A^-1) for model: " |
|---|
| 25 | Prompt qmax "Enter maximum q-value (A^-1) for model: " |
|---|
| 26 | |
|---|
| 27 | Make/O/D/n=(num) xwave_fle,ywave_fle |
|---|
| 28 | xwave_fle = alog(log(qmin) + x*((log(qmax)-log(qmin))/num)) |
|---|
| 29 | Make/O/D coef_fle = {1.,1000,100,20,1.0e-6,6.3e-6,0.0001} |
|---|
| 30 | make/o/t parameters_fle = {"scale","Contour Length (A)","Kuhn Length, b (A)","Radius (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","bkgd (cm^-1)"} |
|---|
| 31 | Edit parameters_fle,coef_fle |
|---|
| 32 | |
|---|
| 33 | Variable/G root:g_fle |
|---|
| 34 | g_fle := FlexExclVolCyl(coef_fle,ywave_fle,xwave_fle) |
|---|
| 35 | Display ywave_fle vs xwave_fle |
|---|
| 36 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
|---|
| 37 | Label bottom "q (A\\S-1\\M)" |
|---|
| 38 | Label left "Intensity (cm\\S-1\\M)" |
|---|
| 39 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
|---|
| 40 | |
|---|
| 41 | AddModelToStrings("FlexExclVolCyl","coef_fle","parameters_fle","fle") |
|---|
| 42 | End |
|---|
| 43 | /////////////////////////////////////////////////////////// |
|---|
| 44 | |
|---|
| 45 | // - sets up a dependency to a wrapper, not the actual SmearedModelFunction |
|---|
| 46 | Proc PlotSmearedFlexExclVolCyl(str) |
|---|
| 47 | String str |
|---|
| 48 | Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) |
|---|
| 49 | |
|---|
| 50 | // if any of the resolution waves are missing => abort |
|---|
| 51 | if(ResolutionWavesMissingDF(str)) //updated to NOT use global strings (in GaussUtils) |
|---|
| 52 | Abort |
|---|
| 53 | endif |
|---|
| 54 | |
|---|
| 55 | SetDataFolder $("root:"+str) |
|---|
| 56 | |
|---|
| 57 | // Setup parameter table for model function |
|---|
| 58 | Make/O/D smear_coef_fle = {1.,1000,100,20,1.0e-6,6.3e-6,0.0001} |
|---|
| 59 | make/o/t smear_parameters_fle = {"scale","Contour Length (A)","Kuhn Length, b (A)","Radius (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","bkgd (cm^-1)"} |
|---|
| 60 | Edit smear_parameters_fle,smear_coef_fle |
|---|
| 61 | |
|---|
| 62 | // output smeared intensity wave, dimensions are identical to experimental QSIG values |
|---|
| 63 | // make extra copy of experimental q-values for easy plotting |
|---|
| 64 | Duplicate/O $(str+"_q") smeared_fle,smeared_qvals |
|---|
| 65 | SetScale d,0,0,"1/cm",smeared_fle |
|---|
| 66 | |
|---|
| 67 | Variable/G gs_fle=0 |
|---|
| 68 | gs_fle := fSmearedFlexExclVolCyl(smear_coef_fle,smeared_fle,smeared_qvals) //this wrapper fills the STRUCT |
|---|
| 69 | |
|---|
| 70 | Display smeared_fle vs smeared_qvals |
|---|
| 71 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
|---|
| 72 | Label bottom "q (A\\S-1\\M)" |
|---|
| 73 | Label left "Intensity (cm\\S-1\\M)" |
|---|
| 74 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
|---|
| 75 | |
|---|
| 76 | SetDataFolder root: |
|---|
| 77 | AddModelToStrings("SmearedFlexExclVolCyl","smear_coef_fle","smear_parameters_fle","fle") |
|---|
| 78 | End |
|---|
| 79 | |
|---|
| 80 | |
|---|
| 81 | //AAO version, uses XOP if available |
|---|
| 82 | // simply calls the original single point calculation with |
|---|
| 83 | // a wave assignment (this will behave nicely if given point ranges) |
|---|
| 84 | Function FlexExclVolCyl(cw,yw,xw) : FitFunc |
|---|
| 85 | Wave cw,yw,xw |
|---|
| 86 | |
|---|
| 87 | #if exists("FlexExclVolCylX") |
|---|
| 88 | yw = FlexExclVolCylX(cw,xw) |
|---|
| 89 | #else |
|---|
| 90 | yw = fFlexExclVolCyl(cw,xw) |
|---|
| 91 | #endif |
|---|
| 92 | return(0) |
|---|
| 93 | End |
|---|
| 94 | |
|---|
| 95 | // |
|---|
| 96 | Function fFlexExclVolCyl(ww,x) |
|---|
| 97 | Wave ww |
|---|
| 98 | Variable x |
|---|
| 99 | |
|---|
| 100 | //nice names to the input params |
|---|
| 101 | //ww[0] = scale |
|---|
| 102 | //ww[1] = L [A] |
|---|
| 103 | //ww[2] = B [A] |
|---|
| 104 | //ww[3] = rad [A] cross-sectional radius |
|---|
| 105 | //ww[4] = sld cylinder |
|---|
| 106 | //ww[5] = sld solvent [A^-2] |
|---|
| 107 | //ww[6] = bkg [cm-1] |
|---|
| 108 | Variable scale,L,B,bkg,rad,qr,cont,sldc,slds |
|---|
| 109 | |
|---|
| 110 | scale = ww[0] |
|---|
| 111 | L = ww[1] |
|---|
| 112 | B = ww[2] |
|---|
| 113 | rad = ww[3] |
|---|
| 114 | sldc = ww[4] |
|---|
| 115 | slds = ww[5] |
|---|
| 116 | bkg = ww[6] |
|---|
| 117 | qr = x*rad //used for cross section contribution only |
|---|
| 118 | |
|---|
| 119 | cont = sldc-slds |
|---|
| 120 | |
|---|
| 121 | Variable flex,crossSect |
|---|
| 122 | flex = Sk_WR(x,L,B) |
|---|
| 123 | |
|---|
| 124 | crossSect = (2*bessJ(1,qr)/qr)^2 |
|---|
| 125 | |
|---|
| 126 | //normalize form factor by multiplying by cylinder volume * cont^2 |
|---|
| 127 | // then convert to cm-1 by multiplying by 10^8 |
|---|
| 128 | // then scale = phi |
|---|
| 129 | |
|---|
| 130 | flex *= crossSect |
|---|
| 131 | flex *= Pi*rad*rad*L |
|---|
| 132 | flex *= cont^2 |
|---|
| 133 | flex *= 1.0e8 |
|---|
| 134 | |
|---|
| 135 | return (scale*flex + bkg) |
|---|
| 136 | End |
|---|
| 137 | |
|---|
| 138 | //////////////////WRC corrected code below |
|---|
| 139 | // main function |
|---|
| 140 | function Sk_WR(q,L,b) |
|---|
| 141 | Variable q,L,b |
|---|
| 142 | // |
|---|
| 143 | Variable p1,p2,p1short,p2short,q0,qconnect |
|---|
| 144 | Variable c,epsilon,ans,q0short,Sexvmodify |
|---|
| 145 | |
|---|
| 146 | p1 = 4.12 |
|---|
| 147 | p2 = 4.42 |
|---|
| 148 | p1short = 5.36 |
|---|
| 149 | p2short = 5.62 |
|---|
| 150 | q0 = 3.1 |
|---|
| 151 | qconnect = q0/b |
|---|
| 152 | // |
|---|
| 153 | q0short = max(1.9/sqrt(Rgsquareshort(q,L,b)),3) |
|---|
| 154 | |
|---|
| 155 | // |
|---|
| 156 | if(L/b > 10) |
|---|
| 157 | C = 3.06/(L/b)^0.44 |
|---|
| 158 | epsilon = 0.176 |
|---|
| 159 | else |
|---|
| 160 | C = 1 |
|---|
| 161 | epsilon = 0.170 |
|---|
| 162 | endif |
|---|
| 163 | // |
|---|
| 164 | |
|---|
| 165 | if( L > 4*b ) // Longer Chains |
|---|
| 166 | if (q*b <= 3.1) |
|---|
| 167 | //Modified by Yun on Oct. 15, |
|---|
| 168 | Sexvmodify = Sexvnew(q, L, b) |
|---|
| 169 | ans = Sexvmodify + C * (4/15 + 7./(15*u_WR(q,L,b)) - (11/15 + 7./(15*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L) |
|---|
| 170 | else //q(i)*b > 3.1 |
|---|
| 171 | ans = a1long(q, L, b, p1, p2, q0)/((q*b)^p1) + a2long(q, L, b, p1, p2, q0)/((q*b)^p2) + pi/(q*L) |
|---|
| 172 | endif |
|---|
| 173 | else //L <= 4*b Shorter Chains |
|---|
| 174 | if (q*b <= max(1.9/sqrt(Rgsquareshort(q,L,b)),3) ) |
|---|
| 175 | if (q*b<=0.01) |
|---|
| 176 | ans = 1 - Rgsquareshort(q,L,b)*(q^2)/3 |
|---|
| 177 | else |
|---|
| 178 | ans = Sdebye1(q,L,b) |
|---|
| 179 | endif |
|---|
| 180 | else //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) |
|---|
| 181 | ans = a1short(q,L,b,p1short,p2short,q0short)/((q*b)^p1short) + a2short(q,L,b,p1short,p2short,q0short)/((q*b)^p2short) + pi/(q*L) |
|---|
| 182 | endif |
|---|
| 183 | endif |
|---|
| 184 | |
|---|
| 185 | return(ans) |
|---|
| 186 | end |
|---|
| 187 | |
|---|
| 188 | //WR named this w (too generic) |
|---|
| 189 | Function w_WR(x) |
|---|
| 190 | Variable x |
|---|
| 191 | |
|---|
| 192 | //C4 = 1.523; |
|---|
| 193 | //C5 = 0.1477; |
|---|
| 194 | Variable yy |
|---|
| 195 | yy = 0.5*(1 + tanh((x - 1.523)/0.1477)) |
|---|
| 196 | |
|---|
| 197 | return (yy) |
|---|
| 198 | end |
|---|
| 199 | |
|---|
| 200 | // |
|---|
| 201 | function u1(q,L,b) |
|---|
| 202 | Variable q,L,b |
|---|
| 203 | Variable yy |
|---|
| 204 | |
|---|
| 205 | yy = Rgsquareshort(q,L,b)*q^2 |
|---|
| 206 | |
|---|
| 207 | return yy |
|---|
| 208 | end |
|---|
| 209 | |
|---|
| 210 | // was named u |
|---|
| 211 | function u_WR(q,L,b) |
|---|
| 212 | Variable q,L,b |
|---|
| 213 | |
|---|
| 214 | Variable yy |
|---|
| 215 | yy = Rgsquare(q,L,b)*(q^2) |
|---|
| 216 | return yy |
|---|
| 217 | end |
|---|
| 218 | |
|---|
| 219 | |
|---|
| 220 | |
|---|
| 221 | // |
|---|
| 222 | function Rgsquarezero(q,L,b) |
|---|
| 223 | Variable q,L,b |
|---|
| 224 | |
|---|
| 225 | Variable yy |
|---|
| 226 | yy = (L*b/6) * (1 - 1.5*(b/L) + 1.5*(b/L)^2 - 0.75*(b/L)^3*(1 - exp(-2*(L/b)))) |
|---|
| 227 | |
|---|
| 228 | return yy |
|---|
| 229 | end |
|---|
| 230 | |
|---|
| 231 | // |
|---|
| 232 | function Rgsquareshort(q,L,b) |
|---|
| 233 | Variable q,L,b |
|---|
| 234 | |
|---|
| 235 | Variable yy |
|---|
| 236 | yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b) |
|---|
| 237 | |
|---|
| 238 | return yy |
|---|
| 239 | end |
|---|
| 240 | |
|---|
| 241 | // |
|---|
| 242 | function Rgsquare(q,L,b) |
|---|
| 243 | Variable q,L,b |
|---|
| 244 | |
|---|
| 245 | Variable yy |
|---|
| 246 | yy = AlphaSquare(L/b)*L*b/6 |
|---|
| 247 | |
|---|
| 248 | return yy |
|---|
| 249 | end |
|---|
| 250 | |
|---|
| 251 | // |
|---|
| 252 | function AlphaSquare(x) |
|---|
| 253 | Variable x |
|---|
| 254 | |
|---|
| 255 | Variable yy |
|---|
| 256 | yy = (1 + (x/3.12)^2 + (x/8.67)^3)^(0.176/3) |
|---|
| 257 | |
|---|
| 258 | return yy |
|---|
| 259 | end |
|---|
| 260 | |
|---|
| 261 | // |
|---|
| 262 | function miu(x) |
|---|
| 263 | Variable x |
|---|
| 264 | |
|---|
| 265 | Variable yy |
|---|
| 266 | yy = (1/8)*(9*x - 2 + 2*log(1 + x)/x)*exp(1/2.565*(1/x + (1 - 1/(x^2))*log(1 + x))) |
|---|
| 267 | |
|---|
| 268 | return yy |
|---|
| 269 | end |
|---|
| 270 | |
|---|
| 271 | // |
|---|
| 272 | function Sdebye(q,L,b) |
|---|
| 273 | Variable q,L,b |
|---|
| 274 | |
|---|
| 275 | Variable yy |
|---|
| 276 | yy = 2*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1)/((u_WR(q,L,b))^2) |
|---|
| 277 | |
|---|
| 278 | return yy |
|---|
| 279 | end |
|---|
| 280 | |
|---|
| 281 | // |
|---|
| 282 | function Sdebye1(q,L,b) |
|---|
| 283 | Variable q,L,b |
|---|
| 284 | |
|---|
| 285 | Variable yy |
|---|
| 286 | yy = 2*(exp(-u1(q,L,b)) + u1(q,L,b) -1)/((u1(q,L,b))^2) |
|---|
| 287 | |
|---|
| 288 | return yy |
|---|
| 289 | end |
|---|
| 290 | |
|---|
| 291 | // |
|---|
| 292 | function Sexv(q,L,b) |
|---|
| 293 | Variable q,L,b |
|---|
| 294 | |
|---|
| 295 | Variable yy,C1,C2,C3,miu,Rg2 |
|---|
| 296 | C1=1.22 |
|---|
| 297 | C2=0.4288 |
|---|
| 298 | C3=-1.651 |
|---|
| 299 | miu = 0.585 |
|---|
| 300 | |
|---|
| 301 | Rg2 = Rgsquare(q,L,b) |
|---|
| 302 | |
|---|
| 303 | yy = (1 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*(q*sqrt(Rg2))^(-1/miu) + C2*(q*sqrt(Rg2))^(-2/miu) + C3*(q*sqrt(Rg2))^(-3/miu)) |
|---|
| 304 | |
|---|
| 305 | return yy |
|---|
| 306 | end |
|---|
| 307 | |
|---|
| 308 | // this must be WR modified version |
|---|
| 309 | function Sexvnew(q,L,b) |
|---|
| 310 | Variable q,L,b |
|---|
| 311 | |
|---|
| 312 | Variable yy,C1,C2,C3,miu |
|---|
| 313 | C1=1.22 |
|---|
| 314 | C2=0.4288 |
|---|
| 315 | C3=-1.651 |
|---|
| 316 | miu = 0.585 |
|---|
| 317 | |
|---|
| 318 | //calculating the derivative to decide on the corection (cutoff) term? |
|---|
| 319 | // I have modified this from WRs original code |
|---|
| 320 | Variable del=1.05,C_star2,Rg2 |
|---|
| 321 | if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0 ) |
|---|
| 322 | C_star2 = 0 |
|---|
| 323 | else |
|---|
| 324 | C_star2 = 1 |
|---|
| 325 | endif |
|---|
| 326 | |
|---|
| 327 | Rg2 = Rgsquare(q,L,b) |
|---|
| 328 | |
|---|
| 329 | yy = (1 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*(q*sqrt(Rg2))^(-1/miu) + C2*(q*sqrt(Rg2))^(-2/miu) + C3*(q*sqrt(Rg2))^(-3/miu)) |
|---|
| 330 | |
|---|
| 331 | return yy |
|---|
| 332 | end |
|---|
| 333 | |
|---|
| 334 | |
|---|
| 335 | |
|---|
| 336 | // these are the messy ones |
|---|
| 337 | function a2short(q, L, b, p1short, p2short, q0) |
|---|
| 338 | Variable q, L, b, p1short, p2short, q0 |
|---|
| 339 | |
|---|
| 340 | Variable yy,Rg2_sh |
|---|
| 341 | Rg2_sh = Rgsquareshort(q,L,b) |
|---|
| 342 | |
|---|
| 343 | Variable t1 |
|---|
| 344 | t1 = ((q0^2*Rg2_sh)/b^2) |
|---|
| 345 | |
|---|
| 346 | //E is the number e |
|---|
| 347 | yy = ((-(1/(L*((p1short - p2short))*Rg2_sh^2)*((b*E^(-t1)*q0^(-4 + p2short)*((8*b^3*L - 8*b^3*E^t1*L - 2*b^3*L*p1short + 2*b^3*E^t1*L*p1short + 4*b*L*q0^2*Rg2_sh + 4*b*E^t1*L*q0^2*Rg2_sh - 2*b*E^t1*L*p1short*q0^2*Rg2_sh - E^t1*pi*q0^3*Rg2_sh^2 + E^t1*p1short*pi*q0^3*Rg2_sh^2))))))) |
|---|
| 348 | |
|---|
| 349 | return yy |
|---|
| 350 | end |
|---|
| 351 | |
|---|
| 352 | // |
|---|
| 353 | function a1short(q, L, b, p1short, p2short, q0) |
|---|
| 354 | Variable q, L, b, p1short, p2short, q0 |
|---|
| 355 | |
|---|
| 356 | Variable yy,Rg2_sh |
|---|
| 357 | Rg2_sh = Rgsquareshort(q,L,b) |
|---|
| 358 | |
|---|
| 359 | Variable t1 |
|---|
| 360 | t1 = ((q0^2*Rg2_sh)/b^2) |
|---|
| 361 | |
|---|
| 362 | yy = ((1/(L*((p1short - p2short))*Rg2_sh^2)*((b*E^(-t1)*q0^((-4) + p1short)*((8*b^3*L - 8*b^3*E^t1*L - 2*b^3*L*p2short + 2*b^3*E^t1*L*p2short + 4*b*L*q0^2*Rg2_sh + 4*b*E^t1*L*q0^2*Rg2_sh - 2*b*E^t1*L*p2short*q0^2*Rg2_sh - E^t1*pi*q0^3*Rg2_sh^2 + E^t1*p2short*pi*q0^3*Rg2_sh^2)))))) |
|---|
| 363 | |
|---|
| 364 | return yy |
|---|
| 365 | end |
|---|
| 366 | |
|---|
| 367 | // this one will be lots of trouble |
|---|
| 368 | function a2long(q, L, b, p1, p2, q0) |
|---|
| 369 | variable q, L, b, p1, p2, q0 |
|---|
| 370 | |
|---|
| 371 | Variable yy,c1,c2,c3,c4,c5,miu,c,Rg2,rRg |
|---|
| 372 | Variable t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13 |
|---|
| 373 | |
|---|
| 374 | if( L/b > 10) |
|---|
| 375 | C = 3.06/(L/b)^0.44 |
|---|
| 376 | else |
|---|
| 377 | C = 1 |
|---|
| 378 | endif |
|---|
| 379 | |
|---|
| 380 | C1 = 1.22 |
|---|
| 381 | C2 = 0.4288 |
|---|
| 382 | C3 = -1.651 |
|---|
| 383 | C4 = 1.523 |
|---|
| 384 | C5 = 0.1477 |
|---|
| 385 | miu = 0.585 |
|---|
| 386 | |
|---|
| 387 | Rg2 = Rgsquare(q,L,b) |
|---|
| 388 | t1 = (1/(b* p1*q0^((-1) - p1 - p2) - b*p2*q0^((-1) - p1 - p2))) |
|---|
| 389 | t2 = (b*C*(((-1*((14*b^3)/(15*q0^3*Rg2))) + (14*b^3*E^(-((q0^2*Rg2)/b^2)))/(15*q0^3*Rg2) + (2*E^(-((q0^2*Rg2)/b^2))*q0*((11/15 + (7*b^2)/(15*q0^2*Rg2)))*Rg2)/b)))/L |
|---|
| 390 | t3 = (sqrt(Rg2)*((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))*sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2)/(2*C5) |
|---|
| 391 | t4 = (b^4*sqrt(Rg2)*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2)/(C5*q0^4*Rg2^2) |
|---|
| 392 | t5 = (2*b^4*(((2*q0*Rg2)/b - (2*E^(-((q0^2*Rg2)/b^2))*q0*Rg2)/b))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q0^4*Rg2^2) |
|---|
| 393 | t6 = (8*b^5*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q0^5*Rg2^2) |
|---|
| 394 | t7 = (((-((3*C3*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 3/miu))/miu)) - (2*C2*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 2/miu))/miu - (C1*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 1/miu))/miu)) |
|---|
| 395 | t8 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))) |
|---|
| 396 | t9 = (b*C*((4/15 - E^(-((q0^2*Rg2)/b^2))*((11/15 + (7*b^2)/(15*q0^2*Rg2))) + (7*b^2)/(15*q0^2*Rg2))))/L |
|---|
| 397 | t10 = (2*b^4*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q0^4*Rg2^2) |
|---|
| 398 | |
|---|
| 399 | |
|---|
| 400 | yy = ((-1*(t1* (((-q0^(-p1))*(((b^2*pi)/(L*q0^2) + t2 + t3 - t4 + t5 - t6 + 1/2*t7*t8)) - b*p1*q0^((-1) - p1)*(((-((b*pi)/(L*q0))) + t9 + t10 + 1/2*((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))) |
|---|
| 401 | |
|---|
| 402 | |
|---|
| 403 | return yy |
|---|
| 404 | end |
|---|
| 405 | |
|---|
| 406 | //need to define this on my own |
|---|
| 407 | Function sech_WR(x) |
|---|
| 408 | variable x |
|---|
| 409 | |
|---|
| 410 | return(1/cosh(x)) |
|---|
| 411 | end |
|---|
| 412 | // |
|---|
| 413 | function a1long(q, L, b, p1, p2, q0) |
|---|
| 414 | Variable q, L, b, p1, p2, q0 |
|---|
| 415 | |
|---|
| 416 | Variable yy,c,c1,c2,c3,c4,c5,miu,Rg2,rRg |
|---|
| 417 | Variable t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16 |
|---|
| 418 | |
|---|
| 419 | if( L/b > 10) |
|---|
| 420 | C = 3.06/(L/b)^0.44 |
|---|
| 421 | else |
|---|
| 422 | C = 1 |
|---|
| 423 | endif |
|---|
| 424 | |
|---|
| 425 | C1 = 1.22 |
|---|
| 426 | C2 = 0.4288 |
|---|
| 427 | C3 = -1.651 |
|---|
| 428 | C4 = 1.523 |
|---|
| 429 | C5 = 0.1477 |
|---|
| 430 | miu = 0.585 |
|---|
| 431 | |
|---|
| 432 | Rg2 = Rgsquare(q,L,b) |
|---|
| 433 | t1 = (b*C*((4/15 - E^(-((q0^2*Rg2)/b^2))*((11/15 + (7*b^2)/(15*q0^2*Rg2))) + (7*b^2)/(15*q0^2*Rg2)))) |
|---|
| 434 | t2 = (2*b^4*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))) |
|---|
| 435 | t3 = ((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu))) |
|---|
| 436 | t4 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))) |
|---|
| 437 | t5 = (1/(b*p1*q0^((-1) - p1 - p2) - b*p2*q0^((-1) - p1 - p2))) |
|---|
| 438 | t6 = (b*C*(((-((14*b^3)/(15*q0^3*Rg2))) + (14*b^3*E^(-((q0^2*Rg2)/b^2)))/(15*q0^3*Rg2) + (2*E^(-((q0^2*Rg2)/b^2))*q0*((11/15 + (7*b^2)/(15*q0^2*Rg2)))*Rg2)/b))) |
|---|
| 439 | t7 = (sqrt(Rg2)*((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))*sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2) |
|---|
| 440 | t8 = (b^4*sqrt(Rg2)*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2) |
|---|
| 441 | t9 = (2*b^4*(((2*q0*Rg2)/b - (2*E^(-((q0^2*Rg2)/b^2))*q0*Rg2)/b))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))) |
|---|
| 442 | t10 = (8*b^5*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))) |
|---|
| 443 | t11 = (((-((3*C3*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 3/miu))/miu)) - (2*C2*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 2/miu))/miu - (C1*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 1/miu))/miu)) |
|---|
| 444 | t12 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))) |
|---|
| 445 | t13 = (b*C*((4/15 - E^(-((q0^2*Rg2)/b^2))*((11/15 + (7*b^2)/(15*q0^2* Rg2))) + (7*b^2)/(15*q0^2*Rg2)))) |
|---|
| 446 | t14 = (2*b^4*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))) |
|---|
| 447 | t15 = ((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu))) |
|---|
| 448 | |
|---|
| 449 | |
|---|
| 450 | yy = (q0^p1*(((-((b*pi)/(L*q0))) +t1/L +t2/(q0^4*Rg2^2) + 1/2*t3*t4)) + (t5*((q0^(p1 - p2)*(((-q0^(-p1))*(((b^2*pi)/(L*q0^2) +t6/L +t7/(2*C5) -t8/(C5*q0^4*Rg2^2) +t9/(q0^4*Rg2^2) -t10/(q0^5*Rg2^2) + 1/2*t11*t12)) - b*p1*q0^((-1) - p1)*(((-((b*pi)/(L*q0))) +t13/L +t14/(q0^4*Rg2^2) + 1/2*t15*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))))) |
|---|
| 451 | |
|---|
| 452 | |
|---|
| 453 | return yy |
|---|
| 454 | end |
|---|
| 455 | |
|---|
| 456 | |
|---|
| 457 | // unused functions copied from WRC Matlab code |
|---|
| 458 | |
|---|
| 459 | |
|---|
| 460 | //function Vmic(phi) |
|---|
| 461 | // Variable phi |
|---|
| 462 | // |
|---|
| 463 | // Variable yy |
|---|
| 464 | // yy = 2.53*10^6*(phi)^(0.5) + 8.35*(phi)^(-2) |
|---|
| 465 | // |
|---|
| 466 | // return yy |
|---|
| 467 | //end |
|---|
| 468 | |
|---|
| 469 | //function Scs(x) |
|---|
| 470 | // Variable x |
|---|
| 471 | // |
|---|
| 472 | // Variable yy |
|---|
| 473 | // yy = (2*bessj(1,x)/x)^2 |
|---|
| 474 | // |
|---|
| 475 | // return yy |
|---|
| 476 | //end |
|---|
| 477 | |
|---|
| 478 | //function Srod(q,L,b) |
|---|
| 479 | // Variable q,L,b |
|---|
| 480 | // |
|---|
| 481 | // Variable yy |
|---|
| 482 | // yy = 2 * Si(q,L,b)/(q*L) - 4 * (sin(q*L/2)^2)/((q*L)^2) |
|---|
| 483 | // |
|---|
| 484 | // return yy |
|---|
| 485 | //end |
|---|
| 486 | // |
|---|
| 487 | //function Si(q,L,b) |
|---|
| 488 | // Variable q,L,b |
|---|
| 489 | // |
|---|
| 490 | // Variable yy |
|---|
| 491 | // |
|---|
| 492 | //// for i=1:length(q) |
|---|
| 493 | //// y(i) = quadl('sin(x)./(x)',10e-8,q(i)*L); |
|---|
| 494 | //// end |
|---|
| 495 | // |
|---|
| 496 | // return yy |
|---|
| 497 | //end |
|---|
| 498 | |
|---|
| 499 | ////////////////// |
|---|
| 500 | |
|---|
| 501 | //wrapper to calculate the smeared model as an AAO-Struct |
|---|
| 502 | // fills the struct and calls the ususal function with the STRUCT parameter |
|---|
| 503 | // |
|---|
| 504 | // used only for the dependency, not for fitting |
|---|
| 505 | // |
|---|
| 506 | Function fSmearedFlexExclVolCyl(coefW,yW,xW) |
|---|
| 507 | Wave coefW,yW,xW |
|---|
| 508 | |
|---|
| 509 | String str = getWavesDataFolder(yW,0) |
|---|
| 510 | String DF="root:"+str+":" |
|---|
| 511 | |
|---|
| 512 | WAVE resW = $(DF+str+"_res") |
|---|
| 513 | |
|---|
| 514 | STRUCT ResSmearAAOStruct fs |
|---|
| 515 | WAVE fs.coefW = coefW |
|---|
| 516 | WAVE fs.yW = yW |
|---|
| 517 | WAVE fs.xW = xW |
|---|
| 518 | WAVE fs.resW = resW |
|---|
| 519 | |
|---|
| 520 | Variable err |
|---|
| 521 | err = SmearedFlexExclVolCyl(fs) |
|---|
| 522 | |
|---|
| 523 | return (0) |
|---|
| 524 | End |
|---|
| 525 | |
|---|
| 526 | // this is all there is to the smeared calculation! |
|---|
| 527 | Function SmearedFlexExclVolCyl(s) :FitFunc |
|---|
| 528 | Struct ResSmearAAOStruct &s |
|---|
| 529 | |
|---|
| 530 | // the name of your unsmeared model (AAO) is the first argument |
|---|
| 531 | Smear_Model_20(FlexExclVolCyl,s.coefW,s.xW,s.yW,s.resW) |
|---|
| 532 | |
|---|
| 533 | return(0) |
|---|
| 534 | End |
|---|