| 1 | #pragma rtGlobals=1 // Use modern global access method. |
|---|
| 2 | #pragma IgorVersion = 6.0 |
|---|
| 3 | |
|---|
| 4 | //////////////////////////////////////////////////// |
|---|
| 5 | // |
|---|
| 6 | // description |
|---|
| 7 | // |
|---|
| 8 | // |
|---|
| 9 | // reference - |
|---|
| 10 | // - errors in paper |
|---|
| 11 | // |
|---|
| 12 | // warning about unphysical regions |
|---|
| 13 | // |
|---|
| 14 | // + how to properly handle if result is unphysical???? |
|---|
| 15 | // |
|---|
| 16 | // |
|---|
| 17 | // S. Kline 15 JUL 2004 |
|---|
| 18 | // |
|---|
| 19 | //////////////////////////////////////////////////// |
|---|
| 20 | |
|---|
| 21 | //this macro sets up all the necessary parameters and waves that are |
|---|
| 22 | //needed to calculate the model function. |
|---|
| 23 | // |
|---|
| 24 | Proc PlotStickyHS_Struct(num,qmin,qmax) |
|---|
| 25 | Variable num=256, qmin=.001, qmax=.5 |
|---|
| 26 | Prompt num "Enter number of data points for model: " |
|---|
| 27 | Prompt qmin "Enter minimum q-value (^1) for model: " |
|---|
| 28 | Prompt qmax "Enter maximum q-value (^1) for model: " |
|---|
| 29 | // |
|---|
| 30 | Make/O/D/n=(num) xwave_shsSQ, ywave_shsSQ |
|---|
| 31 | xwave_shsSQ = alog(log(qmin) + x*((log(qmax)-log(qmin))/num)) |
|---|
| 32 | Make/O/D coef_shsSQ = {50,0.1,0.05,0.2} //CH#2 |
|---|
| 33 | make/o/t parameters_shsSQ = {"Radius","volume fraction","perturbation parameter (0.1)","stickiness, tau"} //CH#3 |
|---|
| 34 | Edit parameters_shsSQ, coef_shsSQ |
|---|
| 35 | Variable/G root:g_shsSQ |
|---|
| 36 | g_shsSQ := StickyHS_Struct(coef_shsSQ, ywave_shsSQ, xwave_shsSQ) |
|---|
| 37 | // ywave_shsSQ := StickyHS_Struct(coef_shsSQ, xwave_shsSQ) |
|---|
| 38 | Display ywave_shsSQ vs xwave_shsSQ |
|---|
| 39 | ModifyGraph marker=29, msize=2, mode=4 |
|---|
| 40 | ModifyGraph log=0 |
|---|
| 41 | Label bottom "q (\\S-1\\M) " |
|---|
| 42 | Label left "S(q)" |
|---|
| 43 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
|---|
| 44 | |
|---|
| 45 | AddModelToStrings("StickyHS_Struct","coef_shsSQ","shsSQ") |
|---|
| 46 | // |
|---|
| 47 | End |
|---|
| 48 | |
|---|
| 49 | //AAO function |
|---|
| 50 | Function StickyHS_Struct(cw,yw,xw) : FitFunc |
|---|
| 51 | Wave cw,yw,xw |
|---|
| 52 | |
|---|
| 53 | #if exists("StickyHS_StructX") |
|---|
| 54 | yw = StickyHS_StructX(cw,xw) |
|---|
| 55 | #else |
|---|
| 56 | yw = fStickyHS_Struct(cw,xw) |
|---|
| 57 | #endif |
|---|
| 58 | return(0) |
|---|
| 59 | End |
|---|
| 60 | |
|---|
| 61 | Function fStickyHS_Struct(w,x) : FitFunc |
|---|
| 62 | Wave w |
|---|
| 63 | Variable x |
|---|
| 64 | // Input (fitting) variables are: |
|---|
| 65 | //radius = w[0] |
|---|
| 66 | //volume fraction = w[1] |
|---|
| 67 | //epsilon (perurbation param) = w[2] |
|---|
| 68 | //tau (stickiness) = w[3] |
|---|
| 69 | Variable rad,phi,eps,tau,eta |
|---|
| 70 | Variable sig,aa,etam1,qa,qb,qc,radic |
|---|
| 71 | Variable lam,lam2,test,mu,alpha,beta |
|---|
| 72 | Variable qv,kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq |
|---|
| 73 | rad = w[0] |
|---|
| 74 | phi = w[1] |
|---|
| 75 | eps = w[2] |
|---|
| 76 | tau = w[3] |
|---|
| 77 | |
|---|
| 78 | eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps) |
|---|
| 79 | |
|---|
| 80 | sig = 2.0 * rad |
|---|
| 81 | aa = sig/(1.0 - eps) |
|---|
| 82 | etam1 = 1.0 - eta |
|---|
| 83 | //C |
|---|
| 84 | //C SOLVE QUADRATIC FOR LAMBDA |
|---|
| 85 | //C |
|---|
| 86 | qa = eta/12.0 |
|---|
| 87 | qb = -1.0*(tau + eta/(etam1)) |
|---|
| 88 | qc = (1.0 + eta/2.0)/(etam1*etam1) |
|---|
| 89 | radic = qb*qb - 4.0*qa*qc |
|---|
| 90 | if(radic<0) |
|---|
| 91 | if(x>0.01 && x<0.015) |
|---|
| 92 | Print "Lambda unphysical - both roots imaginary" |
|---|
| 93 | endif |
|---|
| 94 | return(-1) |
|---|
| 95 | endif |
|---|
| 96 | //C KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL |
|---|
| 97 | lam = (-1.0*qb-sqrt(radic))/(2.0*qa) |
|---|
| 98 | lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa) |
|---|
| 99 | if(lam2<lam) |
|---|
| 100 | lam = lam2 |
|---|
| 101 | endif |
|---|
| 102 | test = 1.0 + 2.0*eta |
|---|
| 103 | mu = lam*eta*etam1 |
|---|
| 104 | if(mu>test) |
|---|
| 105 | if(x>0.01 && x<0.015) |
|---|
| 106 | Print "Lambda unphysical mu>test" |
|---|
| 107 | endif |
|---|
| 108 | return(-1) |
|---|
| 109 | endif |
|---|
| 110 | alpha = (1.0 + 2.0*eta - mu)/(etam1*etam1) |
|---|
| 111 | beta = (mu - 3.0*eta)/(2.0*etam1*etam1) |
|---|
| 112 | |
|---|
| 113 | //C |
|---|
| 114 | //C CALCULATE THE STRUCTURE FACTOR |
|---|
| 115 | //C |
|---|
| 116 | |
|---|
| 117 | qv = x |
|---|
| 118 | kk = qv*aa |
|---|
| 119 | k2 = kk*kk |
|---|
| 120 | k3 = kk*k2 |
|---|
| 121 | ds = sin(kk) |
|---|
| 122 | dc = cos(kk) |
|---|
| 123 | aq1 = ((ds - kk*dc)*alpha)/k3 |
|---|
| 124 | aq2 = (beta*(1.0-dc))/k2 |
|---|
| 125 | aq3 = (lam*ds)/(12.0*kk) |
|---|
| 126 | aq = 1.0 + 12.0*eta*(aq1+aq2-aq3) |
|---|
| 127 | // |
|---|
| 128 | bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3) |
|---|
| 129 | bq2 = beta*(1.0/kk - ds/k2) |
|---|
| 130 | bq3 = (lam/12.0)*((1.0 - dc)/kk) |
|---|
| 131 | bq = 12.0*eta*(bq1+bq2-bq3) |
|---|
| 132 | // |
|---|
| 133 | sq = 1.0/(aq*aq +bq*bq) |
|---|
| 134 | |
|---|
| 135 | Return (sq) |
|---|
| 136 | End |
|---|
| 137 | |
|---|
