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