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

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

I'm not sure this is a great idea, but I'm putting the FFT / Debye sphere work that I have completed into SVN. It's all really rough - the math, I believe is correct, but the interface if really, really rough. But it's not going to develop without help.

File size: 20.9 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       
403        $w=0
404        X_CylindersAtPoints($w,num,rad,len,sobol,periodic)
405       
406        NumberOfPoints()
407end
408
409Proc PutXAxisCylindersSquare(w,rad,len,sep)
410        String w="mat"
411        Variable rad=20,len=300,sep=80
412        Prompt w,"matrix"
413        prompt rad,"radius of cylinders"
414        prompt len,"length of cylinders"
415        prompt sep,"center-to-center separation of cylinders"
416       
417        $w=0
418        X_CylindersSquareGrid($w,rad,len,sep)
419       
420        NumberOfPoints()
421end
422
423Proc PutXAxisCylindersHexagonal(w,rad,len,sep)
424        String w="mat"
425        Variable rad=20,len=300,sep=80
426        Prompt w,"matrix"
427        prompt rad,"radius of cylinders"
428        prompt len,"length of cylinders"
429        prompt sep,"center-to-center separation of cylinders"
430       
431        $w=0
432        X_CylindersHexagonalGrid($w,rad,len,sep)
433       
434        NumberOfPoints()
435end
436
437Function CylindersAtPoints(mat,len,periodic)
438        Wave mat
439        variable len            //length in direction
440        Variable periodic               //==1 if periodic boundaries
441       
442        Make/O/B/N=3 stVec,dirVec,endVec
443
444        ParseMatrix3D_rho(mat)
445        Wave x3d=x3d
446        Wave y3d=y3d
447        Wave z3d=z3d
448       
449        Variable ii=0,num=numpnts(x3d)
450        NVAR  grid=root:FFT_T
451       
452        Variable nAdd
453       
454        for(ii=0;ii<num;ii+=1)
455                dirVec[0] = PickDirection()
456                dirVec[1] = PickDirection()
457                dirVec[2] = PickDirection()
458                stVec[0] = x3d[ii]
459                stVec[1] = y3d[ii]
460                stVec[2] = z3d[ii]
461                nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
462                //FillSphereRadius(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],1)
463        endfor
464       
465        return(0)
466End
467
468///////////////
469               
470Function X_CylindersAtPoints(mat,num,rad,len,sobol,periodic)
471        Wave mat
472        variable num,rad,len            //length in direction
473        Variable periodic               //==1 if periodic boundaries
474        Variable sobol          //==1 if sobol selection of points (2D)
475       
476
477        Variable np
478        np = DimSize(mat,0)                     // assumes that all dimensions are the same
479       
480        // fill a 2D plane with points
481        Make/O/B/N=(np,np) plane
482        plane = 0
483        if(sobol)
484                if(WaveExists($"SobolInit") == 0)
485                        Make/O/D/N=2 SobolInit
486                        SobolX(-1,SobolInit)
487                else
488                        SobolPoints2D(plane,num)
489                endif
490        else
491                RandomPoints2D(plane,num)
492        endif
493        // put it in the proper plane of the matrix
494        mat[np/2][][] = plane[q][r]                     // in the YZ plane
495       
496        ParseMatrix3D_rho(mat)
497        Wave x3d=x3d
498        Wave y3d=y3d
499        Wave z3d=z3d
500       
501        Variable ii=0
502        NVAR  grid=root:FFT_T
503       
504        for(ii=0;ii<num;ii+=1)
505                FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1)               //cylinder 1
506        endfor
507       
508        return(0)
509End
510
511
512Function X_CylindersSquareGrid(mat,rad,len,sep)
513        Wave mat
514        variable rad,len                //length of cylinders
515        Variable sep                    // EDGE separation, in same units as cylinder
516       
517       
518        NVAR grid=root:FFT_T
519        Variable np,spacing
520        np = DimSize(mat,0)                     // assumes that all dimensions are the same
521       
522        // fill a 2D plane with points
523        Make/O/B/N=(np,np) plane
524        plane = 0
525
526
527        spacing = round(sep/grid)               // so it's an integer
528        FillPlaneSquareGrid(plane,spacing)
529       
530        // put it in the proper plane of the matrix
531        mat[np/2][][] = plane[q][r]                     // in the YZ plane
532       
533        ParseMatrix3D_rho(mat)
534        Wave x3d=x3d
535        Wave y3d=y3d
536        Wave z3d=z3d
537       
538        Variable ii=0,num
539        num = numpnts(x3d)     
540       
541        for(ii=0;ii<num;ii+=1)
542                FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1)               //cylinder 1
543        endfor
544       
545        return(0)
546End
547
548Function X_CylindersHexagonalGrid(mat,rad,len,sep)
549        Wave mat
550        variable rad,len                //length of cylinders
551        Variable sep                    // EDGE separation, in same units as cylinder
552       
553       
554        NVAR grid=root:FFT_T
555        Variable np,spacing
556        np = DimSize(mat,0)                     // assumes that all dimensions are the same
557       
558        // fill a 2D plane with points
559        Make/O/B/N=(np,np) plane
560        plane = 0
561
562        spacing = round(sep/grid)               // so it's an integer
563        FillPlaneHexagonal(plane,spacing)
564       
565        // put it in the proper plane of the matrix
566        mat[np/2][][] = plane[q][r]                     // in the YZ plane
567       
568        ParseMatrix3D_rho(mat)
569        Wave x3d=x3d
570        Wave y3d=y3d
571        Wave z3d=z3d
572       
573        Variable ii=0,num
574        num = numpnts(x3d)     
575       
576        for(ii=0;ii<num;ii+=1)
577                FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1)               //cylinder 1
578        endfor
579
580//// makes a crude core-shell cylinder 
581//      for(ii=0;ii<num;ii+=1)
582//              FillXCylinder(mat,grid,rad-10,x3d[ii],y3d[ii],z3d[ii],len,3)            //cylinder 1
583//      endfor
584//     
585        return(0)
586End
587
588
589
590Function FillPlaneSquareGrid(w,sp)
591        Wave w
592        Variable sp             //spacing in pixels
593       
594        Variable np,ii,jj
595       
596        np = DimSize(w,0)
597       
598        for(ii=sp;ii<(np-sp/2);ii+=sp)
599                for(jj=sp;jj<(np-sp/2);jj+=sp)
600                        w[ii][jj] = 1
601                endfor
602        endfor
603       
604       
605        return(0)
606End
607
608// sp is the
609Function FillPlaneHexagonal(w,sp)
610        Wave w
611        Variable sp             //spacing in pixels
612       
613        Variable np,ii,jj,xx,yy,fill
614       
615        np = DimSize(w,0)
616       
617        fill = 1
618        // do even columns, then odds
619        for(ii=sp;ii<(np-sp/2);ii+=(sp*sqrt(3)) )
620                for(jj=sp;jj<(np-sp/2);jj+=sp)
621                        w[ii][jj] = fill
622                endfor
623        endfor
624       
625        // now odds
626        for(ii=(sp+sp*sqrt(3)/2);ii<(np-sp/2);ii+=(sp*sqrt(3)))
627                for(jj=sp/2;jj<(np-sp/2);jj+=sp)
628                        w[ii][jj] = fill
629                endfor
630        endfor
631       
632        return(0)
633End
634
635Proc HexagonalCylinders()
636
637        Variable xc,yc,zc,len,rad,aa,cosa,sina,meanRad,pd,meanLen
638       
639        Variable grid=root:FFT_T
640        Variable npt=root:FFT_N
641       
642        xc=npt/2
643        yc=npt/2
644        zc=npt/2
645       
646        meanRad=20
647        pd=0.01
648        meanLen=400
649        len=meanLen
650       
651//      aa=meanRad/grid*2 +40   //grid units, +2 gives 2 units of separation
652        aa=meanRad/grid*2 +2    //grid units, +2 gives 2 units of separation
653        cosa = round(0.5*aa)
654        sina = round(0.866*aa)
655       
656
657        Variable ii,nPass,x1,x2,x3
658        ii=0
659        nPass=1
660
661
662        x1 = Gauss(meanRad-pd*meanRad,meanRad,pd*meanRad)
663        x2 = Gauss(meanRad,meanRad,pd*meanRad)
664        x3 = Gauss(meanRad+pd*meanRad,meanRad,pd*meanRad)
665       
666        Print x1,x2,x3
667//     
668/////////// average over 3 point distribution --
669//      //calculate and sum the three contributions
670//      FastOp mat=0
671//      rad = meanRad-pd*meanRad
672////    len = meanLen + gnoise(pd*meanLen)
673//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
674//      Execute "DoFFT()"
675//      Duplicate/O iBin iBin_1g_128
676//      iBin_1g_128 *= x1
677//     
678//      FastOp mat=0
679//      rad = meanRad
680////    len = meanLen + gnoise(pd*meanLen)
681//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
682//      Execute "DoFFT()"
683//      iBin_1g_128 += x2*iBin
684//     
685//      FastOp mat=0
686//      rad = meanRad+pd*meanRad
687////    len = meanLen + gnoise(pd*meanLen)
688//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
689//      Execute "DoFFT()"
690//      iBin_1g_128 += x3*iBin
691//     
692//      iBin_1g_128 /= (x1+x2+x3)               //renormalize
693//     
694////// end 3 pt average
695
696//////////////////////
697// paired cylinders
698/////////////////////////
699//      //calculate and sum the three contributions
700//      FastOp mat=0
701//      rad = meanRad-pd*meanRad
702////    len = meanLen + gnoise(pd*meanLen)
703//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
704//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
705//      Execute "DoFFT()"
706//      Duplicate/O iBin iBin_2g_256_sep40
707//      iBin_2g_256_sep40 *= x1
708//     
709//      FastOp mat=0
710//      rad = meanRad
711////    len = meanLen + gnoise(pd*meanLen)
712//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
713//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
714//      Execute "DoFFT()"
715//      iBin_2g_256_sep40 += x2*iBin
716//     
717//      FastOp mat=0
718//      rad = meanRad+pd*meanRad
719////    len = meanLen + gnoise(pd*meanLen)
720//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
721//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
722//      Execute "DoFFT()"
723//      iBin_2g_256_sep40 += x3*iBin
724//     
725//      iBin_2g_256_sep40 /= (x1+x2+x3)         //renormalize
726       
727/////////////////////////       
728//
729//
730////
731//      ii = 0 
732//      do
733//              FastOp mat=0
734//              rad = meanRad + gnoise(pd*meanRad)
735//              print ii,rad
736//              FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
737//              Execute "DoFFT()"
738//              if(ii==0)       
739//                      Duplicate/O iBin iBin_1g_256
740//              else
741//                      iBin_1g_256 += iBin
742//              endif
743//              ii+=1
744//      while(ii<nPass)
745//      iBin_1g_256 /= nPass
746//     
747/////   
748//      ii=0
749//      nPass=1
750//      do
751//              FastOp mat=0
752//              rad = meanRad + gnoise(pd*meanRad)
753//              print ii,rad
754//              FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)              //cylinder 1
755//              rad = meanRad + gnoise(pd*meanRad)
756//              print ii,rad
757//              FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
758//              Execute "DoFFT()"
759//              if(ii==0)       
760//                      Duplicate/O iBin iBin_2g_256
761//              else
762//                      iBin_2g_256 += iBin
763//              endif
764//              ii+=1
765//      while(ii<nPass)
766//      iBin_2g_256 /= nPass
767//     
768        ii=0
769        npass=1
770        do
771                FastOp mat=0
772                rad = meanRad + gnoise(pd*meanRad)
773                print ii,rad
774                FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)              //cylinder 1
775                rad = meanRad + gnoise(pd*meanRad)
776                print ii,rad
777                FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)           //cylinder 2
778                rad = meanRad + gnoise(pd*meanRad)
779                print ii,rad
780                FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,1)            //cylinder 3
781                Execute "DoFFT()"
782                if(ii==0)       
783                        Duplicate/O iBin iBin_3g_256
784                else
785                        iBin_3g_256 += iBin
786                endif
787                ii+=1
788        while(ii<nPass)
789        iBin_3g_256 /= nPass
790//     
791        rad = meanRad + gnoise(pd*meanRad)
792        print rad
793        FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)
794        Execute "DoFFT()"
795        Duplicate/O iBin iBin_2
796///     
797        rad = meanRad + gnoise(pd*meanRad)
798        print rad
799        FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,1)
800        Execute "DoFFT()"
801        Duplicate/O iBin iBin_3
802///     
803        rad = meanRad + gnoise(pd*meanRad)
804        print rad
805        FillZCylinder(mat,grid,rad,xc+sina,yc-cosa,zc,len,1)
806        Execute "DoFFT()"
807        Duplicate/O iBin iBin_4
808///     
809        rad = meanRad + gnoise(pd*meanRad)
810        print rad
811        FillZCylinder(mat,grid,rad,xc,yc-aa,zc,len,1)
812        Execute "DoFFT()"
813        Duplicate/O iBin iBin_5
814///     
815        rad = meanRad + gnoise(pd*meanRad)
816        print rad
817        FillZCylinder(mat,grid,rad,xc-sina,yc-cosa,zc,len,1)
818        Execute "DoFFT()"
819        Duplicate/O iBin iBin_6
820///     
821        rad = meanRad + gnoise(pd*meanRad)
822        print rad
823        FillZCylinder(mat,grid,rad,xc-sina,yc+cosa,zc,len,1)
824        Execute "DoFFT()"
825        Duplicate/O iBin iBin_7
826
827end
828
829
830Proc Vary_N_In_Direction()
831
832        FFTEraseMatrixButtonProc("")
833        ConnectedRodFill(mat,2,30000*4,1)
834        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
835        Duplicate/O ival_XOP ival_XOP_2
836        Duplicate/O qval_XOP qval_XOP_2
837       
838        FFTEraseMatrixButtonProc("")
839        ConnectedRodFill(mat,3,30000*4,1)
840        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
841        Duplicate/O ival_XOP ival_XOP_3
842        Duplicate/O qval_XOP qval_XOP_3
843       
844        FFTEraseMatrixButtonProc("")
845        ConnectedRodFill(mat,4,30000*4,1)
846        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
847        Duplicate/O ival_XOP ival_XOP_4
848        Duplicate/O qval_XOP qval_XOP_4
849       
850        FFTEraseMatrixButtonProc("")
851        ConnectedRodFill(mat,5,30000*4,1)
852        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
853        Duplicate/O ival_XOP ival_XOP_5
854        Duplicate/O qval_XOP qval_XOP_5
855       
856        FFTEraseMatrixButtonProc("")
857        ConnectedRodFill(mat,6,30000*4,1)
858        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
859        Duplicate/O ival_XOP ival_XOP_6
860        Duplicate/O qval_XOP qval_XOP_6
861       
862        FFTEraseMatrixButtonProc("")
863        ConnectedRodFill(mat,7,30000*4,1)
864        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
865        Duplicate/O ival_XOP ival_XOP_7
866        Duplicate/O qval_XOP qval_XOP_7
867
868        FFTEraseMatrixButtonProc("")
869        ConnectedRodFill(mat,8,30000*4,1)
870        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
871        Duplicate/O ival_XOP ival_XOP_8
872        Duplicate/O qval_XOP qval_XOP_8
873       
874        FFTEraseMatrixButtonProc("")
875        ConnectedRodFill(mat,9,30000*4,1)
876        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
877        Duplicate/O ival_XOP ival_XOP_9
878        Duplicate/O qval_XOP qval_XOP_9
879       
880        FFTEraseMatrixButtonProc("")
881        ConnectedRodFill(mat,10,30000*4,1)
882        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
883        Duplicate/O ival_XOP ival_XOP_10
884        Duplicate/O qval_XOP qval_XOP_10
885               
886End
Note: See TracBrowser for help on using the repository browser.