source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/Cylinder_PolyRadius.ipf @ 211

Last change on this file since 211 was 211, checked in by srkline, 15 years ago

Several disjoint changes:

1) Threaded the 1D Cylinder_PolyRadius function
2) in Wrapper, popping the data set list now will auto-pop the function, then the coef to properly update the table of coefficients.
3) tick units in the plots are suppressed (y-axis)
4) Killing data (PlotManager?) now kills the dependent variable so that the rest can be killed
5) Ellipsoid_2D function now has the correct notation of the rotation axis in the parameter names.

File size: 9.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4#include "CylinderForm"
5
6// calculates the form factor of a cylinder with polydispersity of radius
7// the length distribution is a Schulz distribution, and any normalized distribution
8// could be used, as the average is performed numerically
9//
10// since the cylinder form factor is already a numerical integration, the size average is a
11// second integral, and significantly slows the calculation, and smearing adds a third integration.
12//
13//CORRECT! 12/5/2000 - Invariant is now correct vs. monodisperse cylinders
14// + upper limit of integration has been changed to account for skew of
15//Schulz distribution at high (>0.5) polydispersity
16//Requires 20 gauss points for integration of the radius (5 is not enough)
17//Requires either CylinderFit XOP (MacOSX only) or the normal CylinderForm Function
18//
19Proc PlotCyl_PolyRadius(num,qmin,qmax)
20        Variable num=128,qmin=0.001,qmax=0.7
21        Prompt num "Enter number of data points for model: "
22        Prompt qmin "Enter minimum q-value (^-1) for model: "
23        Prompt qmax "Enter maximum q-value (^-1) for model: "
24       
25        make/o/d/n=(num) xwave_cypr,ywave_cypr
26        xwave_cypr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
27        make/o/d coef_cypr = {1.,20.,400,0.2,3.0e-6,0.01}
28        make/o/t parameters_cypr = {"scale","radius (A)","length (A)","polydispersity of Radius","SLD diff (A^-2)","incoh. bkg (cm^-1)"}
29        Edit parameters_cypr,coef_cypr
30       
31        Variable/G root:g_cypr
32        g_cypr := Cyl_PolyRadius(coef_cypr,ywave_cypr,xwave_cypr)
33        Display ywave_cypr vs xwave_cypr
34        ModifyGraph log=1,marker=29,msize=2,mode=4
35        Label bottom "q (\\S-1\\M)"
36        Label left "Intensity (cm\\S-1\\M)"
37        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
38       
39        AddModelToStrings("Cyl_PolyRadius","coef_cypr","cypr")
40End
41
42///////////////////////////////////////////////////////////
43// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
44Proc PlotSmearedCyl_PolyRadius(str)                                                             
45        String str
46        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
47       
48        // if any of the resolution waves are missing => abort
49        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
50                Abort
51        endif
52       
53        SetDataFolder $("root:"+str)
54       
55        // Setup parameter table for model function
56        make/o/D smear_coef_cypr = {1.,20.,400,0.2,3.0e-6,0.01}
57        make/o/t smear_parameters_cypr = {"scale","radius (A)","length (A)","polydispersity of Radius","SLD diff (A^-2)","incoh. bkg (cm^-1)"}
58        Edit smear_parameters_cypr,smear_coef_cypr
59       
60        // output smeared intensity wave, dimensions are identical to experimental QSIG values
61        // make extra copy of experimental q-values for easy plotting
62        Duplicate/O $(str+"_q") smeared_cypr,smeared_qvals
63        SetScale d,0,0,"1/cm",smeared_cypr     
64                                       
65        Variable/G gs_cypr=0
66        gs_cypr := fSmearedCyl_PolyRadius(smear_coef_cypr,smeared_cypr,smeared_qvals)   //this wrapper fills the STRUCT
67       
68        Display smeared_cypr vs smeared_qvals
69        ModifyGraph log=1,marker=29,msize=2,mode=4
70        Label bottom "q (\\S-1\\M)"
71        Label left "Intensity (cm\\S-1\\M)"
72        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
73       
74        SetDataFolder root:
75        AddModelToStrings("SmearedCyl_PolyRadius","smear_coef_cypr","cypr")
76End
77       
78
79// non-threaded version, use the threaded version instead...
80//
81//AAO version, uses XOP if available
82// simply calls the original single point calculation with
83// a wave assignment (this will behave nicely if given point ranges)
84Function xCyl_PolyRadius(cw,yw,xw) : FitFunc
85        Wave cw,yw,xw
86       
87#if exists("Cyl_PolyRadiusX")
88        yw = Cyl_PolyRadiusX(cw,xw)
89#else
90        yw = fCyl_PolyRadius(cw,xw)
91#endif
92        return(0)
93End
94
95Function fCyl_PolyRadius(w,x) :FitFunc
96        Wave w
97        Variable x
98
99        //The input variables are (and output)
100        //[0] scale
101        //[1] avg RADIUS (A)
102        //[2] Length (A)
103        //[3] polydispersity (0<p<1)
104        //[4] contrast (A^-2)
105        //[5] background (cm^-1)
106        Variable scale,radius,pd,delrho,bkg,zz,length
107        scale = w[0]
108        radius = w[1]
109        length = w[2]
110        pd = w[3]
111        delrho = w[4]
112        bkg = w[5]
113       
114        zz = (1/pd)^2-1
115//
116// the OUTPUT form factor is <f^2>/Vavg [cm-1]
117//
118// local variables
119        Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
120        Variable answer,zp1,zp2,zp3,vpoly
121        String weightStr,zStr
122       
123//      nord = 5       
124//      weightStr = "gauss5wt"
125//      zStr = "gauss5z"
126        nord = 20
127        weightStr = "gauss20wt"
128        zStr = "gauss20z"
129
130//      if wt,z waves don't exist, create them
131// 5 Gauss points (not enough for cylinder radius = high q oscillations)
132// use 20 Gauss points
133        if (WaveExists($weightStr) == 0) // wave reference is not valid,
134                Make/D/N=(nord) $weightStr,$zStr
135                Wave wtGau = $weightStr
136                Wave zGau = $zStr               // wave references to pass
137                Make20GaussPoints(wtGau,zGau)   
138                //Make5GaussPoints(wtGau,zGau) 
139//      //                  printf "w[0],z[0] = %g %g\r", wtGau[0],zGau[0]
140        else
141                if(exists(weightStr) > 1)
142                         Abort "wave name is already in use"    // execute if condition is false
143                endif
144                Wave wtGau = $weightStr
145                Wave zGau = $zStr
146//      //          printf "w[0],z[0] = %g %g\r", wtGau[0],zGau[0]     
147        endif
148
149// set up the integration
150// end points and weights
151// limits are technically 0-inf, but wisely choose non-zero region of distribution
152        Variable range=3.4              //multiples of the std. dev. fom the mean
153        a = radius*(1-range*pd)
154        if (a<0)
155                a=0             //otherwise numerical error when pd >= 0.3, making a<0
156        endif
157        If(pd>0.3)
158                range = 3.4 + (pd-0.3)*18
159        Endif
160        b = radius*(1+range*pd) // is this far enough past avg radius?
161//      printf "a,b,ravg = %g %g %g\r", a,b,radius
162        va =a
163        vb =b
164
165// evaluate at Gauss points
166        // remember to index from 0,size-1     
167        qq = x          //current x point is the q-value for evaluation
168        summ = 0.0              // initialize integral
169   ii=0
170   do
171   //printf "top of nord loop, i = %g\r",i
172        // Using 5 Gauss points         
173                zi = ( zGau[ii]*(vb-va) + vb + va )/2.0         
174                yyy = wtGau[ii] * rad_kernel(qq,radius,length,zz,delrho,zi)
175                summ = yyy + summ
176                ii+=1
177        while (ii<nord)                         // end of loop over quadrature points
178//   
179// calculate value of integral to return
180   answer = (vb-va)/2.0*summ
181     
182//  contrast^2 is included in integration rad_kernel
183//      answer *= delrho*delrho
184//normalize by polydisperse volume
185// now volume depends on polydisperse RADIUS - so normalize by the second moment
186// 2nd moment = (zz+2)/(zz+1)
187        vpoly = Pi*(radius)^2*length*(zz+2)/(zz+1)
188//Divide by vol, since volume has been "un-normalized" out
189        answer /= vpoly
190//convert to [cm-1]
191        answer *= 1.0e8
192//scale
193        answer *= scale
194// add in the background
195        answer += bkg
196
197        Return (answer)
198End             //End of function PolyRadCylForm()
199
200Function rad_kernel(qw,ravg,len,zz,delrho,rad)
201        Variable qw,ravg,len,zz,delrho,rad
202       
203        Variable Pq,vcyl,dr
204       
205        //calculate the orientationally averaged P(q) for the input rad
206        //this is correct - see K&C (1983) or Lin &Tsao JACryst (1996)29 170.
207        Make/O/D/n=5 kernpar
208        Wave kp = kernpar
209        kp[0] = 1               //scale fixed at 1
210        kp[1] = rad
211        kp[2] = len
212        kp[3] = delrho
213        kp[4] = 0               //bkg fixed at 0
214       
215#if exists("CylinderFormX")
216        Pq = CylinderFormX(kp,qw)
217#else
218        Pq = fCylinderForm(kp,qw)
219#endif
220       
221        // undo the normalization that CylinderForm does
222        vcyl=Pi*rad*rad*len
223        Pq *= vcyl
224        //un-convert from [cm-1]
225        Pq /= 1.0e8
226       
227        // calculate normalized distribution at len value
228        dr = Schulz_Point_pr(rad,ravg,zz)
229       
230        return (Pq*dr) 
231End
232
233Function Schulz_Point_pr(x,avg,zz)
234        Variable x,avg,zz
235       
236        Variable dr
237       
238        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
239       
240        return (exp(dr))
241End
242
243//wrapper to calculate the smeared model as an AAO-Struct
244// fills the struct and calls the ususal function with the STRUCT parameter
245//
246// used only for the dependency, not for fitting
247//
248Function fSmearedCyl_PolyRadius(coefW,yW,xW)
249        Wave coefW,yW,xW
250       
251        String str = getWavesDataFolder(yW,0)
252        String DF="root:"+str+":"
253       
254        WAVE resW = $(DF+str+"_res")
255       
256        STRUCT ResSmearAAOStruct fs
257        WAVE fs.coefW = coefW   
258        WAVE fs.yW = yW
259        WAVE fs.xW = xW
260        WAVE fs.resW = resW
261       
262        Variable err
263        err = SmearedCyl_PolyRadius(fs)
264       
265        return (0)
266End
267
268// this is all there is to the smeared calculation!
269Function SmearedCyl_PolyRadius(s) :FitFunc
270        Struct ResSmearAAOStruct &s
271
272//      the name of your unsmeared model (AAO) is the first argument
273        Smear_Model_20(Cyl_PolyRadius,s.coefW,s.xW,s.yW,s.resW)
274
275        return(0)
276End
277
278
279
280//// experimental threaded version...
281// don't try to thread the smeared calculation, it's good enough
282// to thread the unsmeared version
283
284//threaded version of the function
285ThreadSafe Function Cyl_PolyRadius_T(cw,yw,xw,p1,p2)
286        WAVE cw,yw,xw
287        Variable p1,p2
288       
289#if exists("Cyl_PolyRadiusX")
290        yw[p1,p2] = Cyl_PolyRadiusX(cw,xw)
291#else
292        yw[p1,p2] = fCyl_PolyRadius(cw,xw)
293#endif
294
295        return 0
296End
297
298//
299//  Fit function that is actually a wrapper to dispatch the calculation to N threads
300//
301// nthreads is 1 or an even number, typically 2
302// it doesn't matter if npt is odd. In this case, fractional point numbers are passed
303// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points
304// and the points "2.5" and "3.5" evaluate correctly as 2 and 3
305//
306Function Cyl_PolyRadius(cw,yw,xw) : FitFunc
307        Wave cw,yw,xw
308       
309        Variable npt=numpnts(yw)
310        Variable i,nthreads= ThreadProcessorCount
311        variable mt= ThreadGroupCreate(nthreads)
312
313//      Variable t1=StopMSTimer(-2)
314       
315        for(i=0;i<nthreads;i+=1)
316        //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
317                ThreadStart mt,i,Cyl_PolyRadius_T(cw,yw,xw,(i*npt/nthreads),((i+1)*npt/nthreads-1))
318        endfor
319
320        do
321                variable tgs= ThreadGroupWait(mt,100)
322        while( tgs != 0 )
323
324        variable dummy= ThreadGroupRelease(mt)
325       
326//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
327       
328        return(0)
329End
Note: See TracBrowser for help on using the repository browser.