source:sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models/Sphere_v40.ipf@338

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

merging structural changes from Dev/trunk into Release/trunk

File size: 6.7 KB
Line
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////
5// this procedure calculates the form factor of a sphere
6//
7// 06 NOV 98 SRK
8//
9// modified June 2007 to calculate all-at-once
10// and to use a structure for resolution information
11//
12////////////////////////////////////////////////
13
14Proc PlotSphereForm(num,qmin,qmax)
15        Variable num=128,qmin=0.001,qmax=0.7
16        Prompt num "Enter number of data points for model: "
17        Prompt qmin "Enter minimum q-value (^-1) for model: "
18        Prompt qmax "Enter maximum q-value (^-1) for model: "
19
20        Make/O/D/n=(num) xwave_sf,ywave_sf
21        xwave_sf = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
22        Make/O/D coef_sf = {1.,60,1e-6,6.3e-6,0.01}
23        make/o/t parameters_sf = {"scale","Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)","bkgd (cm-1)"}
24        Edit parameters_sf,coef_sf
25        Variable/G root:g_sf=0
26        root:g_sf := SphereForm(coef_sf,ywave_sf,xwave_sf)      // AAO calculation, "fake" dependency
27//      ywave_sf := SphereForm(coef_sf,xwave_sf)                //point calculation
28        Display ywave_sf vs xwave_sf
29        ModifyGraph log=1,marker=29,msize=2,mode=4
30        Label bottom "q (\\S-1\\M)"
31        Label left "Intensity (cm\\S-1\\M)"
32        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
33
34        AddModelToStrings("SphereForm","coef_sf","sf")
35End
36
37// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
38Proc PlotSmearedSphereForm(str)
39        String str
40        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
41
42        // if any of the resolution waves are missing => abort
43        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
44                Abort
45        endif
46
47        SetDataFolder \$("root:"+str)
48
49        // Setup parameter table for model function
50        Make/O/D smear_coef_sf = {1.,60,1e-6,6.3e-6,0.01}
51        make/o/t smear_parameters_sf = {"scale","Radius (A)","SLD sphere (A-2)","SLD solvent (A-2)","bkgd (cm-1)"}
52        Edit smear_parameters_sf,smear_coef_sf
53
54        Duplicate/O \$(str+"_q") smeared_sf,smeared_qvals
55        SetScale d,0,0,"1/cm",smeared_sf
56
57        Variable/G gs_sf=0
58        gs_sf := fSmearedSphereForm(smear_coef_sf,smeared_sf,smeared_qvals)     //this wrapper fills the STRUCT
59        Display smeared_sf vs smeared_qvals
60
61        ModifyGraph log=1,marker=29,msize=2,mode=4
62        Label bottom "q (\\S-1\\M)"
63        Label left "Intensity (cm\\S-1\\M)"
64        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
65
66        SetDataFolder root:
67        AddModelToStrings("SmearedSphereForm","smear_coef_sf","sf")
68End
69
70//AAO version, uses XOP if available
71// simply calls the original single point calculation with
72// a wave assignment (this will behave nicely if given point ranges)
73Function SphereForm(cw,yw,xw) : FitFunc
74        Wave cw,yw,xw
75
76#if exists("SphereFormX")
77        yw = SphereFormX(cw,xw)
78#else
79        yw = fSphereForm(cw,xw)
80#endif
81        return(0)
82End
83
84///////////////////////////////////////////////////////////////
85// unsmeared model calculation
86// this is a "regular" calculation at a single q-value
87///////////////////////////
88Function fSphereForm(w,x)
89        Wave w
90        Variable x
91
92        // variables are:
93        //[0] scale
94        //[1] radius ()
95        //[2] sld sphere (-2)
96        //[3] sld solv
97        //[4] background (cm-1)
98
99        Variable scale,radius,delrho,bkg,sldSph,sldSolv
100        scale = w[0]
101        radius = w[1]
102        sldSph = w[2]
103        sldSolv = w[3]
104        bkg = w[4]
105
106        delrho = sldSph - sldSolv
107        // calculates scale * f^2/Vol where f=Vol*3*delrho*((sin(qr)-qrcos(qr))/qr^3
108        // and is rescaled to give [=] cm^-1
109
110        Variable bes,f,vol,f2
111        //
112        //handle q==0 separately
113        If(x==0)
114                f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
115                return(f)
116        Endif
117
118        bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
119        vol = 4*pi/3*radius^3
120        f = vol*bes*delrho              // [=] 
121        // normalize to single particle volume, convert to 1/cm
122        f2 = f * f / vol * 1.0e8                // [=] 1/cm
123
124        return (scale*f2+bkg)   // Scale, then add in the background
125End
126
127
128//wrapper to calculate the smeared model as an AAO-Struct
129// fills the struct
130// calls the ususal function with the STRUCT parameter
131//
132// -- problem- this assumes USANS w/no matrix
133// -- needs experimental q-values to properly calculate the fitted curve?? is this true?
134// -- how to fill the struct if the type of data is unknown?
135// -- where do I find the matrix?
136//
137Function fSmearedSphereForm(coefW,yW,xW)
138        Wave coefW,yW,xW
139
140        String str = getWavesDataFolder(yW,0)
141        String DF="root:"+str+":"
142
143        WAVE resW = \$(DF+str+"_res")
144
145        STRUCT ResSmearAAOStruct fs
146        WAVE fs.coefW = coefW           //is this the proper way to populate? seems redundant...
147        WAVE fs.yW = yW
148        WAVE fs.xW = xW
149        WAVE fs.resW = resW
150
151        Variable err
152        err = SmearedSphereForm(fs)
153
154        return (0)
155End
156
157// smeared calculation, AAO and using a structure...
158// defined as as STRUCT, there can never be a dependency linked directly to this function
159// - so set a dependency to the wrapper
160//
161// like the unsmeared function, AAO is equivalent to a wave assignment to the point calculation
162// - but now the function passed is an AAO function
163//
164// Smear_Model_20() takes care of what calculation is done, depending on the resolution information
165//
166//
167Function SmearedSphereForm(s) : FitFunc
168        Struct ResSmearAAOStruct &s
169
170////the name of your unsmeared model is the first argument
171        Smear_Model_20(SphereForm,s.coefW,s.xW,s.yW,s.resW)
172
173        return(0)
174End
175
176
177// wrapper to do the desired fit
178//
179// str is the data folder for the desired data set
180//
181// -- this looks like something that can be made rather generic rather easily
182//
183//Function SphereFitWrapper(str)
184//      String str
185//
186//      SetDataFolder root:
187//      String DF="root:"+str+":"
188//
189//      Struct ResSmearAAOStruct fs
190//      WAVE resW = \$(DF+str+"_res")
191//      WAVE fs.resW =  resW
192//
193//      WAVE cw=\$(DF+"smear_coef_sf")
194//      WAVE yw=\$(DF+str+"_i")
195//      WAVE xw=\$(DF+str+"_q")
196//      WAVE sw=\$(DF+str+"_s")
197//
198//      Duplicate/O yw \$(DF+"FitYw")
199//
200//      //can't use execute if /STRC is needed since structure instance is not a global!
201//      //don't use the auto-destination with no flag, it doesn't appear to work correctly
202//      //force it to use a wave of identical length, at least - no guarantee that the q-values
203//      // will be the same? be sure that the smearing iterpolates
204//      //
205//      // ?? how can I get the hold string in correctly?? - from a global? - no, as string function
206//      //
207//      FuncFit/H=getHoldStr() /NTHR=0 SmearedSphereForm, cw, yw /X=xw /W=sw /I=1 /STRC=fs
208////    FuncFit/H=getHoldStr() /NTHR=0 SphereForm cw, yw /X=xw /W=sw /I=1 /D=\$(DF+"FitYw")
209//
210////    FuncFit/H="0010"/NTHR=0 SmearedSphereForm cw, yw /X=xw /W=sw /I=1 /D=\$(DF+"FitYw") /STRC=fs
211//      Wave fityw =  \$(DF+"FitYw")
212//      fs.yW = fityw
213//      SmearedSphereForm(fs)
214//      AppendToGraph fityw vs xw
215//
216//      print "V_chisq = ",V_chisq
217//      print cw
218//      WAVE w_sigma
219//      print w_sigma
220//
221//      return(0)
222//End
223//
224//Function/S getHoldStr()
225//
226//      String str="0010"
227//      return str
228//End
Note: See TracBrowser for help on using the repository browser.