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

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

Changes to the real-space modeling to sped up the drawing of cylinders, and to provide a few more examples of scritped, unattended calculations, and saving the 3D matrix.
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                Print "cyl = ",ii,num
619        endfor
620
621//// makes a crude core-shell cylinder 
622//      for(ii=0;ii<num;ii+=1)
623//              FillXCylinder(mat,grid,rad-10,x3d[ii],y3d[ii],z3d[ii],len,3*fill)               //cylinder 1
624//      endfor
625//     
626        return(0)
627End
628
629
630
631Function FillPlaneSquareGrid(w,sp,fill)
632        Wave w
633        Variable sp             //spacing in pixels
634        Variable fill
635       
636        Variable np,ii,jj
637       
638        np = DimSize(w,0)
639       
640        for(ii=sp;ii<(np-sp/2);ii+=sp)
641                for(jj=sp;jj<(np-sp/2);jj+=sp)
642                        w[ii][jj] = fill
643                endfor
644        endfor
645       
646       
647        return(0)
648End
649
650// sp is the
651Function FillPlaneHexagonal(w,sp,fill)
652        Wave w
653        Variable sp             //spacing in pixels
654        Variable fill
655       
656        Variable np,ii,jj,xx,yy
657       
658        np = DimSize(w,0)
659       
660
661        // do even columns, then odds
662        for(ii=sp;ii<(np-sp/2);ii+=(sp*sqrt(3)) )
663                for(jj=sp;jj<(np-sp/2);jj+=sp)
664                        w[ii][jj] = fill
665                endfor
666        endfor
667       
668        // now odds
669        for(ii=(sp+sp*sqrt(3)/2);ii<(np-sp/2);ii+=(sp*sqrt(3)))
670                for(jj=sp/2;jj<(np-sp/2);jj+=sp)
671                        w[ii][jj] = fill
672                endfor
673        endfor
674       
675        return(0)
676End
677
678Proc HexagonalCylinders()
679
680        Variable xc,yc,zc,len,rad,aa,cosa,sina,meanRad,pd,meanLen
681       
682        Variable grid=root:FFT_T
683        Variable npt=root:FFT_N
684       
685        xc=npt/2
686        yc=npt/2
687        zc=npt/2
688       
689        meanRad=20
690        pd=0.01
691        meanLen=400
692        len=meanLen
693       
694//      aa=meanRad/grid*2 +40   //grid units, +2 gives 2 units of separation
695        aa=meanRad/grid*2 +2    //grid units, +2 gives 2 units of separation
696        cosa = round(0.5*aa)
697        sina = round(0.866*aa)
698       
699
700        Variable ii,nPass,x1,x2,x3
701        ii=0
702        nPass=1
703
704
705        x1 = Gauss(meanRad-pd*meanRad,meanRad,pd*meanRad)
706        x2 = Gauss(meanRad,meanRad,pd*meanRad)
707        x3 = Gauss(meanRad+pd*meanRad,meanRad,pd*meanRad)
708       
709        Print x1,x2,x3
710//     
711/////////// average over 3 point distribution --
712//      //calculate and sum the three contributions
713//      FastOp mat=0
714//      rad = meanRad-pd*meanRad
715////    len = meanLen + gnoise(pd*meanLen)
716//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
717//      Execute "DoFFT()"
718//      Duplicate/O iBin iBin_1g_128
719//      iBin_1g_128 *= x1
720//     
721//      FastOp mat=0
722//      rad = meanRad
723////    len = meanLen + gnoise(pd*meanLen)
724//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
725//      Execute "DoFFT()"
726//      iBin_1g_128 += x2*iBin
727//     
728//      FastOp mat=0
729//      rad = meanRad+pd*meanRad
730////    len = meanLen + gnoise(pd*meanLen)
731//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
732//      Execute "DoFFT()"
733//      iBin_1g_128 += x3*iBin
734//     
735//      iBin_1g_128 /= (x1+x2+x3)               //renormalize
736//     
737////// end 3 pt average
738
739//////////////////////
740// paired cylinders
741/////////////////////////
742//      //calculate and sum the three contributions
743//      FastOp mat=0
744//      rad = meanRad-pd*meanRad
745////    len = meanLen + gnoise(pd*meanLen)
746//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
747//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
748//      Execute "DoFFT()"
749//      Duplicate/O iBin iBin_2g_256_sep40
750//      iBin_2g_256_sep40 *= x1
751//     
752//      FastOp mat=0
753//      rad = meanRad
754////    len = meanLen + gnoise(pd*meanLen)
755//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
756//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
757//      Execute "DoFFT()"
758//      iBin_2g_256_sep40 += x2*iBin
759//     
760//      FastOp mat=0
761//      rad = meanRad+pd*meanRad
762////    len = meanLen + gnoise(pd*meanLen)
763//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
764//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
765//      Execute "DoFFT()"
766//      iBin_2g_256_sep40 += x3*iBin
767//     
768//      iBin_2g_256_sep40 /= (x1+x2+x3)         //renormalize
769       
770/////////////////////////       
771//
772//
773////
774//      ii = 0 
775//      do
776//              FastOp mat=0
777//              rad = meanRad + gnoise(pd*meanRad)
778//              print ii,rad
779//              FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
780//              Execute "DoFFT()"
781//              if(ii==0)       
782//                      Duplicate/O iBin iBin_1g_256
783//              else
784//                      iBin_1g_256 += iBin
785//              endif
786//              ii+=1
787//      while(ii<nPass)
788//      iBin_1g_256 /= nPass
789//     
790/////   
791//      ii=0
792//      nPass=1
793//      do
794//              FastOp mat=0
795//              rad = meanRad + gnoise(pd*meanRad)
796//              print ii,rad
797//              FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)              //cylinder 1
798//              rad = meanRad + gnoise(pd*meanRad)
799//              print ii,rad
800//              FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
801//              Execute "DoFFT()"
802//              if(ii==0)       
803//                      Duplicate/O iBin iBin_2g_256
804//              else
805//                      iBin_2g_256 += iBin
806//              endif
807//              ii+=1
808//      while(ii<nPass)
809//      iBin_2g_256 /= nPass
810//     
811        ii=0
812        npass=1
813        do
814                FastOp mat=0
815                rad = meanRad + gnoise(pd*meanRad)
816                print ii,rad
817                FillZCylinder(mat,grid,rad,xc,yc,zc,len,10)             //cylinder 1
818                rad = meanRad + gnoise(pd*meanRad)
819                print ii,rad
820                FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,10)          //cylinder 2
821                rad = meanRad + gnoise(pd*meanRad)
822                print ii,rad
823                FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,10)           //cylinder 3
824                Execute "DoFFT()"
825                if(ii==0)       
826                        Duplicate/O iBin iBin_3g_256
827                else
828                        iBin_3g_256 += iBin
829                endif
830                ii+=1
831        while(ii<nPass)
832        iBin_3g_256 /= nPass
833//     
834        rad = meanRad + gnoise(pd*meanRad)
835        print rad
836        FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,10)
837        Execute "DoFFT()"
838        Duplicate/O iBin iBin_2
839///     
840        rad = meanRad + gnoise(pd*meanRad)
841        print rad
842        FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,10)
843        Execute "DoFFT()"
844        Duplicate/O iBin iBin_3
845///     
846        rad = meanRad + gnoise(pd*meanRad)
847        print rad
848        FillZCylinder(mat,grid,rad,xc+sina,yc-cosa,zc,len,10)
849        Execute "DoFFT()"
850        Duplicate/O iBin iBin_4
851///     
852        rad = meanRad + gnoise(pd*meanRad)
853        print rad
854        FillZCylinder(mat,grid,rad,xc,yc-aa,zc,len,10)
855        Execute "DoFFT()"
856        Duplicate/O iBin iBin_5
857///     
858        rad = meanRad + gnoise(pd*meanRad)
859        print rad
860        FillZCylinder(mat,grid,rad,xc-sina,yc-cosa,zc,len,10)
861        Execute "DoFFT()"
862        Duplicate/O iBin iBin_6
863///     
864        rad = meanRad + gnoise(pd*meanRad)
865        print rad
866        FillZCylinder(mat,grid,rad,xc-sina,yc+cosa,zc,len,10)
867        Execute "DoFFT()"
868        Duplicate/O iBin iBin_7
869
870end
871
872
873Proc Vary_N_In_Direction()
874
875        Variable fill=10
876       
877        FFTEraseMatrixButtonProc("")
878        ConnectedRodFill(mat,2,30000*4,1,fill)
879        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
880        Duplicate/O ival_XOP ival_XOP_2
881        Duplicate/O qval_XOP qval_XOP_2
882       
883        FFTEraseMatrixButtonProc("")
884        ConnectedRodFill(mat,3,30000*4,1,fill)
885        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
886        Duplicate/O ival_XOP ival_XOP_3
887        Duplicate/O qval_XOP qval_XOP_3
888       
889        FFTEraseMatrixButtonProc("")
890        ConnectedRodFill(mat,4,30000*4,1,fill)
891        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
892        Duplicate/O ival_XOP ival_XOP_4
893        Duplicate/O qval_XOP qval_XOP_4
894       
895        FFTEraseMatrixButtonProc("")
896        ConnectedRodFill(mat,5,30000*4,1,fill)
897        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
898        Duplicate/O ival_XOP ival_XOP_5
899        Duplicate/O qval_XOP qval_XOP_5
900       
901        FFTEraseMatrixButtonProc("")
902        ConnectedRodFill(mat,6,30000*4,1,fill)
903        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
904        Duplicate/O ival_XOP ival_XOP_6
905        Duplicate/O qval_XOP qval_XOP_6
906       
907        FFTEraseMatrixButtonProc("")
908        ConnectedRodFill(mat,7,30000*4,1,fill)
909        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
910        Duplicate/O ival_XOP ival_XOP_7
911        Duplicate/O qval_XOP qval_XOP_7
912
913        FFTEraseMatrixButtonProc("")
914        ConnectedRodFill(mat,8,30000*4,1,fill)
915        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
916        Duplicate/O ival_XOP ival_XOP_8
917        Duplicate/O qval_XOP qval_XOP_8
918       
919        FFTEraseMatrixButtonProc("")
920        ConnectedRodFill(mat,9,30000*4,1,fill)
921        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
922        Duplicate/O ival_XOP ival_XOP_9
923        Duplicate/O qval_XOP qval_XOP_9
924       
925        FFTEraseMatrixButtonProc("")
926        ConnectedRodFill(mat,10,30000*4,1,fill)
927        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
928        Duplicate/O ival_XOP ival_XOP_10
929        Duplicate/O qval_XOP qval_XOP_10
930               
931End
Note: See TracBrowser for help on using the repository browser.