source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/FlexCyl_PolyLen_v40.ipf @ 381

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

Merging Dev/trunk revision 374+ into Release/trunk for version 6.004

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","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","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.