source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/FlexCyl_PolyLen_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.0 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4#include "FlexibleCylinder_v40"
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 (A^-1) for model: "
14        Prompt qmax "Enter maximum q-value (A^-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,1e-6,6.3e-6,0.0001}
20        make/o/t parameters_flepl =  {"scale","Contour Length (A)","polydispersity of Contour Length","Kuhn Length, b (A)","Radius (A)","SLD cylinder (A^-2)","SLD solvent (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 (A\\S-1\\M)"
28        Label left "Intensity (cm\\S-1\\M)"
29        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
30       
31        AddModelToStrings("FlexCyl_PolyLen","coef_flepl","parameters_flepl","flepl")
32End
33
34// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
35Proc PlotSmearedFlexCyl_PolyLen(str)                                                           
36        String str
37        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
38       
39        // if any of the resolution waves are missing => abort
40        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
41                Abort
42        endif
43       
44        SetDataFolder $("root:"+str)
45       
46        // Setup parameter table for model function
47        Make/O/D smear_coef_flepl = {1.,1000,0.001,100,20,1e-6,6.3e-6,0.0001}           //CH#4
48        make/o/t smear_parameters_flepl = {"scale","Contour Length (A)","polydispersity of Contour Length","Kuhn Length, b (A)","Radius (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","bkgd (cm^-1)"}
49        Edit smear_parameters_flepl,smear_coef_flepl                                    //display parameters in a table
50       
51        // output smeared intensity wave, dimensions are identical to experimental QSIG values
52        // make extra copy of experimental q-values for easy plotting
53        Duplicate/O $(str+"_q") smeared_flepl,smeared_qvals                             //
54        SetScale d,0,0,"1/cm",smeared_flepl                                                     //
55                                       
56        Variable/G gs_flepl=0
57        gs_flepl := fSmearedFlexCyl_PolyLen(smear_coef_flepl,smeared_flepl,smeared_qvals)       //this wrapper fills the STRUCT
58       
59        Display smeared_flepl vs smeared_qvals                                                                  //
60        ModifyGraph log=1,marker=29,msize=2,mode=4
61        Label bottom "q (A\\S-1\\M)"
62        Label left "I(q) (cm\\S-1\\M)"
63        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
64       
65        SetDataFolder root:
66        AddModelToStrings("SmearedFlexCyl_PolyLen","smear_coef_flepl","smear_parameters_flepl","flepl")
67End
68       
69
70Function Schulz_Point_flepl(x,avg,zz)
71
72        //Wave w
73        Variable x,avg,zz
74        Variable dr
75       
76        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
77       
78        return (exp(dr))
79       
80End
81
82
83//AAO version, uses XOP if available
84// simply calls the original single point calculation with
85// a wave assignment (this will behave nicely if given point ranges)
86Function FlexCyl_PolyLen(cw,yw,xw) : FitFunc
87        Wave cw,yw,xw
88       
89#if exists("FlexCyl_PolyLenX")
90        yw = FlexCyl_PolyLenX(cw,xw)
91#else
92        yw = fFlexCyl_PolyLen(cw,xw)
93#endif
94        return(0)
95End
96
97Function fFlexCyl_PolyLen(w,x) : FitFunc
98        Wave w
99        Variable x
100     
101        Variable scale,radius,pd,delrho,bkg,zz,length,lb,sldc,slds
102        scale = w[0]
103        length = w[1]
104        pd = w[2]
105        lb = w[3]
106        radius = w[4]
107        sldc = w[5]
108        slds = w[6]
109        bkg = w[7]
110       
111        delrho = sldc-slds
112       
113        zz = (1/pd)^2-1
114//
115// the OUTPUT form factor is <f^2>/Vavg [cm-1]
116//
117// local variables
118        Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
119        Variable answer,zp1,zp2,zp3,vpoly
120        String weightStr,zStr
121       
122        nord = 20
123        weightStr = "gauss20wt"
124        zStr = "gauss20z"
125
126// use 20 Gauss points
127        if (WaveExists($weightStr) == 0) // wave reference is not valid,
128                Make/D/N=(nord) $weightStr,$zStr
129                Wave wtGau = $weightStr
130                Wave zGau = $zStr               // wave references to pass
131                Make20GaussPoints(wtGau,zGau)   
132        else
133                if(exists(weightStr) > 1)
134                         Abort "wave name is already in use"    // execute if condition is false
135                endif
136                Wave wtGau = $weightStr
137                Wave zGau = $zStr
138        endif
139
140// set up the integration
141// end points and weights
142// limits are technically 0-inf, but wisely choose non-zero region of distribution
143        Variable range=3.4              //multiples of the std. dev. fom the mean
144        a = length*(1-range*pd)
145        if (a<0)
146                a=0             //otherwise numerical error when pd >= 0.3, making a<0
147        endif
148        If(pd>0.3)
149                range = 3.4 + (pd-0.3)*18
150        Endif
151        b = length*(1+range*pd) // is this far enough past avg length?
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                zi = ( zGau[ii]*(vb-va) + vb + va )/2.0         
162                yyy = wtGau[ii] * fle_kernel(qq,radius,length,lb,zz,sldc,slds,zi)
163                summ = yyy + summ
164                ii+=1
165        while (ii<nord)                         // end of loop over quadrature points
166//   
167// calculate value of integral to return
168   answer = (vb-va)/2.0*summ
169     
170//  contrast^2 is included in integration rad_kernel
171//      answer *= delrho*delrho
172//normalize by polydisperse volume
173// now volume depends on polydisperse Length - so normalize by the FIRST moment
174// 1st moment = volume!
175        vpoly = Pi*(radius)^2*length
176//Divide by vol, since volume has been "un-normalized" out
177        answer /= vpoly
178//convert to [cm-1]
179        answer *= 1.0e8
180//scale
181        answer *= scale
182// add in the background
183        answer += bkg
184
185        Return (answer)
186End             //End of function PolyLenExclVolCyl(w,x)
187
188Function fle_kernel(qw,rad,len_avg,lb,zz,sldc,slds,len_i)
189        Variable qw,rad,len_avg,lb,zz,sldc,slds,len_i
190       
191        //ww[0] = scale
192        //ww[1] = L [A]
193        //ww[2] = B [A]
194        //ww[3] = rad [A] cross-sectional radius
195        //ww[4] = sld cyl [A^-2]
196        //ww[5] = sld solv
197        //ww[6] = bkg [cm-1]
198        Variable Pq,vcyl,dl
199        Make/O/n=7 fle_ker
200        Wave kp = fle_ker
201        kp[0] = 1               //scale fixed at 1
202        kp[1] = len_i
203        kp[2] = lb
204        kp[3] = rad
205        kp[4] = sldc
206        kp[5] = slds
207        kp[6] = 0               //bkg fixed at 0
208       
209#if exists("FlexExclVolCylX")
210        Pq = FlexExclVolCylX(kp,qw)
211#else
212        Pq = fFlexExclVolCyl(kp,qw)
213#endif
214       
215        vcyl=Pi*rad*rad*len_i
216        Pq *= vcyl
217        //un-convert from [cm-1]
218        Pq /= 1.0e8
219       
220        dl = Schulz_Point_flepl(len_i,len_avg,zz)
221        return (Pq*dl) 
222End
223
224//wrapper to calculate the smeared model as an AAO-Struct
225// fills the struct and calls the ususal function with the STRUCT parameter
226//
227// used only for the dependency, not for fitting
228//
229Function fSmearedFlexCyl_PolyLen(coefW,yW,xW)
230        Wave coefW,yW,xW
231       
232        String str = getWavesDataFolder(yW,0)
233        String DF="root:"+str+":"
234       
235        WAVE resW = $(DF+str+"_res")
236       
237        STRUCT ResSmearAAOStruct fs
238        WAVE fs.coefW = coefW   
239        WAVE fs.yW = yW
240        WAVE fs.xW = xW
241        WAVE fs.resW = resW
242       
243        Variable err
244        err = SmearedFlexCyl_PolyLen(fs)
245       
246        return (0)
247End
248
249// this is all there is to the smeared calculation!
250Function SmearedFlexCyl_PolyLen(s) :FitFunc
251        Struct ResSmearAAOStruct &s
252
253//      the name of your unsmeared model (AAO) is the first argument
254        Smear_Model_20(FlexCyl_PolyLen,s.coefW,s.xW,s.yW,s.resW)
255
256        return(0)
257End
Note: See TracBrowser for help on using the repository browser.