source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/FlexCyl_PolyRadius_v40.ipf @ 510

Last change on this file since 510 was 510, checked in by srkline, 14 years ago

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

File size: 7.7 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 (A^-1) for model: "
13        Prompt qmax "Enter maximum q-value (A^-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 (A\\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","parameters_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 (A\\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","smear_parameters_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.