source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/ProlateForm.ipf @ 146

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

Typo in dialog in every model...

File size: 7.5 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,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)","Contrast (core-shell) (A-2)","Constrast (shell-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)
32End
33
34///////////////////////////////////////////////////////////
35// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
36Proc PlotSmearedProlateForm(str)                                                               
37        String str
38        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
39       
40        // if any of the resolution waves are missing => abort
41        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
42                Abort
43        endif
44       
45        SetDataFolder $("root:"+str)
46       
47        // Setup parameter table for model function
48        Make/O/D smear_coef_pef = {1.,100,50,110,60,1e-6,2e-6,0.001}
49        make/o/t smear_parameters_pef = {"scale","major core radius (A)","minor core radius (A)","major shell radius (A)","minor shell radius (A)","Contrast (core-shell) (A-2)","Contrast (shell-solvent) (A-2)","bkg (cm-1)"}
50        Edit smear_parameters_pef,smear_coef_pef
51       
52        // output smeared intensity wave, dimensions are identical to experimental QSIG values
53        // make extra copy of experimental q-values for easy plotting
54        Duplicate/O $(str+"_q") smeared_pef,smeared_qvals                               //**** mod
55        SetScale d,0,0,"1/cm",smeared_pef                                               
56               
57        Variable/G gs_pef=0
58        gs_pef := fSmearedProlateForm(smear_coef_pef,smeared_pef,smeared_qvals) //this wrapper fills the STRUCT
59
60        Display smeared_pef vs smeared_qvals
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:
67End
68
69//AAO version
70Function ProlateForm(cw,yw,xw) : FitFunc
71        Wave cw,yw,xw
72
73#if exists("ProlateFormX")
74        yw = ProlateFormX(cw,xw)
75#else
76        yw = fProlateForm(cw,xw)
77#endif
78        return(0)
79End
80
81///////////////////////////////////////////////////////////////
82// unsmeared model calculation
83///////////////////////////
84Function fProlateForm(w,x) : FitFunc
85        Wave w
86        Variable x
87
88//The input variables are (and output)
89        //[0] scale
90        //[1] crmaj, major radius of core       []
91        //[2] crmin, minor radius of core
92        //[3] trmaj, overall major radius
93        //[4] trmin, overall minor radius
94        //[5] delpc, SLD difference (core-shell)        [-2]
95        //[6] delps, SLD difference (shell-solvent)
96        //[7] bkg [cm-1]
97        Variable scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg
98        scale = w[0]
99        crmaj = w[1]
100        crmin = w[2]
101        trmaj = w[3]
102        trmin = w[4]
103        delpc = w[5]
104        delps = w[6]
105        bkg = w[7]
106
107// local variables
108        Variable yyy,va,vb,ii,nord,zi,qq,summ,nfn,npro,answer,prolatevol
109        String weightStr,zStr
110       
111        weightStr = "gauss76wt"
112        zStr = "gauss76z"
113
114//      if wt,z waves don't exist, create them
115
116        if (WaveExists($weightStr) == 0) // wave reference is not valid,
117                Make/D/N=76 $weightStr,$zStr
118                Wave w76 = $weightStr
119                Wave z76 = $zStr                // wave references to pass
120                Make76GaussPoints(w76,z76)     
121        else
122                if(exists(weightStr) > 1)
123                         Abort "wave name is already in use"    // execute if condition is false
124                endif
125                Wave w76 = $weightStr
126                Wave z76 = $zStr
127        endif
128
129// set up the integration
130        // end points and weights
131        nord = 76
132        nfn = 2         //only <f^2> is calculated
133        npro = 1        // PROLATE ELLIPSOIDS
134        va =0
135        vb =1
136//move this zi(i) evaluation inside other nord loop, since I don't have an array
137//      i=0
138//      do
139//       zi[i] = ( z76[i]*(vb-va) + vb + va )/2.0
140 //       i +=1
141 //     while (i<nord)
142//
143// evaluate at Gauss points
144        // remember to index from 0,size-1
145
146        qq = x          //current x point is the q-value for evaluation
147        summ = 0.0
148        ii=0
149        do
150                //printf "top of nord loop, i = %g\r",i
151                if(nfn ==1) //then              // "f1" required for beta factor
152                        if(npro ==1) //then     // prolate
153                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0 
154//           yyy = w76[ii]*gfn1(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
155                        Endif
156//
157                        if(npro ==0) //then     // oblate 
158                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
159//            yyy = w76[i]*gfn3(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
160                        Endif
161                Endif           //nfn = 1
162          //
163                if(nfn !=1) //then              //calculate"f2" = <f^2> = averaged form factor
164                        if(npro ==1) //then     //prolate
165                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
166                                yyy = w76[ii]*gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
167                                //printf "yyy = %g\r",yyy
168                        Endif
169//
170                        if(npro ==0) //then     //oblate
171                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
172//              yyy = w76[ii]*gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
173                        Endif
174                Endif           //nfn <>1
175         
176                summ = yyy + summ               // get running total of integral
177                ii+=1
178        while (ii<nord)                         // end of loop over quadrature points
179//   
180// calculate value of integral to return
181
182        answer = (vb-va)/2.0*summ
183       
184        //normailze by particle volume
185        prolatevol = 4*Pi/3*trmaj*trmin*trmin
186        answer /= prolatevol
187       
188        // rescale from 1/ to 1/cm
189        answer *= 1.0e8
190        //scale (arb)
191        answer *= scale
192        ////then add in background
193        answer += bkg
194
195        Return (answer)
196End     //prolate form factor
197
198//
199//     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
200//                       BY (53) AND (56,57) IN CHEN AND
201//                       KOTLARCHYK REFERENCE
202//
203//     <PROLATE ELLIPSOIDS>
204//
205Function gfn2(xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
206        Variable xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq
207        // local variables
208        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn2,pi43,gfn
209
210        PI43=4.0/3.0*PI
211        aa = crmaj
212        bb = crmin
213        u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx))
214        ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx))
215        uq = sqrt(u2)*qq
216        ut= sqrt(ut2)*qq
217        vc = PI43*aa*bb*bb
218        vt = PI43*trmaj*trmin*trmin
219        gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc
220        gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps
221        gfn = gfnc+gfnt
222        gfn2 = gfn*gfn
223       
224        return gfn2
225End             //function gfn2 for prolate ellipsoids
226
227
228// this is all there is to the smeared calculation!
229Function SmearedProlateForm(s) :FitFunc
230        Struct ResSmearAAOStruct &s
231
232////the name of your unsmeared model is the first argument
233        s.yW = Smear_Model_20(ProlateForm,s.coefW,s.xW,s.resW)
234
235        return(0)
236End
237
238//wrapper to calculate the smeared model as an AAO-Struct
239// fills the struct and calls the ususal function with the STRUCT parameter
240//
241// used only for the dependency, not for fitting
242//
243Function fSmearedProlateForm(coefW,yW,xW)
244        Wave coefW,yW,xW
245       
246        String str = getWavesDataFolder(yW,0)
247        String DF="root:"+str+":"
248       
249        WAVE resW = $(DF+str+"_res")
250       
251        STRUCT ResSmearAAOStruct fs
252        WAVE fs.coefW = coefW   
253        WAVE fs.yW = yW
254        WAVE fs.xW = xW
255        WAVE fs.resW = resW
256       
257        Variable err
258        err = SmearedProlateForm(fs)
259       
260        return (0)
261End
Note: See TracBrowser for help on using the repository browser.