source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Cylinder_Fills.ipf @ 800

Last change on this file since 800 was 800, checked in by srkline, 12 years ago

Changed the 2D resolution calculation to include gravity. This seems to work, but will still require some testing. It requires no modifications to the smearing calculation, still using parallel and perpendicular directions.

Added a Fake Update feature to the RealTime? update. There are specific, separate instructions for how to set this up. Usef (possibly) for summer schools or demos.

Adjusted the 2D MonteCarlo? simulation to a corrected beam center. The Gauss Peak was not symmetric around the old beam center.

Corrected the AAO resolution smearing functions to work with cursors.

File size: 21.0 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3// VERY incomplete methods for filling volumes with connected cylinders
4//
5// most is still being worked out and verified in 2D
6//
7
8
9
10//These macros may or may not work - designed to help automate the
11//process of calculating a series of configurations.
12
13
14
15Proc ConnectDots3D(maxNumConn)
16        Variable maxNumConn=1
17       
18        Variable num=numpnts(x3d)
19        Make/O/N=(num) numConnection3D=0
20        Make/O/T/N=(num) connectedTo3D=""
21        fConnectDots3D(maxNumConn)
22end
23
24// connect the dots, no periodic boundary conditions
25//thickness is set to 1
26Function fConnectDots3D(maxNumConn)
27        Variable maxNumConn
28       
29       
30        WAVE x3d=x3d
31        WAVE y3d=y3d
32        WAVE z3d=z3d
33        WAVE numConnection3D = numConnection3D
34        WAVE/T connectedTo3D = connectedTo3D
35        WAVE mat=mat
36        Variable num=numpnts(x3d),ii=0
37        Variable matSize=DimSize(mat,0),jj
38
39        Variable nnX,nnY,nnZ,nnInd,nnDist               //return values of the nearest neighbor
40        Variable thick = 1
41        do
42                for(ii=0;ii<num;ii+=1)
43                        if(numConnection3D[ii] < maxNumConn)
44                                FindNearestNeighbor3D(x3d,y3d,z3d,ii,maxNumConn,nnX,nnY,nnZ,nnInd,nnDist)
45                               
46                                numConnection3D[nnInd] += 1
47                                numConnection3D[ii] += 1
48                                connectedTo3D[nnInd] += num2str(ii)+","
49                                connectedTo3D[ii] += num2str(nnInd)+","
50                               
51                                ConnectPoints3D(mat, x3d[ii],y3d[ii],z3d[ii],x3d[nnInd],y3d[nnInd],z3d[nnInd],thick,1)          //always fill
52                        endif
53                endfor
54                //without rotation, a very closed network structure tends to form (bias?)
55                // rotate (10?) points so that we're not starting from the same place every time
56                // -- but the "connectedTo" information is LOST...
57                //Rotate 10, x2d,y2d, numConnection,connectedTo
58                //Print sum(numConnection,-inf,inf), num*maxNumConn
59        while(sum(numConnection3D,-inf,inf) < num*maxNumConn)
60
61End
62
63Function FindNearestNeighbor3D(xw,yw,zw,startPt,maxConn,nnX,nnY,nnZ,nnInd,nnDist)
64        Wave xw,yw,zw
65        Variable startPt,maxConn,&nnX,&nnY,&nnZ,&nnInd,&nnDist
66       
67        Variable num,matSize,numToFind
68       
69        //make waves to hold the answers
70        num=numpnts(xw)
71        Make/O/N=(num) distWave3D,distIndex3D
72       
73        //now calculate the distances
74        Variable ii,dist,testPt
75       
76        //this is probably the slowest step by far in the process...(tiny XOP?)
77        for(ii=0;ii<num;ii+=1)
78                distWave3D[ii] = (xw[ii]-xw[startPt])^2 +  (yw[ii]-yw[startPt])^2        + (zw[ii]-zw[startPt])^2               //d^2
79        endfor
80
81        MakeIndex distWave3D distIndex3D                //then distWave[distIndex[ii]] will loop through in order of distance
82        WAVE numConnection3D = numConnection3D
83        WAVE/T connectedTo3D = connectedTo3D
84        WAVE mat = mat
85       
86        for(ii=1;ii<num;ii+=1)          //[0] point is the test point, dist == 0
87                testPt =  distIndex3D[ii]               //index
88                if(numConnection3D[testPt] < maxConn )          //can test pt accept another connection?
89                        if(WhichListItem(num2str(testPt),connectedTo3D[startPt],",",0) == -1) // not already connected to the starting point?
90        //                      Print ii,testPt
91        //                      Printf "nearest point is (%d,%d), dist^2 = %g\r",xw[testPt],yw[testPt],distWave[testPt]
92                                nnX = xw[testPt]                //return values as pass-by-reference
93                                nnY = yw[testPt]
94                                nnZ = zw[testPt]
95                                nnInd = testPt
96                                nnDist = distWave3D[testPt]
97
98                                return(0)               //found a point, return
99                        endif
100                endif
101        endfor
102
103        return (1)              //did not find a point, return error
104End
105
106///!!! fails if the points are along an axis (or even close) - then the box becomes infinitely thin
107Function ConnectPoints3D_old(w, x1,y1,z1,x2,y2,z2,thick,fill)
108        Wave w
109        Variable  x1,y1,z1,x2,y2,z2,thick,fill
110       
111        //find all of the points within dist=thick of the line defined by p1 and p2
112        // and give these points (in w) the value 1
113       
114        // search through the box defined by p1 and p2
115        Variable ii,jj,kk,dist
116       
117        for(ii=min(x1,x2);ii<=max(x1,x2);ii+=1)
118                for(jj=min(y1,y2);jj<=max(y1,y2);jj+=1)
119                        for(kk=min(z1,z2);kk<=max(z1,z2);kk+=1)
120                                //find distance between current point and P1P2
121                                dist = Pt_Line_Distance3D(ii,jj,kk, x1,y1,z1,x2,y2,z2)
122                                //print dist
123                                if(dist<=thick)
124                                        w[ii][jj][kk] = fill            //add the point
125                                endif
126                        endfor
127                endfor
128        endfor
129       
130End
131
132//
133// thickness is given in voxels
134//
135// the search box is expanded by the thickness in each direction (within the box)
136// to account for connections that are along an axis direction
137//
138Function ConnectPoints3D(w, x1,y1,z1,x2,y2,z2,thick,fill)
139        Wave w
140        Variable  x1,y1,z1,x2,y2,z2,thick,fill
141       
142        //find all of the points within dist=thick of the line defined by p1 and p2
143        // and give these points (in w) the value 1
144       
145        // search through the box defined by p1 and p2
146        Variable ii,jj,kk,dist,minX,minY,minZ,maxX,maxY,maxZ
147        Variable boxMax,slop
148        boxMax = DimSize(w,0) - 1               //max index of box
149        slop = thick + 0                // to make the bounding box large enough the get it all
150       
151        minX=min(x1,x2)                 //give me the smaller x
152        minX=max(minX-slop,0)   //give me zero if it's negative
153        minY=min(y1,y2)
154        minY=max(minY-slop,0)
155        minZ=min(z1,z2)
156        minZ=max(minZ-slop,0)
157       
158        maxX=max(x1,x2)                                 //give me the larger x
159        maxX=min(maxX+slop,boxMax)              //give the the box dimension if it's too big
160        maxY=max(y1,y2)
161        maxY=min(maxY+slop,boxMax)
162        maxZ=max(z1,z2)
163        maxZ=min(maxZ+slop,boxMax)
164               
165        for(ii=minX;ii<=maxX;ii+=1)
166                for(jj=minY;jj<=maxY;jj+=1)
167                        for(kk=minZ;kk<=maxZ;kk+=1)
168                                //find distance between current point and P1P2
169                                dist = Pt_Line_Distance3D(ii,jj,kk, x1,y1,z1,x2,y2,z2)
170                                //print dist
171                                if(dist<=thick)
172                                        w[ii][jj][kk] = fill            //add the point
173                                endif
174                        endfor
175                endfor
176        endfor
177       
178End
179
180//solution courtesy of the math forum @ Drexel "Ask Dr. Math"
181Function Pt_Line_Distance3D(xx,yy,zz, x1,y1,z1,x2,y2,z2)
182        Variable xx,yy,zz, x1,y1,z1,x2,y2,z2
183       
184        Variable dist,x3,y3,z3,x4,y4,z4,aa,ab
185       
186        //AB
187        x3 = x2-x1
188        y3 = y2-y1
189        z3 = z2-z1
190        //AP
191        x4 = xx-x1
192        y4 = yy-y1
193        z4 = zz-z1
194       
195        //aa = length of ABxAP
196        aa = sqrt( (y3*z4-z3*y4)^2 + (x3*z4-x4*z3)^2 + (x3*y4-x4*y3)^2 )
197       
198        //ab = length of AB
199        ab = sqrt(x3^2+y3^2+z3^2)
200       
201        dist = aa/ab
202       
203        return(dist)
204End
205
206//////////////////////////
207//// random rod fills - just a voxel at each point, so use T=diameter
208
209// fill a volume with "rods" that are made up of a numer of spheres
210// lined up in a random direction
211//
212// all of the spheres/rods that are added must be connected = "adjacent" to
213// an occupied voxel.
214//
215// this implementation is a random walk, end-to-end connection
216//
217Function ConnectedRodFill(mat,len,num,periodic)
218        Wave mat
219        variable len,num                //length in direction, number of spheres to add
220        Variable periodic               //==1 if periodic boundaries
221       
222        Variable nptx,npty,nptz,ii,fail=0,nAdd
223        Make/O/I/N=3 stVec,dirVec,endVec
224       
225        nptx=DimSize(mat,0)
226        npty=DimSize(mat,1)
227        nptz=DimSize(mat,2)
228       
229        //place first rod randomly
230        stVec[0]=trunc(abs(enoise(nptx)))               //random integer distr betw (0,npt)
231        stVec[1]=trunc(abs(enoise(npty)))
232        stVec[2]=trunc(abs(enoise(nptz)))
233        dirVec[0] = PickDirection()
234        dirVec[1] = PickDirection()
235        dirVec[2] = PickDirection()
236
237        nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
238       
239        stVec = endVec          //end point is the new starting point
240
241        //loop over remaining
242        ii=nAdd
243        do
244
245                        dirVec[0] = PickDirection()
246                        dirVec[1] = PickDirection()
247                        dirVec[2] = PickDirection()
248                        nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
249                        stVec = endVec          //end point is the new starting point
250
251                        ii+=nAdd
252
253                        if(mod(ii,100000) == 0)
254                                print "Point = ",ii
255                        endif
256                       
257        while(ii<num)
258       
259       
260        return(0)
261End
262
263// fill a volume with "rods" that are made up of a numer of spheres
264// lined up in a random direction
265//
266//
267Function UnConnectedRodFill(mat,len,num,periodic)
268        Wave mat
269        variable len,num                //length in direction, number of spheres to add
270        Variable periodic               //==1 if periodic boundaries
271       
272        Variable nptx,npty,nptz,ii,fail=0,nAdd
273        Make/O/I/N=3 stVec,dirVec,endVec
274
275        nptx=DimSize(mat,0)
276        npty=DimSize(mat,1)
277        nptz=DimSize(mat,2)
278       
279        //place first rod randomly
280        stVec[0]=trunc(abs(enoise(nptx)))               //random integer distr betw (0,npt)
281        stVec[1]=trunc(abs(enoise(npty)))
282        stVec[2]=trunc(abs(enoise(nptz)))
283        dirVec[0] = PickDirection()
284        dirVec[1] = PickDirection()
285        dirVec[2] = PickDirection()
286
287        nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
288        //loop over remaining
289        ii=nAdd
290        do
291                //pick a new point
292                stVec[0]=trunc(abs(enoise(nptx)))               //random integer distr betw (0,npt)
293                stVec[1]=trunc(abs(enoise(npty)))
294                stVec[2]=trunc(abs(enoise(nptz)))
295                dirVec[0] = PickDirection()
296                dirVec[1] = PickDirection()
297                dirVec[2] = PickDirection()
298                //add as many spheres as possible in that direction
299                nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
300                ii+=nAdd
301                //Print "point ",ii
302
303        while(ii<num)
304        Print "failures = ",fail
305       
306//      ParseMatrix3D_rho(mat)          //convert to XYZ
307       
308        return(0)
309End
310
311
312Function PickDirection()
313        return( trunc(abs(enoise(3))) - 1       )       //integer -1,0,1
314End
315
316//need to add periodic boundary conditions (done)
317//
318// need to insist on connectivity
319//
320//
321Function AddSpheresInDirection(m,stVec,dirVec,endVec,nTry,periodic)
322        WAVE m,stVec,dirVec,endVec
323        Variable nTry,periodic
324       
325        Variable added,nx,ny,nz,nDimx,nDimy,nDimz
326       
327        nDimx=DimSize(m,0)-1
328        nDimy=DimSize(m,1)-1
329        nDimz=DimSize(m,2)-1
330        added=0
331        nx=stVec[0]
332        ny=stVec[1]
333        nz=stVec[2]
334       
335        do
336                nx+=dirVec[0]           //coordinates on next test point to try to add in d(xyz) direction
337                ny+=dirVec[1]
338                nz+=dirVec[2]
339                //check for out-of bounds - if periodic ==1, reflect to opposite face of cube, else return
340                if( (nx>nDimx) || (nx<0) )
341                        if(periodic)
342                                nx = abs(nDimx+1 - abs(nx))
343                        else
344                                break
345                        endif
346                Endif
347                if( (ny>nDimy) || (ny<0) )
348                        if(periodic)
349                                ny = abs(nDimy+1 - abs(ny))
350                        else
351                                break
352                        endif
353                Endif
354                if( (nz>nDimz) || (nz<0) )
355                        if(periodic)
356                                nz = abs(nDimz+1 - abs(nz))
357                        else
358                                break
359                        endif
360                Endif
361                //add the point or exit
362                if(m[nx][ny][nz])
363                        break           //point is already occupied
364                else
365                        m[nx][ny][nz] = 1
366//                      nx +=1          //we'll back this off before returning
367//                      ny +=1
368//                      nz +=1
369                        added += 1
370                endif
371        while(added<nTry)
372       
373        endVec[0] = nx - dirVec[0]                      //back up to the last point actually added
374        endVec[1] = ny - dirVec[1]     
375        endVec[2] = nz - dirVec[2]     
376        return added
377End
378
379Proc PutRandomCylindersAtPoints(w,num,len,periodic)
380        String w="mat"
381        Variable num=1000,len=10,periodic=1
382        Prompt w,"matrix"
383        Prompt num,"number of starting points"
384        prompt len,"length of cylinders"
385        prompt periodic"1=periodic, 0=non-periodic boundaries"
386       
387        $w=0
388        RandomFill3DMat($w,num)
389        CylindersAtPoints($w,len,periodic)
390       
391        NumberOfPoints()
392end
393
394Proc PutXAxisCylindersAtPoints(w,num,rad,len,periodic,sobol)
395        String w="mat"
396        Variable num=100,rad=20,len=300,periodic=1,sobol=1
397        Prompt w,"matrix"
398        Prompt num,"number of starting points"
399        prompt rad,"radius of cylinders"
400        prompt len,"length of cylinders"
401        prompt periodic,"1=periodic, 0=non-periodic boundaries"
402        Prompt sobol,"1=Sobol, 0=random"
403       
404        $w=0
405        X_CylindersAtPoints($w,num,rad,len,sobol,periodic)
406       
407        NumberOfPoints()
408end
409
410Proc PutXAxisCylindersSquare(w,rad,len,sep)
411        String w="mat"
412        Variable rad=20,len=300,sep=80
413        Prompt w,"matrix"
414        prompt rad,"radius of cylinders"
415        prompt len,"length of cylinders"
416        prompt sep,"center-to-center separation of cylinders"
417       
418        $w=0
419        X_CylindersSquareGrid($w,rad,len,sep)
420       
421        NumberOfPoints()
422end
423
424Proc PutXAxisCylindersHexagonal(w,rad,len,sep)
425        String w="mat"
426        Variable rad=20,len=300,sep=80
427        Prompt w,"matrix"
428        prompt rad,"radius of cylinders"
429        prompt len,"length of cylinders"
430        prompt sep,"center-to-center separation of cylinders"
431       
432        $w=0
433        X_CylindersHexagonalGrid($w,rad,len,sep)
434       
435        NumberOfPoints()
436end
437
438Function CylindersAtPoints(mat,len,periodic)
439        Wave mat
440        variable len            //length in direction
441        Variable periodic               //==1 if periodic boundaries
442       
443        Make/O/B/N=3 stVec,dirVec,endVec
444
445        ParseMatrix3D_rho(mat)
446        Wave x3d=x3d
447        Wave y3d=y3d
448        Wave z3d=z3d
449       
450        Variable ii=0,num=numpnts(x3d)
451        NVAR  grid=root:FFT_T
452       
453        Variable nAdd
454       
455        for(ii=0;ii<num;ii+=1)
456                dirVec[0] = PickDirection()
457                dirVec[1] = PickDirection()
458                dirVec[2] = PickDirection()
459                stVec[0] = x3d[ii]
460                stVec[1] = y3d[ii]
461                stVec[2] = z3d[ii]
462                nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
463                //FillSphereRadius(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],1)
464        endfor
465       
466        return(0)
467End
468
469///////////////
470               
471Function X_CylindersAtPoints(mat,num,rad,len,sobol,periodic)
472        Wave mat
473        variable num,rad,len            //length in direction
474        Variable periodic               //==1 if periodic boundaries
475        Variable sobol          //==1 if sobol selection of points (2D)
476       
477
478        Variable np
479        np = DimSize(mat,0)                     // assumes that all dimensions are the same
480       
481        // fill a 2D plane with points
482        Make/O/B/N=(np,np) plane
483        plane = 0
484        if(sobol)
485                if(WaveExists($"SobolInit") == 0)
486                        Make/O/D/N=2 SobolInit
487                        SobolX(-1,SobolInit)
488                else
489                        SobolPoints2D(plane,num)
490                endif
491        else
492                RandomPoints2D(plane,num)
493        endif
494        // put it in the proper plane of the matrix
495        mat[np/2][][] = plane[q][r]                     // in the YZ plane
496       
497        ParseMatrix3D_rho(mat)
498        Wave x3d=x3d
499        Wave y3d=y3d
500        Wave z3d=z3d
501       
502        Variable ii=0
503        NVAR  grid=root:FFT_T
504       
505        for(ii=0;ii<num;ii+=1)
506                FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1)               //cylinder 1
507        endfor
508       
509        return(0)
510End
511
512
513Function X_CylindersSquareGrid(mat,rad,len,sep)
514        Wave mat
515        variable rad,len                //length of cylinders
516        Variable sep                    // EDGE separation, in same units as cylinder
517       
518       
519        NVAR grid=root:FFT_T
520        Variable np,spacing
521        np = DimSize(mat,0)                     // assumes that all dimensions are the same
522       
523        // fill a 2D plane with points
524        Make/O/B/N=(np,np) plane
525        plane = 0
526
527
528        spacing = round(sep/grid)               // so it's an integer
529        FillPlaneSquareGrid(plane,spacing)
530       
531        // put it in the proper plane of the matrix
532        mat[np/2][][] = plane[q][r]                     // in the YZ plane
533       
534        ParseMatrix3D_rho(mat)
535        Wave x3d=x3d
536        Wave y3d=y3d
537        Wave z3d=z3d
538       
539        Variable ii=0,num
540        num = numpnts(x3d)     
541       
542        for(ii=0;ii<num;ii+=1)
543                FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1)               //cylinder 1
544        endfor
545       
546        return(0)
547End
548
549Function X_CylindersHexagonalGrid(mat,rad,len,sep)
550        Wave mat
551        variable rad,len                //length of cylinders
552        Variable sep                    // EDGE separation, in same units as cylinder
553       
554       
555        NVAR grid=root:FFT_T
556        Variable np,spacing
557        np = DimSize(mat,0)                     // assumes that all dimensions are the same
558       
559        // fill a 2D plane with points
560        Make/O/B/N=(np,np) plane
561        plane = 0
562
563        spacing = round(sep/grid)               // so it's an integer
564        FillPlaneHexagonal(plane,spacing)
565       
566        // put it in the proper plane of the matrix
567        mat[np/2][][] = plane[q][r]                     // in the YZ plane
568       
569        ParseMatrix3D_rho(mat)
570        Wave x3d=x3d
571        Wave y3d=y3d
572        Wave z3d=z3d
573       
574        Variable ii=0,num
575        num = numpnts(x3d)     
576       
577        for(ii=0;ii<num;ii+=1)
578                FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1)               //cylinder 1
579        endfor
580
581//// makes a crude core-shell cylinder 
582//      for(ii=0;ii<num;ii+=1)
583//              FillXCylinder(mat,grid,rad-10,x3d[ii],y3d[ii],z3d[ii],len,3)            //cylinder 1
584//      endfor
585//     
586        return(0)
587End
588
589
590
591Function FillPlaneSquareGrid(w,sp)
592        Wave w
593        Variable sp             //spacing in pixels
594       
595        Variable np,ii,jj
596       
597        np = DimSize(w,0)
598       
599        for(ii=sp;ii<(np-sp/2);ii+=sp)
600                for(jj=sp;jj<(np-sp/2);jj+=sp)
601                        w[ii][jj] = 1
602                endfor
603        endfor
604       
605       
606        return(0)
607End
608
609// sp is the
610Function FillPlaneHexagonal(w,sp)
611        Wave w
612        Variable sp             //spacing in pixels
613       
614        Variable np,ii,jj,xx,yy,fill
615       
616        np = DimSize(w,0)
617       
618        fill = 1
619        // do even columns, then odds
620        for(ii=sp;ii<(np-sp/2);ii+=(sp*sqrt(3)) )
621                for(jj=sp;jj<(np-sp/2);jj+=sp)
622                        w[ii][jj] = fill
623                endfor
624        endfor
625       
626        // now odds
627        for(ii=(sp+sp*sqrt(3)/2);ii<(np-sp/2);ii+=(sp*sqrt(3)))
628                for(jj=sp/2;jj<(np-sp/2);jj+=sp)
629                        w[ii][jj] = fill
630                endfor
631        endfor
632       
633        return(0)
634End
635
636Proc HexagonalCylinders()
637
638        Variable xc,yc,zc,len,rad,aa,cosa,sina,meanRad,pd,meanLen
639       
640        Variable grid=root:FFT_T
641        Variable npt=root:FFT_N
642       
643        xc=npt/2
644        yc=npt/2
645        zc=npt/2
646       
647        meanRad=20
648        pd=0.01
649        meanLen=400
650        len=meanLen
651       
652//      aa=meanRad/grid*2 +40   //grid units, +2 gives 2 units of separation
653        aa=meanRad/grid*2 +2    //grid units, +2 gives 2 units of separation
654        cosa = round(0.5*aa)
655        sina = round(0.866*aa)
656       
657
658        Variable ii,nPass,x1,x2,x3
659        ii=0
660        nPass=1
661
662
663        x1 = Gauss(meanRad-pd*meanRad,meanRad,pd*meanRad)
664        x2 = Gauss(meanRad,meanRad,pd*meanRad)
665        x3 = Gauss(meanRad+pd*meanRad,meanRad,pd*meanRad)
666       
667        Print x1,x2,x3
668//     
669/////////// average over 3 point distribution --
670//      //calculate and sum the three contributions
671//      FastOp mat=0
672//      rad = meanRad-pd*meanRad
673////    len = meanLen + gnoise(pd*meanLen)
674//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
675//      Execute "DoFFT()"
676//      Duplicate/O iBin iBin_1g_128
677//      iBin_1g_128 *= x1
678//     
679//      FastOp mat=0
680//      rad = meanRad
681////    len = meanLen + gnoise(pd*meanLen)
682//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
683//      Execute "DoFFT()"
684//      iBin_1g_128 += x2*iBin
685//     
686//      FastOp mat=0
687//      rad = meanRad+pd*meanRad
688////    len = meanLen + gnoise(pd*meanLen)
689//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
690//      Execute "DoFFT()"
691//      iBin_1g_128 += x3*iBin
692//     
693//      iBin_1g_128 /= (x1+x2+x3)               //renormalize
694//     
695////// end 3 pt average
696
697//////////////////////
698// paired cylinders
699/////////////////////////
700//      //calculate and sum the three contributions
701//      FastOp mat=0
702//      rad = meanRad-pd*meanRad
703////    len = meanLen + gnoise(pd*meanLen)
704//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
705//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
706//      Execute "DoFFT()"
707//      Duplicate/O iBin iBin_2g_256_sep40
708//      iBin_2g_256_sep40 *= x1
709//     
710//      FastOp mat=0
711//      rad = meanRad
712////    len = meanLen + gnoise(pd*meanLen)
713//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
714//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
715//      Execute "DoFFT()"
716//      iBin_2g_256_sep40 += x2*iBin
717//     
718//      FastOp mat=0
719//      rad = meanRad+pd*meanRad
720////    len = meanLen + gnoise(pd*meanLen)
721//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
722//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
723//      Execute "DoFFT()"
724//      iBin_2g_256_sep40 += x3*iBin
725//     
726//      iBin_2g_256_sep40 /= (x1+x2+x3)         //renormalize
727       
728/////////////////////////       
729//
730//
731////
732//      ii = 0 
733//      do
734//              FastOp mat=0
735//              rad = meanRad + gnoise(pd*meanRad)
736//              print ii,rad
737//              FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
738//              Execute "DoFFT()"
739//              if(ii==0)       
740//                      Duplicate/O iBin iBin_1g_256
741//              else
742//                      iBin_1g_256 += iBin
743//              endif
744//              ii+=1
745//      while(ii<nPass)
746//      iBin_1g_256 /= nPass
747//     
748/////   
749//      ii=0
750//      nPass=1
751//      do
752//              FastOp mat=0
753//              rad = meanRad + gnoise(pd*meanRad)
754//              print ii,rad
755//              FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)              //cylinder 1
756//              rad = meanRad + gnoise(pd*meanRad)
757//              print ii,rad
758//              FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
759//              Execute "DoFFT()"
760//              if(ii==0)       
761//                      Duplicate/O iBin iBin_2g_256
762//              else
763//                      iBin_2g_256 += iBin
764//              endif
765//              ii+=1
766//      while(ii<nPass)
767//      iBin_2g_256 /= nPass
768//     
769        ii=0
770        npass=1
771        do
772                FastOp mat=0
773                rad = meanRad + gnoise(pd*meanRad)
774                print ii,rad
775                FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)              //cylinder 1
776                rad = meanRad + gnoise(pd*meanRad)
777                print ii,rad
778                FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)           //cylinder 2
779                rad = meanRad + gnoise(pd*meanRad)
780                print ii,rad
781                FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,1)            //cylinder 3
782                Execute "DoFFT()"
783                if(ii==0)       
784                        Duplicate/O iBin iBin_3g_256
785                else
786                        iBin_3g_256 += iBin
787                endif
788                ii+=1
789        while(ii<nPass)
790        iBin_3g_256 /= nPass
791//     
792        rad = meanRad + gnoise(pd*meanRad)
793        print rad
794        FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)
795        Execute "DoFFT()"
796        Duplicate/O iBin iBin_2
797///     
798        rad = meanRad + gnoise(pd*meanRad)
799        print rad
800        FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,1)
801        Execute "DoFFT()"
802        Duplicate/O iBin iBin_3
803///     
804        rad = meanRad + gnoise(pd*meanRad)
805        print rad
806        FillZCylinder(mat,grid,rad,xc+sina,yc-cosa,zc,len,1)
807        Execute "DoFFT()"
808        Duplicate/O iBin iBin_4
809///     
810        rad = meanRad + gnoise(pd*meanRad)
811        print rad
812        FillZCylinder(mat,grid,rad,xc,yc-aa,zc,len,1)
813        Execute "DoFFT()"
814        Duplicate/O iBin iBin_5
815///     
816        rad = meanRad + gnoise(pd*meanRad)
817        print rad
818        FillZCylinder(mat,grid,rad,xc-sina,yc-cosa,zc,len,1)
819        Execute "DoFFT()"
820        Duplicate/O iBin iBin_6
821///     
822        rad = meanRad + gnoise(pd*meanRad)
823        print rad
824        FillZCylinder(mat,grid,rad,xc-sina,yc+cosa,zc,len,1)
825        Execute "DoFFT()"
826        Duplicate/O iBin iBin_7
827
828end
829
830
831Proc Vary_N_In_Direction()
832
833        FFTEraseMatrixButtonProc("")
834        ConnectedRodFill(mat,2,30000*4,1)
835        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
836        Duplicate/O ival_XOP ival_XOP_2
837        Duplicate/O qval_XOP qval_XOP_2
838       
839        FFTEraseMatrixButtonProc("")
840        ConnectedRodFill(mat,3,30000*4,1)
841        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
842        Duplicate/O ival_XOP ival_XOP_3
843        Duplicate/O qval_XOP qval_XOP_3
844       
845        FFTEraseMatrixButtonProc("")
846        ConnectedRodFill(mat,4,30000*4,1)
847        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
848        Duplicate/O ival_XOP ival_XOP_4
849        Duplicate/O qval_XOP qval_XOP_4
850       
851        FFTEraseMatrixButtonProc("")
852        ConnectedRodFill(mat,5,30000*4,1)
853        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
854        Duplicate/O ival_XOP ival_XOP_5
855        Duplicate/O qval_XOP qval_XOP_5
856       
857        FFTEraseMatrixButtonProc("")
858        ConnectedRodFill(mat,6,30000*4,1)
859        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
860        Duplicate/O ival_XOP ival_XOP_6
861        Duplicate/O qval_XOP qval_XOP_6
862       
863        FFTEraseMatrixButtonProc("")
864        ConnectedRodFill(mat,7,30000*4,1)
865        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
866        Duplicate/O ival_XOP ival_XOP_7
867        Duplicate/O qval_XOP qval_XOP_7
868
869        FFTEraseMatrixButtonProc("")
870        ConnectedRodFill(mat,8,30000*4,1)
871        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
872        Duplicate/O ival_XOP ival_XOP_8
873        Duplicate/O qval_XOP qval_XOP_8
874       
875        FFTEraseMatrixButtonProc("")
876        ConnectedRodFill(mat,9,30000*4,1)
877        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
878        Duplicate/O ival_XOP ival_XOP_9
879        Duplicate/O qval_XOP qval_XOP_9
880       
881        FFTEraseMatrixButtonProc("")
882        ConnectedRodFill(mat,10,30000*4,1)
883        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
884        Duplicate/O ival_XOP ival_XOP_10
885        Duplicate/O qval_XOP qval_XOP_10
886               
887End
Note: See TracBrowser for help on using the repository browser.