source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/Cylinder_PolyRadius_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: 9.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4#include "Cylinder_v40"
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 (A^-1) for model: "
23        Prompt qmax "Enter maximum q-value (A^-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,1e-6,6.3e-6,0.01}
28        make/o/t parameters_cypr = {"scale","radius (A)","length (A)","polydispersity of Radius","SLD cylinder (A^-2)","SLD solvent (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 (A\\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","parameters_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,1e-6,6.3e-6,0.01}
57        make/o/t smear_parameters_cypr = {"scale","radius (A)","length (A)","polydispersity of Radius","SLD cylinder (A^-2)","SLD solvent (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 (A\\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","smear_parameters_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)
84//Function Cyl_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)
93//End
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] sld cylinder (A^-2)
105        //[5] sld solvent
106        //[6] background (cm^-1)
107        Variable scale,radius,pd,delrho,bkg,zz,length,sldc,slds
108        scale = w[0]
109        radius = w[1]
110        length = w[2]
111        pd = w[3]
112        sldc = w[4]
113        slds = w[5]
114        bkg = w[6]
115       
116        delrho = sldc - slds
117        zz = (1/pd)^2-1
118//
119// the OUTPUT form factor is <f^2>/Vavg [cm-1]
120//
121// local variables
122        Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq
123        Variable answer,zp1,zp2,zp3,vpoly
124        String weightStr,zStr
125       
126//      nord = 5       
127//      weightStr = "gauss5wt"
128//      zStr = "gauss5z"
129        nord = 20
130        weightStr = "gauss20wt"
131        zStr = "gauss20z"
132
133//      if wt,z waves don't exist, create them
134// 5 Gauss points (not enough for cylinder radius = high q oscillations)
135// use 20 Gauss points
136        if (WaveExists($weightStr) == 0) // wave reference is not valid,
137                Make/D/N=(nord) $weightStr,$zStr
138                Wave wtGau = $weightStr
139                Wave zGau = $zStr               // wave references to pass
140                Make20GaussPoints(wtGau,zGau)   
141                //Make5GaussPoints(wtGau,zGau) 
142//      //                  printf "w[0],z[0] = %g %g\r", wtGau[0],zGau[0]
143        else
144                if(exists(weightStr) > 1)
145                         Abort "wave name is already in use"    // execute if condition is false
146                endif
147                Wave wtGau = $weightStr
148                Wave zGau = $zStr
149//      //          printf "w[0],z[0] = %g %g\r", wtGau[0],zGau[0]     
150        endif
151
152// set up the integration
153// end points and weights
154// limits are technically 0-inf, but wisely choose non-zero region of distribution
155        Variable range=3.4              //multiples of the std. dev. fom the mean
156        a = radius*(1-range*pd)
157        if (a<0)
158                a=0             //otherwise numerical error when pd >= 0.3, making a<0
159        endif
160        If(pd>0.3)
161                range = 3.4 + (pd-0.3)*18
162        Endif
163        b = radius*(1+range*pd) // is this far enough past avg radius?
164//      printf "a,b,ravg = %g %g %g\r", a,b,radius
165        va =a
166        vb =b
167
168// evaluate at Gauss points
169        // remember to index from 0,size-1     
170        qq = x          //current x point is the q-value for evaluation
171        summ = 0.0              // initialize integral
172   ii=0
173   do
174   //printf "top of nord loop, i = %g\r",i
175        // Using 5 Gauss points         
176                zi = ( zGau[ii]*(vb-va) + vb + va )/2.0         
177                yyy = wtGau[ii] * rad_kernel(qq,radius,length,zz,sldc,slds,zi)
178                summ = yyy + summ
179                ii+=1
180        while (ii<nord)                         // end of loop over quadrature points
181//   
182// calculate value of integral to return
183   answer = (vb-va)/2.0*summ
184     
185//  contrast^2 is included in integration rad_kernel
186//      answer *= delrho*delrho
187//normalize by polydisperse volume
188// now volume depends on polydisperse RADIUS - so normalize by the second moment
189// 2nd moment = (zz+2)/(zz+1)
190        vpoly = Pi*(radius)^2*length*(zz+2)/(zz+1)
191//Divide by vol, since volume has been "un-normalized" out
192        answer /= vpoly
193//convert to [cm-1]
194        answer *= 1.0e8
195//scale
196        answer *= scale
197// add in the background
198        answer += bkg
199
200        Return (answer)
201End             //End of function PolyRadCylForm()
202
203Function rad_kernel(qw,ravg,len,zz,sldc,slds,rad)
204        Variable qw,ravg,len,zz,sldc,slds,rad
205       
206        Variable Pq,vcyl,dr
207       
208        //calculate the orientationally averaged P(q) for the input rad
209        //this is correct - see K&C (1983) or Lin &Tsao JACryst (1996)29 170.
210        Make/O/D/n=6 kernpar
211        Wave kp = kernpar
212        kp[0] = 1               //scale fixed at 1
213        kp[1] = rad
214        kp[2] = len
215        kp[3] = sldc
216        kp[4] = slds
217        kp[5] = 0               //bkg fixed at 0
218       
219#if exists("CylinderFormX")
220        Pq = CylinderFormX(kp,qw)
221#else
222        Pq = fCylinderForm(kp,qw)
223#endif
224       
225        // undo the normalization that CylinderForm does
226        vcyl=Pi*rad*rad*len
227        Pq *= vcyl
228        //un-convert from [cm-1]
229        Pq /= 1.0e8
230       
231        // calculate normalized distribution at len value
232        dr = Schulz_Point_pr(rad,ravg,zz)
233       
234        return (Pq*dr) 
235End
236
237Function Schulz_Point_pr(x,avg,zz)
238        Variable x,avg,zz
239       
240        Variable dr
241       
242        dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1))
243       
244        return (exp(dr))
245End
246
247//wrapper to calculate the smeared model as an AAO-Struct
248// fills the struct and calls the ususal function with the STRUCT parameter
249//
250// used only for the dependency, not for fitting
251//
252Function fSmearedCyl_PolyRadius(coefW,yW,xW)
253        Wave coefW,yW,xW
254       
255        String str = getWavesDataFolder(yW,0)
256        String DF="root:"+str+":"
257       
258        WAVE resW = $(DF+str+"_res")
259       
260        STRUCT ResSmearAAOStruct fs
261        WAVE fs.coefW = coefW   
262        WAVE fs.yW = yW
263        WAVE fs.xW = xW
264        WAVE fs.resW = resW
265       
266        Variable err
267        err = SmearedCyl_PolyRadius(fs)
268       
269        return (0)
270End
271
272// this is all there is to the smeared calculation!
273Function SmearedCyl_PolyRadius(s) :FitFunc
274        Struct ResSmearAAOStruct &s
275
276//      the name of your unsmeared model (AAO) is the first argument
277        Smear_Model_20(Cyl_PolyRadius,s.coefW,s.xW,s.yW,s.resW)
278
279        return(0)
280End
281
282//
283//  Fit function that is actually a wrapper to dispatch the calculation to N threads
284//
285// nthreads is 1 or an even number, typically 2
286// it doesn't matter if npt is odd. In this case, fractional point numbers are passed
287// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points
288// and the points "2.5" and "3.5" evaluate correctly as 2 and 3
289//
290Function Cyl_PolyRadius(cw,yw,xw) : FitFunc
291        Wave cw,yw,xw
292
293//      Variable t1=StopMSTimer(-2)
294
295///////NON-THREADED VERSION ///////
296//#if exists("Cyl_PolyRadiusX")
297//      yw = Cyl_PolyRadiusX(cw,xw)
298//#else
299//      yw = fCyl_PolyRadius(cw,xw)
300//#endif
301
302///// THREADED VERSION NEEDS Igor 6.10B04 or higher to avoid crashes //////     
303#if exists("Cyl_PolyRadiusX")
304
305        Variable npt=numpnts(yw)
306        Variable i,nthreads= ThreadProcessorCount
307        variable mt= ThreadGroupCreate(nthreads)
308
309        for(i=0;i<nthreads;i+=1)
310        //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
311                ThreadStart mt,i,Cyl_PolyRadius_T(cw,yw,xw,(i*npt/nthreads),((i+1)*npt/nthreads-1))
312        endfor
313
314        do
315                variable tgs= ThreadGroupWait(mt,100)
316        while( tgs != 0 )
317
318        variable dummy= ThreadGroupRelease(mt)
319       
320#else
321                yw = fCyl_PolyRadius(cw,xw)             //the Igor, non-XOP, non-threaded calculation, messy to make ThreadSafe
322#endif
323
324//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
325
326        return(0)
327End
328
329//// experimental threaded version...
330// don't try to thread the smeared calculation, it's good enough
331// to thread the unsmeared version
332
333//threaded version of the function
334ThreadSafe Function Cyl_PolyRadius_T(cw,yw,xw,p1,p2)
335        WAVE cw,yw,xw
336        Variable p1,p2
337       
338#if exists("Cyl_PolyRadiusX")                   //this check is done in the calling function, simply hide from compiler
339        yw[p1,p2] = Cyl_PolyRadiusX(cw,xw)
340#else
341        yw[p1,p2] = fCyl_PolyRadius(cw,xw)
342#endif
343
344        return 0
345End
Note: See TracBrowser for help on using the repository browser.