source: sans/Dev/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/EllipticalCylinder_v40.ipf @ 379

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

Removed Angstrom symbol (Mac code) and replaced it with simply a capital "A" so strange symbols won't appear anymore on Win.

File size: 5.6 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","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","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.