source: sans/Dev/trunk/NCNR_User_Procedures/Common/Smear_2D.ipf @ 891

Last change on this file since 891 was 838, checked in by srkline, 11 years ago

adding two new help files for real space modeling and description of the 2D resolution function

cleaning up the FFT routines for addition as a beta operation
File size: 16.5 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4
5
6// there are utlity functions here (commented out) for plotting the 2D resolution function
7// at any point on the detector
8
9
10// I can only calculate (nord) points in an AAO fashion, so threading the function calculation
11// doesn't help. Even if I can rewrite this to calculate nord*nord AAO, that will typically be 10*10=100
12// and that is not enough points to thread with any benefit
13
14// but the calculation of each q value is independent, so I can split the total number of q-values between processors
15//
16// -- but I must be careful to pass this a function that is not already threaded!
17// --- BUT - I can't pass either a function reference, OR a structure to a thread!
18// ---- So sadly, each 2D resolution calculation must be threaded by hand.
19//
20// See the Sphere_2D function for an example of how to thread the smearing
21//
22//
23// SRK May 2010
24
25//
26// NOTE: there is a definition of FindTheta = -1* FindPhi() that is a duplicate of what is in RawWindowHook.ipf
27// and should eventually be in a common location for analysis and reduction packages
28//
29
30
31//// this is the completely generic 2D smearing, not threaded, but takes FUNCREF and STRUCT parameters
32//
33// uses resolution ellipse defined perpendicular "Y" and parallel "X",
34// rotating the ellipse into its proper orientaiton based on qx,qy
35//
36// 5 gauss points is not enough - it gives artifacts as a function of phi
37// 10 gauss points is minimally sufficient
38// 20 gauss points are needed if lots of oscillations (just like in 1D)
39// even more may be necessary for highly peaked functions
40//
41//
42Function Smear_2DModel_PP(fcn,s,nord)
43        FUNCREF SANS_2D_ModelAAO_proto fcn
44        Struct ResSmear_2D_AAOStruct &s
45        Variable nord
46       
47        String weightStr,zStr
48
49        Variable ii,jj,kk,num
50        Variable qx,qy,qz,qval,fs
51        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
52       
53        Variable a,b,c,normFactor,phi,theta,maxSig,numStdDev=3
54       
55        switch(nord)   
56                case 5:         
57                        weightStr="gauss5wt"
58                        zStr="gauss5z"
59                        if (WaveExists($weightStr) == 0)
60                                Make/O/D/N=(nord) $weightStr,$zStr
61                                Make5GaussPoints($weightStr,$zStr)     
62                        endif
63                        break                           
64                case 10:               
65                        weightStr="gauss10wt"
66                        zStr="gauss10z"
67                        if (WaveExists($weightStr) == 0)
68                                Make/O/D/N=(nord) $weightStr,$zStr
69                                Make10GaussPoints($weightStr,$zStr)     
70                        endif
71                        break                           
72                case 20:               
73                        weightStr="gauss20wt"
74                        zStr="gauss20z"
75                        if (WaveExists($weightStr) == 0)
76                                Make/O/D/N=(nord) $weightStr,$zStr
77                                Make20GaussPoints($weightStr,$zStr)     
78                        endif
79                        break
80                default:                                                       
81                        Abort "Smear_2DModel_PP called with invalid nord value"                                 
82        endswitch
83       
84        Wave/Z wt = $weightStr
85        Wave/Z xi = $zStr
86       
87/// keep these waves local
88//      Make/O/D/N=1 yPtW
89        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW
90       
91        // now just loop over the points as specified
92        num=numpnts(s.xw[0])
93       
94        answer=0
95       
96        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt
97        Variable qperp_pt,phi_prime,q_prime
98
99        Variable t1=StopMSTimer(-2)
100
101        //loop over q-values
102        for(ii=0;ii<num;ii+=1)
103       
104//              if(mod(ii, 1000 ) == 0)
105//                      Print "ii= ",ii
106//              endif
107               
108                qx = s.xw[0][ii]
109                qy = s.xw[1][ii]
110                qz = s.qz[ii]
111                qval = sqrt(qx^2+qy^2+qz^2)
112                spl = s.sQpl[ii]
113                spp = s.sQpp[ii]
114                fs = s.fs[ii]
115               
116                normFactor = 2*pi*spl*spp
117               
118                phi = -1*FindTheta(qx,qy)               //Findtheta is an exact duplicate of FindPhi() * -1
119               
120                apl = -numStdDev*spl + qval             //parallel = q integration limits
121                bpl = numStdDev*spl + qval
122///             app = -numStdDev*spp + phi              //perpendicular = phi integration limits (WRONG)
123///             bpp = numStdDev*spp + phi
124                app = -numStdDev*spp + 0                //q_perp = 0
125                bpp = numStdDev*spp + 0
126               
127                //make sure the limits are reasonable.
128                if(apl < 0)
129                        apl = 0
130                endif
131                // do I need to specially handle limits when phi ~ 0?
132       
133               
134                sumOut = 0
135                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
136///                     phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2
137                        qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2         //this is now q_perp
138
139                        sumIn=0
140                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
141
142                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2
143                               
144///                             FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
145
146                                // find QxQy given Qpl and Qperp on the grid
147                                //
148                                q_prime = sqrt(qpl_pt^2+qperp_pt^2)
149                                phi_prime = phi + qperp_pt/qpl_pt
150                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt)
151                               
152                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not
153                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
154
155                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) )
156///                             res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) )
157                                res_tot[kk] /= normFactor
158//                              res_tot[kk] *= fs
159
160                        endfor
161                       
162                        fcn(s.coefW,fcnRet,xptw,yptw)                   //calculate nord pts at a time
163                       
164                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0]
165                        fcnRet *= wt[jj]*wt*res_tot
166                        //
167                        answer += (bpl-apl)/2.0*sum(fcnRet)             //
168                endfor
169
170                answer *= (bpp-app)/2.0
171                s.zw[ii] = answer
172        endfor
173       
174        Variable elap = (StopMSTimer(-2) - t1)/1e6
175        Print "elapsed time = ",elap
176       
177        return(0)
178end
179
180
181// this is generic, but I need to declare the Cylinder2D function threadsafe
182// and this calculation is significantly slower than the manually threaded calculation
183// if the function is fast to calculate. Once the function has polydispersity on 2 or more parameters
184// then this AAO calculation and the manual threading are both painfully slow, and more similar in execution time
185//
186// For 128x128 data, and 10x10 smearing, using 25 points for each polydispersity (using DANSE code)
187//      Successively making more things polydisperse gives the following timing (in seconds):
188//
189//                              monodisp                sigR            sigTheta                sigPhi
190// manual THR           1.1             3.9                     74                      1844
191//      this AAO                        8.6             11.4                    104             1930
192//
193// and using 5 points for each polydispersity: (I see no visual difference)
194// manual THR           1.1             1.6                     3.9             16.1
195//
196// so clearly -- use 5 points for the polydispersities, unless there's a good reason not to - and
197// certainly for a survey, it's the way to go.
198//
199Function Smear_2DModel_PP_AAO(fcn,s,nord)
200        FUNCREF SANS_2D_ModelAAO_proto fcn
201        Struct ResSmear_2D_AAOStruct &s
202        Variable nord
203       
204        String weightStr,zStr
205
206        Variable ii,jj,kk,num
207        Variable qx,qy,qz,qval,fs
208        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut
209       
210        Variable a,b,c,normFactor,phi,theta,maxSig,numStdDev=3
211       
212        switch(nord)   
213                case 5:         
214                        weightStr="gauss5wt"
215                        zStr="gauss5z"
216                        if (WaveExists($weightStr) == 0)
217                                Make/O/D/N=(nord) $weightStr,$zStr
218                                Make5GaussPoints($weightStr,$zStr)     
219                        endif
220                        break                           
221                case 10:               
222                        weightStr="gauss10wt"
223                        zStr="gauss10z"
224                        if (WaveExists($weightStr) == 0)
225                                Make/O/D/N=(nord) $weightStr,$zStr
226                                Make10GaussPoints($weightStr,$zStr)     
227                        endif
228                        break                           
229                case 20:               
230                        weightStr="gauss20wt"
231                        zStr="gauss20z"
232                        if (WaveExists($weightStr) == 0)
233                                Make/O/D/N=(nord) $weightStr,$zStr
234                                Make20GaussPoints($weightStr,$zStr)     
235                        endif
236                        break
237                default:                                                       
238                        Abort "Smear_2DModel_PP called with invalid nord value"                                 
239        endswitch
240       
241        Wave/Z wt = $weightStr
242        Wave/Z xi = $zStr
243       
244/// keep these waves local
245//      Make/O/D/N=1 yPtW
246        Make/O/D/N=(nord*nord) fcnRet,xptW,res_tot,yptW,wts
247        Make/O/D/N=(nord) phi_pt,qpl_pt,qperp_pt
248               
249        // now just loop over the points as specified
250        num=numpnts(s.xw[0])
251       
252        answer=0
253       
254        Variable spl,spp,apl,app,bpl,bpp
255        Variable phi_prime,q_prime
256
257        Variable t1=StopMSTimer(-2)
258
259        //loop over q-values
260        for(ii=0;ii<num;ii+=1)
261       
262//              if(mod(ii, 1000 ) == 0)
263//                      Print "ii= ",ii
264//              endif
265               
266                qx = s.xw[0][ii]
267                qy = s.xw[1][ii]
268                qz = s.qz[ii]
269                qval = sqrt(qx^2+qy^2+qz^2)
270                spl = s.sQpl[ii]
271                spp = s.sQpp[ii]
272                fs = s.fs[ii]
273               
274                normFactor = 2*pi*spl*spp
275               
276                phi = -1*FindTheta(qx,qy)               //Findtheta is an exact duplicate of FindPhi() * -1
277               
278                apl = -numStdDev*spl + qval             //parallel = q integration limits
279                bpl = numStdDev*spl + qval
280///             app = -numStdDev*spp + phi              //perpendicular = phi integration limits (WRONG)
281///             bpp = numStdDev*spp + phi
282                app = -numStdDev*spp + 0                //q_perp = 0
283                bpp = numStdDev*spp + 0
284               
285                //make sure the limits are reasonable.
286                if(apl < 0)
287                        apl = 0
288                endif
289                // do I need to specially handle limits when phi ~ 0?
290       
291               
292//              sumOut = 0
293                for(jj=0;jj<nord;jj+=1)         // call phi the "outer'
294//                      phi_pt[jj] = (xi[jj]*(bpp-app)+app+bpp)/2
295                        qperp_pt[jj] = (xi[jj]*(bpp-app)+app+bpp)/2             //this is now q_perp
296                       
297//                      sumIn=0
298                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl
299
300                                qpl_pt[kk] = (xi[kk]*(bpl-apl)+apl+bpl)/2
301                               
302///                             FindQxQy(qpl_pt[kk],phi_pt[jj],qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi
303                               
304                                // find QxQy given Qpl and Qperp on the grid
305                                //
306                                q_prime = sqrt(qpl_pt[kk]^2+qperp_pt[jj]^2)
307                                phi_prime = phi + qperp_pt[jj]/qpl_pt[kk]
308                                FindQxQy(q_prime,phi_prime,qx_pt,qy_pt)
309                                                               
310                                yPtw[nord*jj+kk] = qy_pt                                        //phi is the same in this loop, but qy is not
311                                xPtW[nord*jj+kk] = qx_pt                                        //qx is different here too, as we're varying Qpl
312                               
313                                res_tot[nord*jj+kk] = exp(-0.5*( (qpl_pt[kk]-qval)^2/spl/spl + (qperp_pt[jj])^2/spp/spp ) )
314//                              res_tot[nord*jj+kk] = exp(-0.5*( (qpl_pt[kk]-qval)^2/spl/spl + (phi_pt[jj]-phi)^2/spp/spp ) )
315                                res_tot[nord*jj+kk] /= normFactor
316//                              res_tot[kk] *= fs
317
318                                //weighting
319                                wts[nord*jj+kk] = wt[jj]*wt[kk]
320                        endfor
321                       
322                endfor
323               
324                fcn(s.coefW,fcnRet,xptw,yptw)                   //calculate nord*nord pts at a time
325               
326                fcnRet *= wts*res_tot
327                //
328                answer = (bpl-apl)/2.0*sum(fcnRet)              // get the sum, normalize to parallel direction
329                answer *= (bpp-app)/2.0                                         // and normalize to perpendicular direction
330               
331                s.zw[ii] = answer
332        endfor
333       
334        Variable elap = (StopMSTimer(-2) - t1)/1e6
335        Print "elapsed time = ",elap
336       
337        return(0)
338end
339
340
341
342
343
344
345//phi is defined from +x axis, proceeding CCW around [0,2Pi]
346//rotate the resolution function by theta,  = -phi
347//
348// this is only different by (-1) from FindPhi
349// I'd just call FindPhi, but it's awkward to include
350//
351Threadsafe Function FindTheta(vx,vy)
352        variable vx,vy
353
354       
355        return(-1 * FindPhi(vx,vy))
356end
357
358
359
360       
361//////////////////////////////////////////////////////////
362//////////////////////////////////////////////////////////
363/// utility functions to test the calculation of the resolution function and
364//       plotting of it for visualization
365//
366// -- these need to have the reduction package included for it to compile
367//
368// this works with the raw 2D data, calculating resolution in place
369// type is "RAW"
370// xx,yy are in detector coordinates
371//
372// call as:
373// PlotResolution_atPixel("RAW",pcsr(A),qcsr(A))
374//
375// then:
376// Display;AppendImage res;AppendMatrixContour res;ModifyContour res labels=0,autoLevels={*,*,3}
377//
378
379Function PlotResolution_atPixel(type,xx,yy)
380        String type
381        Variable xx,yy
382       
383///from QxQyExport
384        String destStr="",typeStr=""
385        Variable step=1,refnum
386        destStr = "root:Packages:NIST:"+type
387       
388        //must select the linear_data to export
389        NVAR isLog = $(destStr+":gIsLogScale")
390        if(isLog==1)
391                typeStr = ":linear_data"
392        else
393                typeStr = ":data"
394        endif
395       
396        NVAR pixelsX = root:myGlobals:gNPixelsX
397        NVAR pixelsY = root:myGlobals:gNPixelsY
398       
399        Wave data=$(destStr+typeStr)
400        WAVE intw=$(destStr + ":integersRead")
401        WAVE rw=$(destStr + ":realsRead")
402        WAVE/T textw=$(destStr + ":textRead")
403       
404//      Duplicate/O data,qx_val,qy_val,z_val,qval,qz_val,phi,r_dist
405        Variable qx_val,qy_val,z_val,qval,qz_val,phi,r_dist
406
407        Variable xctr,yctr,sdd,lambda,pixSize
408        xctr = rw[16]
409        yctr = rw[17]
410        sdd = rw[18]
411        lambda = rw[26]
412        pixSize = rw[13]/10             //convert mm to cm (x and y are the same size pixels)
413       
414        qx_val = CalcQx(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)                //+1 converts to detector coordinate system
415        qy_val = CalcQy(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
416       
417        Variable L2 = rw[18]
418        Variable BS = rw[21]
419        Variable S1 = rw[23]
420        Variable S2 = rw[24]
421        Variable L1 = rw[25]
422        Variable lambdaWidth = rw[27]   
423        Variable usingLenses = rw[28]           //new 2007
424       
425        Variable vz_1 = 3.956e5         //velocity [cm/s] of 1 A neutron
426        Variable g = 981.0                              //gravity acceleration [cm/s^2]
427        Variable m_h    = 252.8                 // m/h [=] s/cm^2
428
429        Variable acc,ssd,lambda0,yg_d,qstar
430               
431        G = 981.  //!   ACCELERATION OF GRAVITY, CM/SEC^2
432        acc = vz_1              //      3.956E5 //!     CONVERT WAVELENGTH TO VELOCITY CM/SEC
433        SDD = L2        *100    //1317
434        SSD = L1        *100    //1627          //cm
435        lambda0 = lambda                //              15
436        YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2
437        Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG_d
438                Print "Gravity q* = ",-2*pi/lambda0*2*yg_d/sdd
439        qstar = -2*pi/lambda0*2*yg_d/sdd
440       
441
442// the gravity center is not the resolution center
443// gravity center = beam center
444// resolution center = offset y = dy + (2)*yg_d
445
446
447        qval = CalcQval(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
448        qz_val = CalcQz(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)
449        phi = FindPhi( pixSize*((xx+1)-xctr) , pixSize*((yy+1)-yctr)+(2)*yg_d)          //(dx,dy+yg_d)
450        r_dist = sqrt(  (pixSize*((xx+1)-xctr))^2 +  (pixSize*((yy+1)-yctr)+(2)*yg_d)^2 )               //radial distance from ctr to pt
451       
452        Print pixSize*((yy+1)-yctr),pixSize*((yy+1)-yctr)+(2)*yg_d
453       
454//      Redimension/N=(pixelsX*pixelsY) qz_val,qval,phi,r_dist
455        //everything in 1D now
456//      Duplicate/O qval SigmaQX,SigmaQY,fsubS
457        Variable SigmaQX,SigmaQY,fsubS
458
459        //Two parameters DDET and APOFF are instrument dependent.  Determine
460        //these from the instrument name in the header.
461        //From conversation with JB on 01.06.99 these are the current good values
462        Variable DDet
463        NVAR apOff = root:myGlobals:apOff               //in cm
464        DDet = rw[10]/10                        // header value (X) is in mm, want cm here
465
466        Variable ret1,ret2,ret3,del_r
467        del_r = rw[10]
468       
469        get2DResolution(qval,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,ret1,ret2,ret3)
470        SigmaQX = ret1 
471        SigmaQY = ret2 
472        fsubs = ret3   
473/////
474
475//      Variable theta,phi,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0
476        Variable theta,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0,maxSig,nStdDev=3,normFactor
477        Variable qx_ret,qy_ret
478       
479//      theta = FindPhi(qx_val,qy_val)
480// need to rotate properly - theta is defined a starting from +y axis, moving CW
481// we define phi starting from +x and moving CCW
482        theta = -phi                    //seems to give the right behavior...
483       
484       
485        Print qx_val,qy_val,qval
486        Print "phi, theta",phi,theta
487       
488//      FindQxQy(qval,phi,qx_ret,qy_ret)
489       
490        sx = SigmaQx
491        sy = sigmaQy
492        x0 = qx_val
493        y0 = qy_val
494       
495        a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy)
496        b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy)
497        c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy)
498       
499        normFactor = pi/sqrt(a*c-b*b)
500
501        Make/O/D/N=(num,num) res
502        // so the resolution function 'looks right' on a 2D plot - otherwise it always looks like a circle
503        maxSig = max(sx,sy)
504        Setscale/I x -nStdDev*maxSig+x0,nStdDev*maxSig+x0,res
505        Setscale/I y -nStdDev*maxSig+y0,nStdDev*maxSig+y0,res
506/////   Setscale/I x -nStdDev*sx+x0,nStdDev*sx+x0,res
507/////   Setscale/I y -nStdDev*sy+y0,nStdDev*sy+y0,res
508       
509        Variable xPt,yPt,delx,dely,offx,offy
510        delx = DimDelta(res,0)
511        dely = DimDelta(res,1)
512        offx = DimOffset(res,0)
513        offy = DimOffset(res,1)
514
515        Print "sx,sy = ",sx,sy
516        for(ii=0;ii<num;ii+=1)
517                xPt = offx + ii*delx
518                for(jj=0;jj<num;jj+=1)
519                        yPt = offy + jj*dely
520                        res[ii][jj] = exp(-1*(a*(xPt-x0)^2 + 2*b*(xPt-x0)*(yPt-y0) + c*(yPt-y0)^2))
521                endfor
522        endfor 
523        res /= normFactor
524       
525        //Print sum(res,-inf,inf)*delx*dely
526        if(WaveExists($"coef")==0)
527                Make/O/D/N=6 coef
528        endif
529        Wave coef=coef
530        coef[0] = 1
531        coef[1] = qx_val
532        coef[2] = qy_val
533        coef[3] = sx
534        coef[4] = sy
535        coef[5] = theta
536
537//      Variable t1=StopMSTimer(-2)
538
539//
540        do2dIntegrationGauss(-nStdDev*maxSig+x0,nStdDev*maxSig+x0,-nStdDev*maxSig+y0,nStdDev*maxSig+y0)
541//
542
543//      Variable elap = (StopMSTimer(-2) - t1)/1e6
544//      Print "elapsed time = ",elap
545//      Print "time for 16384 = (minutes)",16384*elap/60
546        return(0)
547End
548
549
550// this is called each time to integrate the gaussian
551Function do2dIntegrationGauss(xMin,xMax,yMin,yMax)
552        Variable xMin,xMax,yMin,yMax
553       
554        Variable/G globalXmin=xMin
555        Variable/G globalXmax=xMax
556        Variable/G globalY
557                       
558        Variable result=Integrate1d(Gauss2DFuncOuter,yMin,yMax,2,5)       
559        KillVariables/z globalXmax,globalXmin,globalY
560        print "integration of 2D = ",result
561End
562
563Function Gauss2DFuncOuter(inY)
564        Variable inY
565       
566        NVAR globalXmin,globalXmax,globalY
567        globalY=inY
568       
569        return integrate1D(Gauss2DFuncInner,globalXmin,globalXmax,2,5)         
570End
571
572Function Gauss2DFuncInner(inX)
573        Variable inX
574       
575        NVAR globalY
576        Wave coef=coef
577       
578        return Gauss2D_theta(coef,inX,GlobalY)
579End
580
581Function Gauss2D_theta(w,x,y)
582        Wave w
583        Variable x,y
584       
585        Variable val,a,b,c
586        Variable scale,x0,y0,sx,sy,theta,normFactor
587       
588        scale = w[0]
589        x0 = w[1]
590        y0 = w[2]
591        sx = w[3]
592        sy = w[4]
593        theta = w[5]
594       
595        a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy)
596        b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy)
597        c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy)
598       
599        val = exp(-1*(a*(x-x0)^2 + 2*b*(x-x0)*(y-y0) + c*(y-y0)^2))
600       
601        normFactor = pi/sqrt(a*c-b*b)
602       
603        return(scale*val/normFactor)
604end
605
606////////////////////////////////////////////////////////////
607////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.