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