source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/FlexCyl_PolyLen.ipf @ 131

Last change on this file since 131 was 131, checked in by srkline, 15 years ago

Typo in dialog in every model...

File size: 6.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4#include "FlexibleCylinder"
5//uses the function FlexibleCylinder(.ipf) as basic function
6//
7// code has been updated with WRC's changes (located in FlexibleCylinder.ipf)
8// JULY 2006
9//
10Proc PlotFlexCyl_PolyLength(num,qmin,qmax)
11        Variable num=128,qmin=0.001,qmax=0.7
12        Prompt num "Enter number of data points for model: "
13        Prompt qmin "Enter minimum q-value (^-1) for model: "
14        Prompt qmax "Enter maximum q-value (^-1) for model: "
15       
16        // Setup parameter table for model function
17        Make/O/D/n=(num) xwave_flepl,ywave_flepl
18        xwave_flepl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
19        Make/O/D coef_flepl = {1.,1000,0.001,100,20,3.0e-6,0.0001}
20        make/o/t parameters_flepl =  {"scale","Contour Length (A)","polydispersity of Contour Length","Kuhn Length, b (A)","Radius (A)","contrast (A^-2)","bkgd (cm^-1)"}
21        Edit parameters_flepl,coef_flepl
22       
23        Variable/G root:g_flepl
24        g_flepl := FlexCyl_PolyLen(coef_flepl,ywave_flepl,xwave_flepl)
25        Display ywave_flepl vs xwave_flepl
26        ModifyGraph log=1,marker=29,msize=2,mode=4
27        Label bottom "q (\\S-1\\M)"
28        Label left "Intensity (cm\\S-1\\M)"
29        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
30End
31
32// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
33Proc PlotSmeared_FlexCyl_PolyLen(str)                                                           
34        String str
35        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
36       
37        // if any of the resolution waves are missing => abort
38        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
39                Abort
40        endif
41       
42        SetDataFolder $("root:"+str)
43       
44        // Setup parameter table for model function
45        Make/O/D smear_coef_flepl = {1.,1000,0.001,100,20,3.0e-6,0.0001}                //CH#4
46        make/o/t smear_parameters_flepl = {"scale","Contour Length (A)","polydispersity of Contour Length","Kuhn Length, b (A)","Radius (A)","contrast (A^-2)","bkgd (cm^-1)"}
47        Edit smear_parameters_flepl,smear_coef_flepl                                    //display parameters in a table
48       
49        // output smeared intensity wave, dimensions are identical to experimental QSIG values
50        // make extra copy of experimental q-values for easy plotting
51        Duplicate/O $(str+"_q") smeared_flepl,smeared_qvals                             //
52        SetScale d,0,0,"1/cm",smeared_flepl                                                     //
53                                       
54        Variable/G gs_flepl=0
55        gs_flepl := fSmearedFlexCyl_PolyLen(smear_coef_flepl,smeared_flepl,smeared_qvals)       //this wrapper fills the STRUCT
56       
57        Display smeared_flepl vs smeared_qvals                                                                  //
58        ModifyGraph log=1,marker=29,msize=2,mode=4
59        Label bottom "q (\\S-1\\M)"
60        Label left "I(q) (cm\\S-1\\M)"
61        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
62       
63        SetDataFolder root:
64End
65       
66
67Function Schulz_Point_flepl(x,avg,zz)
68
69        //Wave w
70        Variable x,avg,zz
71        Variable dr
72       
73        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
74       
75        return (exp(dr))
76       
77End
78
79
80//AAO version, uses XOP if available
81// simply calls the original single point calculation with
82// a wave assignment (this will behave nicely if given point ranges)
83Function FlexCyl_PolyLen(cw,yw,xw) : FitFunc
84        Wave cw,yw,xw
85       
86#if exists("FlexCyl_PolyLenX")
87        yw = FlexCyl_PolyLenX(cw,xw)
88#else
89        yw = fFlexCyl_PolyLen(cw,xw)
90#endif
91        return(0)
92End
93
94Function fFlexCyl_PolyLen(w,x) : FitFunc
95        Wave w
96        Variable x
97     
98        Variable scale,radius,pd,delrho,bkg,zz,length,lb
99        scale = w[0]
100        length = w[1]
101        pd = w[2]
102        lb = w[3]
103        radius = w[4]
104        delrho = w[5]
105        bkg = w[6]
106       
107        zz = (1/pd)^2-1
108//
109// the OUTPUT form factor is <f^2>/Vavg [cm-1]
110//
111// local variables
112        Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
113        Variable answer,zp1,zp2,zp3,vpoly
114        String weightStr,zStr
115       
116        nord = 20
117        weightStr = "gauss20wt"
118        zStr = "gauss20z"
119
120// use 20 Gauss points
121        if (WaveExists($weightStr) == 0) // wave reference is not valid,
122                Make/D/N=(nord) $weightStr,$zStr
123                Wave wtGau = $weightStr
124                Wave zGau = $zStr               // wave references to pass
125                Make20GaussPoints(wtGau,zGau)   
126        else
127                if(exists(weightStr) > 1)
128                         Abort "wave name is already in use"    // execute if condition is false
129                endif
130                Wave wtGau = $weightStr
131                Wave zGau = $zStr
132        endif
133
134// set up the integration
135// end points and weights
136// limits are technically 0-inf, but wisely choose non-zero region of distribution
137        Variable range=3.4              //multiples of the std. dev. fom the mean
138        a = length*(1-range*pd)
139        if (a<0)
140                a=0             //otherwise numerical error when pd >= 0.3, making a<0
141        endif
142        If(pd>0.3)
143                range = 3.4 + (pd-0.3)*18
144        Endif
145        b = length*(1+range*pd) // is this far enough past avg length?
146        va =a
147        vb =b
148
149// evaluate at Gauss points
150        // remember to index from 0,size-1     
151        qq = x          //current x point is the q-value for evaluation
152        summ = 0.0              // initialize integral
153   ii=0
154   do
155                zi = ( zGau[ii]*(vb-va) + vb + va )/2.0         
156                yyy = wtGau[ii] * fle_kernel(qq,radius,length,lb,zz,delrho,zi)
157                summ = yyy + summ
158                ii+=1
159        while (ii<nord)                         // end of loop over quadrature points
160//   
161// calculate value of integral to return
162   answer = (vb-va)/2.0*summ
163     
164//  contrast^2 is included in integration rad_kernel
165//      answer *= delrho*delrho
166//normalize by polydisperse volume
167// now volume depends on polydisperse Length - so normalize by the FIRST moment
168// 1st moment = volume!
169        vpoly = Pi*(radius)^2*length
170//Divide by vol, since volume has been "un-normalized" out
171        answer /= vpoly
172//convert to [cm-1]
173        answer *= 1.0e8
174//scale
175        answer *= scale
176// add in the background
177        answer += bkg
178
179        Return (answer)
180End             //End of function PolyLenExclVolCyl(w,x)
181
182Function fle_kernel(qw,rad,len_avg,lb,zz,delrho,len_i)
183        Variable qw,rad,len_avg,lb,zz,delrho,len_i
184       
185        //ww[0] = scale
186        //ww[1] = L [A]
187        //ww[2] = B [A]
188        //ww[3] = rad [A] cross-sectional radius
189        //ww[4] = contrast [A^-2]
190        //ww[5] = bkg [cm-1]
191        Variable Pq,vcyl,dl
192        Make/O/n=6 fle_ker
193        Wave kp = fle_ker
194        kp[0] = 1               //scale fixed at 1
195        kp[1] = len_i
196        kp[2] = lb
197        kp[3] = rad
198        kp[4] = delrho
199        kp[5] = 0               //bkg fixed at 0
200       
201#if exists("FlexExclVolCylX")
202        Pq = FlexExclVolCylX(kp,qw)
203#else
204        Pq = fFlexExclVolCyl(kp,qw)
205#endif
206       
207        vcyl=Pi*rad*rad*len_i
208        Pq *= vcyl
209        //un-convert from [cm-1]
210        Pq /= 1.0e8
211       
212        dl = Schulz_Point_flepl(len_i,len_avg,zz)
213        return (Pq*dl) 
214End
215
216//wrapper to calculate the smeared model as an AAO-Struct
217// fills the struct and calls the ususal function with the STRUCT parameter
218//
219// used only for the dependency, not for fitting
220//
221Function fSmearedFlexCyl_PolyLen(coefW,yW,xW)
222        Wave coefW,yW,xW
223       
224        String str = getWavesDataFolder(yW,0)
225        String DF="root:"+str+":"
226       
227        WAVE resW = $(DF+str+"_res")
228       
229        STRUCT ResSmearAAOStruct fs
230        WAVE fs.coefW = coefW   
231        WAVE fs.yW = yW
232        WAVE fs.xW = xW
233        WAVE fs.resW = resW
234       
235        Variable err
236        err = SmearedFlexCyl_PolyLen(fs)
237       
238        return (0)
239End
240
241// this is all there is to the smeared calculation!
242Function SmearedFlexCyl_PolyLen(s) :FitFunc
243        Struct ResSmearAAOStruct &s
244
245//      the name of your unsmeared model (AAO) is the first argument
246        s.yW = Smear_Model_20(FlexCyl_PolyLen,s.coefW,s.xW,s.resW)
247
248        return(0)
249End
Note: See TracBrowser for help on using the repository browser.