source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/EllipticalCylinder_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: 5.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4///////////////////////////////////////////
5// this function is for the form factor of an cylinder with an ellipsoidal cross-section
6// and a uniform scattering length density
7//
8// 06 NOV 98 SRK
9//
10// re-written to not use MacOS XOP for calculation
11// now requires the "new" version of GaussUtils that includes the generic quadrature routines
12//
13// 09 SEP 03 SRK
14////////////////////////////////////////////////
15
16Proc PlotEllipticalCylinder(num,qmin,qmax)
17        Variable num=50,qmin=0.001,qmax=0.7
18        Prompt num "Enter number of data points for model: "
19        Prompt qmin "Enter minimum q-value (A^-1) for model: "
20        Prompt qmax "Enter maximum q-value (A^-1) for model: "
21       
22//      //constants needed for the integration if qtrap is used (in a separate procedure file!)
23//      Variable/G root:gNumPoints=200
24//      Variable/G root:gTol=1e-5
25//      Variable/G root:gMaxIter=20
26        //
27        Make/O/D/n=(num) xwave_ecf,ywave_ecf
28        xwave_ecf =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
29        Make/O/D coef_ecf = {1.,20.,1.5,400,1.0e-6,6.3e-6,0.0}
30        make/o/t parameters_ecf = {"scale","minor radius (A)","nu = major/minor (-)","length (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
31        Edit parameters_ecf,coef_ecf
32       
33        Variable/G root:g_ecf
34        g_ecf := EllipticalCylinder(coef_ecf,ywave_ecf,xwave_ecf)
35        Display ywave_ecf vs xwave_ecf
36        ModifyGraph log=1,marker=29,msize=2,mode=4
37        Label bottom "q (A\\S-1\\M)"
38        Label left "Intensity (cm\\S-1\\M)"
39        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
40       
41        AddModelToStrings("EllipticalCylinder","coef_ecf","parameters_ecf","ecf")
42End
43///////////////////////////////////////////////////////////
44
45// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
46Proc PlotSmearedEllipticalCylinder(str)                                                         
47        String str
48        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
49       
50        // if any of the resolution waves are missing => abort
51        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
52                Abort
53        endif
54       
55        SetDataFolder $("root:"+str)
56       
57        // Setup parameter table for model function
58        Make/O/D smear_coef_ecf =  {1.,20.,1.5,400,1.0e-6,6.3e-6,0.0}
59        make/o/t smear_parameters_ecf = {"scale","minor radius (A)","nu = major/minor (-)","length (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
60        Edit smear_parameters_ecf,smear_coef_ecf
61       
62        // output smeared intensity wave, dimensions are identical to experimental QSIG values
63        // make extra copy of experimental q-values for easy plotting
64        Duplicate/O $(str+"_q") smeared_ecf,smeared_qvals
65        SetScale d,0,0,"1/cm",smeared_ecf       
66                                       
67        Variable/G gs_ecf=0
68        gs_ecf := fSmearedEllipticalCylinder(smear_coef_ecf,smeared_ecf,smeared_qvals)  //this wrapper fills the STRUCT
69       
70        Display smeared_ecf vs smeared_qvals
71        ModifyGraph log=1,marker=29,msize=2,mode=4
72        Label bottom "q (A\\S-1\\M)"
73        Label left "Intensity (cm\\S-1\\M)"
74        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
75       
76        SetDataFolder root:
77        AddModelToStrings("SmearedEllipticalCylinder","smear_coef_ecf","smear_parameters_ecf","ecf")
78End
79
80
81
82//AAO version, uses XOP if available
83// simply calls the original single point calculation with
84// a wave assignment (this will behave nicely if given point ranges)
85Function EllipticalCylinder(cw,yw,xw) : FitFunc
86        Wave cw,yw,xw
87       
88#if exists("EllipticalCylinderX")
89        yw = EllipticalCylinderX(cw,xw)
90#else
91        yw = fEllipticalCylinder(cw,xw)
92#endif
93        return(0)
94End
95
96// the main function that calculates the form factor of an elliptical cylinder
97// integrates EllipCyl_integrand, which itself is an integral function
98// 20 points of quadrature seems to be sufficient for both integrals
99//
100Function fEllipticalCylinder(w,x)       : FitFunc
101        Wave w
102        Variable x
103       
104        Variable inten,scale,rad,nu,len,contr,bkg,ii,sldc,slds
105        scale = w[0]
106        rad = w[1]
107        nu = w[2]
108        len = w[3]
109        sldc = w[4]
110        slds = w[5]
111        contr = sldc - slds
112        bkg = w[6]
113       
114        inten = IntegrateFn20(EllipCyl_Integrand,0,1,w,x)
115       
116        //multiply by volume
117        inten *= Pi*rad*rad*nu*len
118        inten *= 1e8    //convert to 1/cm
119        inten *= contr*contr
120        inten *= scale
121        inten += bkg
122       
123        return(inten)
124End
125
126//the outer integral
127Function EllipCyl_Integrand(w,x,dum)
128        Wave w
129        Variable x,dum
130       
131        Variable val,rad,arg,len
132        rad = w[1]
133        len = w[3]
134       
135        arg = rad*sqrt(1-dum^2)
136        duplicate/O w temp_w
137        Wave temp_w=temp_w
138        temp_w[1] = arg         //replace radius with transformed variable
139        val = (1/pi)*IntegrateFn20(Phi_EC,0,Pi,temp_w,x)
140       
141        // equivalent to the 20-pt quadrature
142//      val = (1/pi)*qtrap(Phi_EC,temp_w,x,0,Pi,1e-3,20)
143       
144        arg = x*len*dum/2
145        if(arg==0)
146                val *= 1
147        else
148                val *= sin(arg)*sin(arg)/arg/arg
149        endif
150        //Print "val=",val
151        return(val)
152End
153
154//the inner integral
155Function Phi_EC(w,x,dum)
156        Wave w
157        Variable x,dum
158       
159        Variable ans,arg,aa,nu
160        aa = w[1]               // = rad*sqrt(1-dum^2)
161        nu = w[2]
162        arg = x*aa*sqrt( (1+nu^2)/2 + (1-nu^2)/2*cos(dum) )
163        if(arg==0)
164                ans = (2*0.5)^2         // == 1
165        else
166                ans = 2*2*bessJ(1,arg)*bessJ(1,arg)/arg/arg
167        endif
168        return(ans)
169End
170
171//wrapper to calculate the smeared model as an AAO-Struct
172// fills the struct and calls the ususal function with the STRUCT parameter
173//
174// used only for the dependency, not for fitting
175//
176Function fSmearedEllipticalCylinder(coefW,yW,xW)
177        Wave coefW,yW,xW
178       
179        String str = getWavesDataFolder(yW,0)
180        String DF="root:"+str+":"
181       
182        WAVE resW = $(DF+str+"_res")
183       
184        STRUCT ResSmearAAOStruct fs
185        WAVE fs.coefW = coefW   
186        WAVE fs.yW = yW
187        WAVE fs.xW = xW
188        WAVE fs.resW = resW
189       
190        Variable err
191        err = SmearedEllipticalCylinder(fs)
192       
193        return (0)
194End
195
196// this is all there is to the smeared calculation!
197Function SmearedEllipticalCylinder(s) :FitFunc
198        Struct ResSmearAAOStruct &s
199
200//      the name of your unsmeared model (AAO) is the first argument
201        Smear_Model_20(EllipticalCylinder,s.coefW,s.xW,s.yW,s.resW)
202
203        return(0)
204End
205       
Note: See TracBrowser for help on using the repository browser.