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

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

fixing the shape filling functions to explicitly take fill as a parameter, and make the default fill value = 10 == 10 e-7 A-2

More additions to the RealSpeca? help file

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