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

Last change on this file since 92 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// this function is for the form factor of a prolate ellipsoid
8//
9// 06 NOV 98 SRK
10////////////////////////////////////////////////
11
12Proc PlotProlateForm(num,qmin,qmax)
13        Variable num=128,qmin=0.001,qmax=0.7
14        Prompt num "Enter number of data points for model: "
15        Prompt qmin "Enter minimum q-value (^-1) for model: "
16        Prompt qmax "Enter maximum q-value (^-1) for model: "
17
18        Make/O/D/n=(num) xwave_pef,ywave_pef
19        xwave_pef =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
20        Make/O/D coef_pef = {1.,100,50,110,60,1e-6,2e-6,0.001}
22        Edit parameters_pef,coef_pef
23        ywave_pef := ProlateForm(coef_pef,xwave_pef)
24        Display ywave_pef vs xwave_pef
25        ModifyGraph log=1,marker=29,msize=2,mode=4
26        Label bottom "q (\\S-1\\M)"
27        Label left "Intensity (cm\\S-1\\M)"
28        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
29End
30///////////////////////////////////////////////////////////
31
32Proc PlotSmearedProlateForm()                                                           //**** name of your function
33        //no input parameters necessary, it MUST use the experimental q-values
34        // from the experimental data read in from an AVE/QSIG data file
35
36        // if no gQvals wave, data must not have been loaded => abort
37        if(ResolutionWavesMissing())
38                Abort
39        endif
40
41        // Setup parameter table for model function
42        Make/O/D smear_coef_pef = {1.,100,50,110,60,1e-6,2e-6,0.001}
44        Edit smear_parameters_pef,smear_coef_pef
45
46        // output smeared intensity wave, dimensions are identical to experimental QSIG values
47        // make extra copy of experimental q-values for easy plotting
48        Duplicate/O \$gQvals smeared_pef,smeared_qvals                           //**** mod
49        SetScale d,0,0,"1/cm",smeared_pef                                                       //**** mod
50
51        smeared_pef := SmearedProlateForm(smear_coef_pef,\$gQvals)
52        Display smeared_pef vs smeared_qvals
53        ModifyGraph log=1,marker=29,msize=2,mode=4
54        Label bottom "q (\\S-1\\M)"
55        Label left "Intensity (cm\\S-1\\M)"
56        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
57End
58
59///////////////////////////////////////////////////////////////
60// unsmeared model calculation
61///////////////////////////
62Function ProlateForm(w,x) : FitFunc
63        Wave w
64        Variable x
65
66//The input variables are (and output)
67        //[0] scale
68        //[1] crmaj, major radius of core       []
69        //[2] crmin, minor radius of core
70        //[3] trmaj, overall major radius
71        //[4] trmin, overall minor radius
72        //[5] delpc, SLD difference (core-shell)        [-2]
73        //[6] delps, SLD difference (shell-solvent)
74        //[7] bkg [cm-1]
75        Variable scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg
76        scale = w[0]
77        crmaj = w[1]
78        crmin = w[2]
79        trmaj = w[3]
80        trmin = w[4]
81        delpc = w[5]
82        delps = w[6]
83        bkg = w[7]
84
85// local variables
87        String weightStr,zStr
88
89        weightStr = "gauss76wt"
90        zStr = "gauss76z"
91
92//      if wt,z waves don't exist, create them
93
94        if (WaveExists(\$weightStr) == 0) // wave reference is not valid,
95                Make/D/N=76 \$weightStr,\$zStr
96                Wave w76 = \$weightStr
97                Wave z76 = \$zStr                // wave references to pass
98                Make76GaussPoints(w76,z76)
99        else
100                if(exists(weightStr) > 1)
101                         Abort "wave name is already in use"    // execute if condition is false
102                endif
103                Wave w76 = \$weightStr
104                Wave z76 = \$zStr
105        endif
106
107// set up the integration
108        // end points and weights
109        nord = 76
110        nfn = 2         //only <f^2> is calculated
111        npro = 1        // PROLATE ELLIPSOIDS
112        va =0
113        vb =1
114//move this zi(i) evaluation inside other nord loop, since I don't have an array
115//      i=0
116//      do
117//       zi[i] = ( z76[i]*(vb-va) + vb + va )/2.0
118 //       i +=1
119 //     while (i<nord)
120//
121// evaluate at Gauss points
122        // remember to index from 0,size-1
123
124        qq = x          //current x point is the q-value for evaluation
125        summ = 0.0
126        ii=0
127        do
128                //printf "top of nord loop, i = %g\r",i
129                if(nfn ==1) //then              // "f1" required for beta factor
130                        if(npro ==1) //then     // prolate
131                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
132//           yyy = w76[ii]*gfn1(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
133                        Endif
134//
135                        if(npro ==0) //then     // oblate
136                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
137//            yyy = w76[i]*gfn3(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
138                        Endif
139                Endif           //nfn = 1
140          //
141                if(nfn !=1) //then              //calculate"f2" = <f^2> = averaged form factor
142                        if(npro ==1) //then     //prolate
143                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
144                                yyy = w76[ii]*gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
145                                //printf "yyy = %g\r",yyy
146                        Endif
147//
148                        if(npro ==0) //then     //oblate
149                                zi = ( z76[ii]*(vb-va) + vb + va )/2.0
150//              yyy = w76[ii]*gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
151                        Endif
152                Endif           //nfn <>1
153
154                summ = yyy + summ               // get running total of integral
155                ii+=1
156        while (ii<nord)                         // end of loop over quadrature points
157//
158// calculate value of integral to return
159
161
162        //normailze by particle volume
163        prolatevol = 4*Pi/3*trmaj*trmin*trmin
165
166        // rescale from 1/ to 1/cm
168        //scale (arb)
172
174End     //prolate form factor
175
176//
177//     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
178//                       BY (53) AND (56,57) IN CHEN AND
179//                       KOTLARCHYK REFERENCE
180//
181//     <PROLATE ELLIPSOIDS>
182//
183Function gfn2(xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq)
184        Variable xx,crmaj,crmin,trmaj,trmin,delpc,delps,qq
185        // local variables
186        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn2,pi43,gfn
187
188        PI43=4.0/3.0*PI
189        aa = crmaj
190        bb = crmin
191        u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx))
192        ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx))
193        uq = sqrt(u2)*qq
194        ut= sqrt(ut2)*qq
195        vc = PI43*aa*bb*bb
196        vt = PI43*trmaj*trmin*trmin
197        gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc
198        gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps
199        gfn = gfnc+gfnt
200        gfn2 = gfn*gfn
201
202        return gfn2
203End             //function gfn2 for prolate ellipsoids
204
205
206// this is all there is to the smeared calculation!
207Function SmearedProlateForm(w,x) :FitFunc
208        Wave w
209        Variable x
210
211        Variable ans
212        SVAR sq = gSig_Q
213        SVAR qb = gQ_bar