source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/FlexCyl_PolyRadius_v40.ipf @ 345

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

Merging NewModels_2006 back into where it belongs

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