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

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

Changed Plot* and PlotSmeared?* naming schemes to be all consistent prefixes for the actual function name, so that the macros can be constructed from the function name, or vice versa.

also some tweaks to the wrapper to make sure that plot and append really work

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_PolyLen(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 PlotSmearedFlexCyl_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        Smear_Model_20(FlexCyl_PolyLen,s.coefW,s.xW,s.yW,s.resW)
247
248        return(0)
249End
Note: See TracBrowser for help on using the repository browser.