source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models/OblateCoreShell_v40.ipf @ 381

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

Merging Dev/trunk revision 374+ into Release/trunk for version 6.004

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