source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/OblateCoreShell_v40.ipf @ 510

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

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

File size: 7.9 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","parameters_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","smear_parameters_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.