source: sans/Analysis/trunk/Put in User Procedures/SANS_Models_v3.00/OblateForm.ipf @ 56

Last change on this file since 56 was 42, checked in by srkline, 16 years ago

initial checkin of Analysis v.3.00 files

File size: 6.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3////////////////////////////////////////////////
4// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile
5// Adopting these into the experiment will insure that they are always present
6////////////////////////////////////////////////
7//
8// this function is for the form factor of an oblate ellipsoid with a core-shell structure
9//
10// 06 NOV 98 SRK
11////////////////////////////////////////////////
12
13Proc PlotOblateForm(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_oef,ywave_oef
20        xwave_oef =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
21        Make/O/D coef_oef = {1.,200,20,250,30,1e-6,1e-6,0.001}
22        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)"}
23        Edit parameters_oef,coef_oef
24        ywave_oef := OblateForm(coef_oef,xwave_oef)
25        Display ywave_oef vs xwave_oef
26        ModifyGraph log=1,marker=29,msize=2,mode=4
27        Label bottom "q (\\S-1\\M)"
28        Label left "Intensity (cm\\S-1\\M)"
29        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
30End
31///////////////////////////////////////////////////////////
32
33Proc PlotSmearedOblateForm()                                                            //**** name of your function
34        //no input parameters necessary, it MUST use the experimental q-values
35        // from the experimental data read in from an AVE/QSIG data file
36       
37        // if no gQvals wave, data must not have been loaded => abort
38        if(ResolutionWavesMissing())
39                Abort
40        endif
41       
42        // Setup parameter table for model function
43        Make/O/D smear_coef_oef = {1.,200,20,250,30,1e-6,1e-6,0.001}
44        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)"}
45        Edit smear_parameters_oef,smear_coef_oef
46       
47        // output smeared intensity wave, dimensions are identical to experimental QSIG values
48        // make extra copy of experimental q-values for easy plotting
49        Duplicate/O $gQvals smeared_oef,smeared_qvals                           
50        SetScale d,0,0,"1/cm",smeared_oef                                       
51
52        smeared_oef := SmearedOblateForm(smear_coef_oef,$gQvals)
53        Display smeared_oef vs smeared_qvals
54        ModifyGraph log=1,marker=29,msize=2,mode=4
55        Label bottom "q (\\S-1\\M)"
56        Label left "Intensity (cm\\S-1\\M)"
57        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
58End
59
60///////////////////////////////////////////////////////////////
61// unsmeared model calculation
62///////////////////////////
63Function OblateForm(w,x) : FitFunc
64        Wave w
65        Variable x
66
67//The input variables are (and output)
68        //[0] scale
69        //[1] crmaj, major radius of core       []
70        //[2] crmin, minor radius of core
71        //[3] trmaj, overall major radius
72        //[4] trmin, overall minor radius
73        //[5] delpc, SLD difference (core-shell) [-2]
74        //[6] delps, SLD difference (shell-solvent)
75        //[7] bkg, [cm-1]
76        Variable scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg
77        scale = w[0]
78        crmaj = w[1]
79        crmin = w[2]
80        trmaj = w[3]
81        trmin = w[4]
82        delpc = w[5]
83        delps = w[6]
84        bkg = w[7]
85
86// local variables
87        Variable yyy,va,vb,ii,nord,zi,qq,summ,nfn,npro,answer,oblatevol
88        String weightStr,zStr
89       
90        weightStr = "gauss76wt"
91        zStr = "gauss76z"
92
93       
94//      if wt,z waves don't exist, create them
95
96        if (WaveExists($weightStr) == 0) // wave reference is not valid,
97                Make/D/N=76 $weightStr,$zStr
98                Wave w76 = $weightStr
99                Wave z76 = $zStr                // wave references to pass
100                Make76GaussPoints(w76,z76)     
101        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
102        else
103                if(exists(weightStr) > 1)
104                         Abort "wave name is already in use"    // execute if condition is false
105                endif
106                Wave w76 = $weightStr
107                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
108        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
109        endif
110
111// set up the integration
112        // end points and weights
113        nord = 76
114        nfn = 2         //only <f^2> is calculated
115        npro = 0        // OBLATE ELLIPSOIDS
116        va =0
117        vb =1
118
119        qq = x          //current x point is the q-value for evaluation
120      summ = 0.0
121      ii=0
122      do
123                //printf "top of nord loop, i = %g\r",i
124        if(nfn ==1) //then              // "f1" required for beta factor
125          if(npro ==1) //then   // prolate
126                 zi = ( z76[ii]*(vb-va) + vb + va )/2.0
127//            yyy = w76[ii]*gfn1(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
128          Endif
129//
130          if(npro ==0) //then   // oblate 
131                 zi = ( z76[ii]*(vb-va) + vb + va )/2.0
132//            yyy = w76[ii]*gfn3(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
133          Endif
134        Endif           //nfn = 1
135        //
136        if(nfn !=1) //then              //calculate"f2" = <f^2> = averaged form factor
137          if(npro ==1) //then   //prolate
138             zi = ( z76[ii]*(vb-va) + vb + va )/2.0
139//            yyy = w76[ii]*gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
140          //printf "yyy = %g\r",yyy
141          Endif
142//
143          if(npro ==0) //then   //oblate
144                 zi = ( z76[ii]*(vb-va) + vb + va )/2.0
145                yyy = w76[ii]*gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
146          Endif
147        Endif           //nfn <>1
148       
149        summ = yyy + summ               // get running total of integral
150        ii+=1
151        while (ii<nord)                         // end of loop over quadrature points
152//   
153// calculate value of integral to return
154
155      answer = (vb-va)/2.0*summ
156     
157      // normalize by particle volume
158      oblatevol = 4*Pi/3*trmaj*trmaj*trmin
159      answer /= oblatevol
160     
161      //convert answer [-1] to [cm-1]
162      answer *= 1.0e8 
163      //scale
164      answer *= scale
165      // //then add background
166      answer += bkg
167
168        Return (answer)
169End
170//
171//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
172//                       BY (53) & (58-59) IN CHEN AND
173//                       KOTLARCHYK REFERENCE
174//
175//       <OBLATE ELLIPSOID>
176
177Function gfn4(xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
178        Variable xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq
179        // local variables
180        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43
181       
182        PI43=4.0/3.0*PI
183        aa = crmaj
184        bb = crmin
185        u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx))
186        ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx))
187        uq = sqrt(u2)*qq
188        ut= sqrt(ut2)*qq
189        vc = PI43*aa*aa*bb
190        vt = PI43*trmaj*trmaj*trmin
191        gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc
192        gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps
193        tgfn = gfnc+gfnt
194        gfn4 = tgfn*tgfn
195       
196        return gfn4
197       
198End             // function gfn4 for oblate ellipsoids
199
200// this is all there is to the smeared calculation!
201Function SmearedOblateForm(w,x) :FitFunc
202        Wave w
203        Variable x
204       
205        Variable ans
206        SVAR sq = gSig_Q
207        SVAR qb = gQ_bar
208        SVAR sh = gShadow
209        SVAR gQ = gQVals
210       
211        //the name of your unsmeared model is the first argument
212        ans = Smear_Model_20(OblateForm,$sq,$qb,$sh,$gQ,w,x)
213
214        return(ans)
215End
Note: See TracBrowser for help on using the repository browser.