source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/FlexCyl_PolyRadius.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: 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_PolyRadius(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        s.yW = Smear_Model_20(FlexCyl_PolyRad,s.coefW,s.xW,s.resW)
265
266        return(0)
267End
Note: See TracBrowser for help on using the repository browser.