source: sans/Dev/trunk/NCNR_User_Procedures/SANS_Analysis/Models/ProlateCoreShell_v40.ipf @ 325

Last change on this file since 325 was 325, checked in by ajj, 15 years ago

Adding SANS Analysis Models to new Dev tree

File size: 7.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////
5// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile
6// Adopting these into the experiment will insure that they are always present
7////////////////////////////////////////////////
8// this function is for the form factor of a prolate ellipsoid
9//
10// 06 NOV 98 SRK
11////////////////////////////////////////////////
12
13Proc PlotProlateForm(num,qmin,qmax)
14        Variable num=128,qmin=0.001,qmax=0.7
15        Prompt num "Enter number of data points for model: "
16        Prompt qmin "Enter minimum q-value (^-1) for model: "
17        Prompt qmax "Enter maximum q-value (^-1) for model: "
18       
19        Make/O/D/n=(num) xwave_pef,ywave_pef
20        xwave_pef =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
21        Make/O/D coef_pef = {1.,100,50,110,60,1e-6,2e-6,6.3e-6,0.001}
22        make/o/t parameters_pef = {"scale","major core radius (A)","minor core radius (A)","major shell radius (A)","minor shell radius (A)","SLD core (A-2)","SLD shell (A-2)","SLD solvent (A-2)","bkg (cm-1)"}
23        Edit parameters_pef,coef_pef
24        Variable/G root:g_pef
25        g_pef := ProlateForm(coef_pef,ywave_pef,xwave_pef)
26//      ywave_pef := ProlateForm(coef_pef,xwave_pef)
27        Display ywave_pef vs xwave_pef
28        ModifyGraph log=1,marker=29,msize=2,mode=4
29        Label bottom "q (\\S-1\\M)"
30        Label left "Intensity (cm\\S-1\\M)"
31        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
32       
33        AddModelToStrings("ProlateForm","coef_pef","pef")
34End
35
36///////////////////////////////////////////////////////////
37// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
38Proc PlotSmearedProlateForm(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_pef = {1.,100,50,110,60,1e-6,2e-6,6.3e-6,0.001}
51        make/o/t smear_parameters_pef = {"scale","major core radius (A)","minor core radius (A)","major shell radius (A)","minor shell radius (A)","SLD core (A-2)","SLD shell (A-2)","SLD solvent (A-2)","bkg (cm-1)"}
52        Edit smear_parameters_pef,smear_coef_pef
53       
54        // output smeared intensity wave, dimensions are identical to experimental QSIG values
55        // make extra copy of experimental q-values for easy plotting
56        Duplicate/O $(str+"_q") smeared_pef,smeared_qvals                               //**** mod
57        SetScale d,0,0,"1/cm",smeared_pef                                               
58               
59        Variable/G gs_pef=0
60        gs_pef := fSmearedProlateForm(smear_coef_pef,smeared_pef,smeared_qvals) //this wrapper fills the STRUCT
61
62        Display smeared_pef vs smeared_qvals
63        ModifyGraph log=1,marker=29,msize=2,mode=4
64        Label bottom "q (\\S-1\\M)"
65        Label left "Intensity (cm\\S-1\\M)"
66        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) 
67       
68        SetDataFolder root:
69        AddModelToStrings("SmearedProlateForm","smear_coef_pef","pef")
70End
71
72//AAO version
73Function ProlateForm(cw,yw,xw) : FitFunc
74        Wave cw,yw,xw
75
76#if exists("ProlateFormX")
77        yw = ProlateFormX(cw,xw)
78#else
79        yw = fProlateForm(cw,xw)
80#endif
81        return(0)
82End
83
84///////////////////////////////////////////////////////////////
85// unsmeared model calculation
86///////////////////////////
87Function fProlateForm(w,x) : FitFunc
88        Wave w
89        Variable x
90
91//The input variables are (and output)
92        //[0] scale
93        //[1] crmaj, major radius of core       []
94        //[2] crmin, minor radius of core
95        //[3] trmaj, overall major radius
96        //[4] trmin, overall minor radius
97        //[5] sld core, [A^-2]
98        //[6] sld shell,
99        //[7] sld solvent
100        //[8] bkg [cm-1]
101        Variable scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg,sldc,slds,sld
102        scale = w[0]
103        crmaj = w[1]
104        crmin = w[2]
105        trmaj = w[3]
106        trmin = w[4]
107        sldc = w[5]
108        slds = w[6]
109        sld = w[7]
110        bkg = w[8]
111
112        delpc = sldc - slds                     //core - shell
113        delps = slds - sld                      //shell - solvent
114// local variables
115        Variable yyy,va,vb,ii,nord,zi,qq,summ,nfn,npro,answer,prolatevol
116        String weightStr,zStr
117       
118        weightStr = "gauss76wt"
119        zStr = "gauss76z"
120
121//      if wt,z waves don't exist, create them
122
123        if (WaveExists($weightStr) == 0) // wave reference is not valid,
124                Make/D/N=76 $weightStr,$zStr
125                Wave w76 = $weightStr
126                Wave z76 = $zStr                // wave references to pass
127                Make76GaussPoints(w76,z76)     
128        else
129                if(exists(weightStr) > 1)
130                         Abort "wave name is already in use"    // execute if condition is false
131                endif
132                Wave w76 = $weightStr
133                Wave z76 = $zStr
134        endif
135
136// set up the integration
137        // end points and weights
138        nord = 76
139        nfn = 2         //only <f^2> is calculated
140        npro = 1        // PROLATE ELLIPSOIDS
141        va =0
142        vb =1
143//move this zi(i) evaluation inside other nord loop, since I don't have an array
144//      i=0
145//      do
146//       zi[i] = ( z76[i]*(vb-va) + vb + va )/2.0
147 //       i +=1
148 //     while (i<nord)
149//
150// evaluate at Gauss points
151        // remember to index from 0,size-1
152
153        qq = x          //current x point is the q-value for evaluation
154        summ = 0.0
155        ii=0
156        do
157                //printf "top of nord loop, i = %g\r",i
158                if(nfn ==1) //then              // "f1" required for beta factor
159                        if(npro ==1) //then     // prolate
160                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0 
161//           yyy = w76[ii]*gfn1(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
162                        Endif
163//
164                        if(npro ==0) //then     // oblate 
165                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
166//            yyy = w76[i]*gfn3(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
167                        Endif
168                Endif           //nfn = 1
169          //
170                if(nfn !=1) //then              //calculate"f2" = <f^2> = averaged form factor
171                        if(npro ==1) //then     //prolate
172                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
173                                yyy = w76[ii]*gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
174                                //printf "yyy = %g\r",yyy
175                        Endif
176//
177                        if(npro ==0) //then     //oblate
178                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
179//              yyy = w76[ii]*gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
180                        Endif
181                Endif           //nfn <>1
182         
183                summ = yyy + summ               // get running total of integral
184                ii+=1
185        while (ii<nord)                         // end of loop over quadrature points
186//   
187// calculate value of integral to return
188
189        answer = (vb-va)/2.0*summ
190       
191        //normailze by particle volume
192        prolatevol = 4*Pi/3*trmaj*trmin*trmin
193        answer /= prolatevol
194       
195        // rescale from 1/ to 1/cm
196        answer *= 1.0e8
197        //scale (arb)
198        answer *= scale
199        ////then add in background
200        answer += bkg
201
202        Return (answer)
203End     //prolate form factor
204
205//
206//     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
207//                       BY (53) AND (56,57) IN CHEN AND
208//                       KOTLARCHYK REFERENCE
209//
210//     <PROLATE ELLIPSOIDS>
211//
212Function gfn2(xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
213        Variable xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq
214        // local variables
215        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn2,pi43,gfn
216
217        PI43=4.0/3.0*PI
218        aa = crmaj
219        bb = crmin
220        u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx))
221        ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx))
222        uq = sqrt(u2)*qq
223        ut= sqrt(ut2)*qq
224        vc = PI43*aa*bb*bb
225        vt = PI43*trmaj*trmin*trmin
226        gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc
227        gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps
228        gfn = gfnc+gfnt
229        gfn2 = gfn*gfn
230       
231        return gfn2
232End             //function gfn2 for prolate ellipsoids
233
234
235// this is all there is to the smeared calculation!
236Function SmearedProlateForm(s) :FitFunc
237        Struct ResSmearAAOStruct &s
238
239////the name of your unsmeared model is the first argument
240        Smear_Model_20(ProlateForm,s.coefW,s.xW,s.yW,s.resW)
241
242        return(0)
243End
244
245//wrapper to calculate the smeared model as an AAO-Struct
246// fills the struct and calls the ususal function with the STRUCT parameter
247//
248// used only for the dependency, not for fitting
249//
250Function fSmearedProlateForm(coefW,yW,xW)
251        Wave coefW,yW,xW
252       
253        String str = getWavesDataFolder(yW,0)
254        String DF="root:"+str+":"
255       
256        WAVE resW = $(DF+str+"_res")
257       
258        STRUCT ResSmearAAOStruct fs
259        WAVE fs.coefW = coefW   
260        WAVE fs.yW = yW
261        WAVE fs.xW = xW
262        WAVE fs.resW = resW
263       
264        Variable err
265        err = SmearedProlateForm(fs)
266       
267        return (0)
268End
Note: See TracBrowser for help on using the repository browser.