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

#pragma rtGlobals=1             // Use modern global access method.

// VERY incomplete methods for filling volumes with connected cylinders
//
// most is still being worked out and verified in 2D
//



//These macros may or may not work - designed to help automate the 
//process of calculating a series of configurations.



Proc ConnectDots3D(maxNumConn)
        Variable maxNumConn=1
        
        Variable num=numpnts(x3d)
        Make/O/N=(num) numConnection3D=0
        Make/O/T/N=(num) connectedTo3D=""
        fConnectDots3D(maxNumConn)
end

// connect the dots, no periodic boundary conditions
//thickness is set to 1
Function fConnectDots3D(maxNumConn)
        Variable maxNumConn
        
        
        WAVE x3d=x3d
        WAVE y3d=y3d
        WAVE z3d=z3d
        WAVE numConnection3D = numConnection3D
        WAVE/T connectedTo3D = connectedTo3D
        WAVE mat=mat
        Variable num=numpnts(x3d),ii=0
        Variable matSize=DimSize(mat,0),jj

        Variable nnX,nnY,nnZ,nnInd,nnDist               //return values of the nearest neighbor
        Variable thick = 1
        do
                for(ii=0;ii<num;ii+=1)
                        if(numConnection3D[ii] < maxNumConn)
                                FindNearestNeighbor3D(x3d,y3d,z3d,ii,maxNumConn,nnX,nnY,nnZ,nnInd,nnDist)
                                
                                numConnection3D[nnInd] += 1
                                numConnection3D[ii] += 1
                                connectedTo3D[nnInd] += num2str(ii)+","
                                connectedTo3D[ii] += num2str(nnInd)+","
                                
                                ConnectPoints3D(mat, x3d[ii],y3d[ii],z3d[ii],x3d[nnInd],y3d[nnInd],z3d[nnInd],thick,1)          //always fill
                        endif
                endfor
                //without rotation, a very closed network structure tends to form (bias?)
                // rotate (10?) points so that we're not starting from the same place every time
                // -- but the "connectedTo" information is LOST...
                //Rotate 10, x2d,y2d, numConnection,connectedTo
                //Print sum(numConnection,-inf,inf), num*maxNumConn
        while(sum(numConnection3D,-inf,inf) < num*maxNumConn)

End

Function FindNearestNeighbor3D(xw,yw,zw,startPt,maxConn,nnX,nnY,nnZ,nnInd,nnDist)
        Wave xw,yw,zw
        Variable startPt,maxConn,&nnX,&nnY,&nnZ,&nnInd,&nnDist
        
        Variable num,matSize,numToFind
        
        //make waves to hold the answers
        num=numpnts(xw)
        Make/O/N=(num) distWave3D,distIndex3D
        
        //now calculate the distances
        Variable ii,dist,testPt
        
        //this is probably the slowest step by far in the process...(tiny XOP?)
        for(ii=0;ii<num;ii+=1)
                distWave3D[ii] = (xw[ii]-xw[startPt])^2 +  (yw[ii]-yw[startPt])^2        + (zw[ii]-zw[startPt])^2               //d^2
        endfor

        MakeIndex distWave3D distIndex3D                //then distWave[distIndex[ii]] will loop through in order of distance
        WAVE numConnection3D = numConnection3D
        WAVE/T connectedTo3D = connectedTo3D
        WAVE mat = mat
        
        for(ii=1;ii<num;ii+=1)          //[0] point is the test point, dist == 0
                testPt =  distIndex3D[ii]               //index
                if(numConnection3D[testPt] < maxConn )          //can test pt accept another connection?
                        if(WhichListItem(num2str(testPt),connectedTo3D[startPt],",",0) == -1) // not already connected to the starting point?
        //                      Print ii,testPt
        //                      Printf "nearest point is (%d,%d), dist^2 = %g\r",xw[testPt],yw[testPt],distWave[testPt]
                                nnX = xw[testPt]                //return values as pass-by-reference
                                nnY = yw[testPt]
                                nnZ = zw[testPt]
                                nnInd = testPt
                                nnDist = distWave3D[testPt]

                                return(0)               //found a point, return
                        endif
                endif
        endfor

        return (1)              //did not find a point, return error
End

///!!! fails if the points are along an axis (or even close) - then the box becomes infinitely thin
Function ConnectPoints3D_old(w, x1,y1,z1,x2,y2,z2,thick,fill)
        Wave w
        Variable  x1,y1,z1,x2,y2,z2,thick,fill
        
        //find all of the points within dist=thick of the line defined by p1 and p2
        // and give these points (in w) the value 1
        
        // search through the box defined by p1 and p2
        Variable ii,jj,kk,dist
        
        for(ii=min(x1,x2);ii<=max(x1,x2);ii+=1)
                for(jj=min(y1,y2);jj<=max(y1,y2);jj+=1)
                        for(kk=min(z1,z2);kk<=max(z1,z2);kk+=1)
                                //find distance between current point and P1P2
                                dist = Pt_Line_Distance3D(ii,jj,kk, x1,y1,z1,x2,y2,z2)
                                //print dist
                                if(dist<=thick)
                                        w[ii][jj][kk] = fill            //add the point
                                endif
                        endfor
                endfor
        endfor
        
End

//
// thickness is given in voxels
//
// the search box is expanded by the thickness in each direction (within the box)
// to account for connections that are along an axis direction
// 
Function ConnectPoints3D(w, x1,y1,z1,x2,y2,z2,thick,fill)
        Wave w
        Variable  x1,y1,z1,x2,y2,z2,thick,fill
        
        //find all of the points within dist=thick of the line defined by p1 and p2
        // and give these points (in w) the value 1
        
        // search through the box defined by p1 and p2
        Variable ii,jj,kk,dist,minX,minY,minZ,maxX,maxY,maxZ
        Variable boxMax,slop
        boxMax = DimSize(w,0) - 1               //max index of box
        slop = thick + 0                // to make the bounding box large enough the get it all
        
        minX=min(x1,x2)                 //give me the smaller x
        minX=max(minX-slop,0)   //give me zero if it's negative
        minY=min(y1,y2)
        minY=max(minY-slop,0)
        minZ=min(z1,z2)
        minZ=max(minZ-slop,0)
        
        maxX=max(x1,x2)                                 //give me the larger x
        maxX=min(maxX+slop,boxMax)              //give the the box dimension if it's too big
        maxY=max(y1,y2)
        maxY=min(maxY+slop,boxMax)
        maxZ=max(z1,z2)
        maxZ=min(maxZ+slop,boxMax)
                
        for(ii=minX;ii<=maxX;ii+=1)
                for(jj=minY;jj<=maxY;jj+=1)
                        for(kk=minZ;kk<=maxZ;kk+=1)
                                //find distance between current point and P1P2
                                dist = Pt_Line_Distance3D(ii,jj,kk, x1,y1,z1,x2,y2,z2)
                                //print dist
                                if(dist<=thick)
                                        w[ii][jj][kk] = fill            //add the point
                                endif
                        endfor
                endfor
        endfor
        
End

//solution courtesy of the math forum @ Drexel "Ask Dr. Math"
Function Pt_Line_Distance3D(xx,yy,zz, x1,y1,z1,x2,y2,z2)
        Variable xx,yy,zz, x1,y1,z1,x2,y2,z2
        
        Variable dist,x3,y3,z3,x4,y4,z4,aa,ab
        
        //AB
        x3 = x2-x1
        y3 = y2-y1
        z3 = z2-z1
        //AP
        x4 = xx-x1
        y4 = yy-y1
        z4 = zz-z1
        
        //aa = length of ABxAP
        aa = sqrt( (y3*z4-z3*y4)^2 + (x3*z4-x4*z3)^2 + (x3*y4-x4*y3)^2 )
        
        //ab = length of AB
        ab = sqrt(x3^2+y3^2+z3^2)
        
        dist = aa/ab
        
        return(dist)
End

//////////////////////////
//// random rod fills - just a voxel at each point, so use T=diameter

// fill a volume with "rods" that are made up of a numer of spheres
// lined up in a random direction
//
// all of the spheres/rods that are added must be connected = "adjacent" to
// an occupied voxel.
//
// this implementation is a random walk, end-to-end connection
//
Function ConnectedRodFill(mat,len,num,periodic)
        Wave mat
        variable len,num                //length in direction, number of spheres to add
        Variable periodic               //==1 if periodic boundaries
        
        Variable nptx,npty,nptz,ii,fail=0,nAdd
        Make/O/I/N=3 stVec,dirVec,endVec
        
        nptx=DimSize(mat,0)
        npty=DimSize(mat,1)
        nptz=DimSize(mat,2)
        
        //place first rod randomly
        stVec[0]=trunc(abs(enoise(nptx)))               //random integer distr betw (0,npt)
        stVec[1]=trunc(abs(enoise(npty)))
        stVec[2]=trunc(abs(enoise(nptz)))
        dirVec[0] = PickDirection() 
        dirVec[1] = PickDirection() 
        dirVec[2] = PickDirection() 

        nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
        
        stVec = endVec          //end point is the new starting point

        //loop over remaining
        ii=nAdd
        do

                        dirVec[0] = PickDirection() 
                        dirVec[1] = PickDirection() 
                        dirVec[2] = PickDirection() 
                        nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
                        stVec = endVec          //end point is the new starting point

                        ii+=nAdd

                        if(mod(ii,100000) == 0)
                                print "Point = ",ii
                        endif
                        
        while(ii<num)
        
        
        return(0)
End

// fill a volume with "rods" that are made up of a numer of spheres
// lined up in a random direction
//
//
Function UnConnectedRodFill(mat,len,num,periodic)
        Wave mat
        variable len,num                //length in direction, number of spheres to add
        Variable periodic               //==1 if periodic boundaries
        
        Variable nptx,npty,nptz,ii,fail=0,nAdd
        Make/O/I/N=3 stVec,dirVec,endVec

        nptx=DimSize(mat,0)
        npty=DimSize(mat,1)
        nptz=DimSize(mat,2)
        
        //place first rod randomly
        stVec[0]=trunc(abs(enoise(nptx)))               //random integer distr betw (0,npt)
        stVec[1]=trunc(abs(enoise(npty)))
        stVec[2]=trunc(abs(enoise(nptz)))
        dirVec[0] = PickDirection() 
        dirVec[1] = PickDirection() 
        dirVec[2] = PickDirection() 

        nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
        //loop over remaining
        ii=nAdd
        do
                //pick a new point
                stVec[0]=trunc(abs(enoise(nptx)))               //random integer distr betw (0,npt)
                stVec[1]=trunc(abs(enoise(npty)))
                stVec[2]=trunc(abs(enoise(nptz)))
                dirVec[0] = PickDirection() 
                dirVec[1] = PickDirection() 
                dirVec[2] = PickDirection() 
                //add as many spheres as possible in that direction
                nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
                ii+=nAdd
                //Print "point ",ii

        while(ii<num)
        Print "failures = ",fail
        
//      ParseMatrix3D_rho(mat)          //convert to XYZ
        
        return(0)
End


Function PickDirection()
        return( trunc(abs(enoise(3))) - 1       )       //integer -1,0,1
End

//need to add periodic boundary conditions (done)
//
// need to insist on connectivity
//
//
Function AddSpheresInDirection(m,stVec,dirVec,endVec,nTry,periodic)
        WAVE m,stVec,dirVec,endVec
        Variable nTry,periodic
        
        Variable added,nx,ny,nz,nDimx,nDimy,nDimz
        
        nDimx=DimSize(m,0)-1
        nDimy=DimSize(m,1)-1
        nDimz=DimSize(m,2)-1
        added=0
        nx=stVec[0]
        ny=stVec[1]
        nz=stVec[2]
        
        do
                nx+=dirVec[0]           //coordinates on next test point to try to add in d(xyz) direction
                ny+=dirVec[1]
                nz+=dirVec[2]
                //check for out-of bounds - if periodic ==1, reflect to opposite face of cube, else return
                if( (nx>nDimx) || (nx<0) )
                        if(periodic)
                                nx = abs(nDimx+1 - abs(nx))
                        else
                                break
                        endif
                Endif
                if( (ny>nDimy) || (ny<0) )
                        if(periodic)
                                ny = abs(nDimy+1 - abs(ny))
                        else
                                break
                        endif
                Endif
                if( (nz>nDimz) || (nz<0) )
                        if(periodic)
                                nz = abs(nDimz+1 - abs(nz))
                        else
                                break
                        endif
                Endif
                //add the point or exit
                if(m[nx][ny][nz])
                        break           //point is already occupied
                else
                        m[nx][ny][nz] = 1
//                      nx +=1          //we'll back this off before returning
//                      ny +=1
//                      nz +=1
                        added += 1
                endif
        while(added<nTry)
        
        endVec[0] = nx - dirVec[0]                      //back up to the last point actually added
        endVec[1] = ny - dirVec[1]      
        endVec[2] = nz - dirVec[2]      
        return added
End

Proc PutRandomCylindersAtPoints(w,num,len,periodic)
        String w="mat"
        Variable num=1000,len=10,periodic=1
        Prompt w,"matrix"
        Prompt num,"number of starting points"
        prompt len,"length of cylinders"
        prompt periodic"1=periodic, 0=non-periodic boundaries"
        
        $w=0
        RandomFill3DMat($w,num)
        CylindersAtPoints($w,len,periodic)
        
        NumberOfPoints()
end

Proc PutXAxisCylindersAtPoints(w,num,rad,len,periodic,sobol)
        String w="mat"
        Variable num=100,rad=20,len=300,periodic=1,sobol=1
        Prompt w,"matrix"
        Prompt num,"number of starting points"
        prompt rad,"radius of cylinders"
        prompt len,"length of cylinders"
        prompt periodic,"1=periodic, 0=non-periodic boundaries"
        Prompt sobol,"1=Sobol, 0=random"
        
        $w=0
        X_CylindersAtPoints($w,num,rad,len,sobol,periodic)
        
        NumberOfPoints()
end

Proc PutXAxisCylindersSquare(w,rad,len,sep)
        String w="mat"
        Variable rad=20,len=300,sep=80
        Prompt w,"matrix"
        prompt rad,"radius of cylinders"
        prompt len,"length of cylinders"
        prompt sep,"center-to-center separation of cylinders"
        
        $w=0
        X_CylindersSquareGrid($w,rad,len,sep)
        
        NumberOfPoints()
end

Proc PutXAxisCylindersHexagonal(w,rad,len,sep)
        String w="mat"
        Variable rad=20,len=300,sep=80
        Prompt w,"matrix"
        prompt rad,"radius of cylinders"
        prompt len,"length of cylinders"
        prompt sep,"center-to-center separation of cylinders"
        
        $w=0
        X_CylindersHexagonalGrid($w,rad,len,sep)
        
        NumberOfPoints()
end

Function CylindersAtPoints(mat,len,periodic)
        Wave mat
        variable len            //length in direction
        Variable periodic               //==1 if periodic boundaries
        
        Make/O/B/N=3 stVec,dirVec,endVec

        ParseMatrix3D_rho(mat)
        Wave x3d=x3d
        Wave y3d=y3d
        Wave z3d=z3d
        
        Variable ii=0,num=numpnts(x3d)
        NVAR  grid=root:FFT_T
        
        Variable nAdd
        
        for(ii=0;ii<num;ii+=1)
                dirVec[0] = PickDirection() 
                dirVec[1] = PickDirection() 
                dirVec[2] = PickDirection() 
                stVec[0] = x3d[ii]
                stVec[1] = y3d[ii]
                stVec[2] = z3d[ii]
                nAdd = AddSpheresInDirection(mat,stVec,dirVec,endVec,len,periodic)
                //FillSphereRadius(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],1)
        endfor
        
        return(0)
End

///////////////
                
Function X_CylindersAtPoints(mat,num,rad,len,sobol,periodic)
        Wave mat
        variable num,rad,len            //length in direction
        Variable periodic               //==1 if periodic boundaries
        Variable sobol          //==1 if sobol selection of points (2D)
        

        Variable np
        np = DimSize(mat,0)                     // assumes that all dimensions are the same
        
        // fill a 2D plane with points
        Make/O/B/N=(np,np) plane
        plane = 0
        if(sobol)
                if(WaveExists($"SobolInit") == 0)
                        Make/O/D/N=2 SobolInit
                        SobolX(-1,SobolInit)
                else
                        SobolPoints2D(plane,num)
                endif
        else
                RandomPoints2D(plane,num)
        endif
        // put it in the proper plane of the matrix
        mat[np/2][][] = plane[q][r]                     // in the YZ plane
        
        ParseMatrix3D_rho(mat)
        Wave x3d=x3d
        Wave y3d=y3d
        Wave z3d=z3d
        
        Variable ii=0
        NVAR  grid=root:FFT_T
        
        for(ii=0;ii<num;ii+=1)
                FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1)               //cylinder 1
        endfor
        
        return(0)
End


Function X_CylindersSquareGrid(mat,rad,len,sep)
        Wave mat
        variable rad,len                //length of cylinders
        Variable sep                    // EDGE separation, in same units as cylinder
        
        
        NVAR grid=root:FFT_T
        Variable np,spacing
        np = DimSize(mat,0)                     // assumes that all dimensions are the same
        
        // fill a 2D plane with points
        Make/O/B/N=(np,np) plane
        plane = 0


        spacing = round(sep/grid)               // so it's an integer
        FillPlaneSquareGrid(plane,spacing)
        
        // put it in the proper plane of the matrix
        mat[np/2][][] = plane[q][r]                     // in the YZ plane
        
        ParseMatrix3D_rho(mat)
        Wave x3d=x3d
        Wave y3d=y3d
        Wave z3d=z3d
        
        Variable ii=0,num
        num = numpnts(x3d)      
        
        for(ii=0;ii<num;ii+=1)
                FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1)               //cylinder 1
        endfor
        
        return(0)
End

Function X_CylindersHexagonalGrid(mat,rad,len,sep)
        Wave mat
        variable rad,len                //length of cylinders
        Variable sep                    // EDGE separation, in same units as cylinder
        
        
        NVAR grid=root:FFT_T
        Variable np,spacing
        np = DimSize(mat,0)                     // assumes that all dimensions are the same
        
        // fill a 2D plane with points
        Make/O/B/N=(np,np) plane
        plane = 0

        spacing = round(sep/grid)               // so it's an integer
        FillPlaneHexagonal(plane,spacing)
        
        // put it in the proper plane of the matrix
        mat[np/2][][] = plane[q][r]                     // in the YZ plane
        
        ParseMatrix3D_rho(mat)
        Wave x3d=x3d
        Wave y3d=y3d
        Wave z3d=z3d
        
        Variable ii=0,num
        num = numpnts(x3d)      
        
        for(ii=0;ii<num;ii+=1)
                FillXCylinder(mat,grid,rad,x3d[ii],y3d[ii],z3d[ii],len,1)               //cylinder 1
        endfor

//// makes a crude core-shell cylinder  
//      for(ii=0;ii<num;ii+=1)
//              FillXCylinder(mat,grid,rad-10,x3d[ii],y3d[ii],z3d[ii],len,3)            //cylinder 1
//      endfor
//      
        return(0)
End



Function FillPlaneSquareGrid(w,sp)
        Wave w
        Variable sp             //spacing in pixels
        
        Variable np,ii,jj
        
        np = DimSize(w,0)
        
        for(ii=sp;ii<(np-sp/2);ii+=sp)
                for(jj=sp;jj<(np-sp/2);jj+=sp)
                        w[ii][jj] = 1
                endfor
        endfor
        
        
        return(0)
End

// sp is the 
Function FillPlaneHexagonal(w,sp)
        Wave w
        Variable sp             //spacing in pixels
        
        Variable np,ii,jj,xx,yy,fill
        
        np = DimSize(w,0)
        
        fill = 1
        // do even columns, then odds
        for(ii=sp;ii<(np-sp/2);ii+=(sp*sqrt(3)) )
                for(jj=sp;jj<(np-sp/2);jj+=sp)
                        w[ii][jj] = fill
                endfor
        endfor
        
        // now odds
        for(ii=(sp+sp*sqrt(3)/2);ii<(np-sp/2);ii+=(sp*sqrt(3)))
                for(jj=sp/2;jj<(np-sp/2);jj+=sp)
                        w[ii][jj] = fill
                endfor
        endfor
        
        return(0)
End

Proc HexagonalCylinders()

        Variable xc,yc,zc,len,rad,aa,cosa,sina,meanRad,pd,meanLen
        
        Variable grid=root:FFT_T
        Variable npt=root:FFT_N
        
        xc=npt/2
        yc=npt/2
        zc=npt/2
        
        meanRad=20
        pd=0.01
        meanLen=400
        len=meanLen
        
//      aa=meanRad/grid*2 +40   //grid units, +2 gives 2 units of separation
        aa=meanRad/grid*2 +2    //grid units, +2 gives 2 units of separation
        cosa = round(0.5*aa)
        sina = round(0.866*aa)
        

        Variable ii,nPass,x1,x2,x3
        ii=0
        nPass=1


        x1 = Gauss(meanRad-pd*meanRad,meanRad,pd*meanRad)
        x2 = Gauss(meanRad,meanRad,pd*meanRad)
        x3 = Gauss(meanRad+pd*meanRad,meanRad,pd*meanRad)
        
        Print x1,x2,x3
//      
/////////// average over 3 point distribution --
//      //calculate and sum the three contributions
//      FastOp mat=0
//      rad = meanRad-pd*meanRad
////    len = meanLen + gnoise(pd*meanLen)
//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
//      Execute "DoFFT()"
//      Duplicate/O iBin iBin_1g_128
//      iBin_1g_128 *= x1
//      
//      FastOp mat=0
//      rad = meanRad
////    len = meanLen + gnoise(pd*meanLen)
//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
//      Execute "DoFFT()"
//      iBin_1g_128 += x2*iBin
//      
//      FastOp mat=0
//      rad = meanRad+pd*meanRad
////    len = meanLen + gnoise(pd*meanLen)
//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
//      Execute "DoFFT()"
//      iBin_1g_128 += x3*iBin
//      
//      iBin_1g_128 /= (x1+x2+x3)               //renormalize
//      
////// end 3 pt average

//////////////////////
// paired cylinders
/////////////////////////
//      //calculate and sum the three contributions
//      FastOp mat=0
//      rad = meanRad-pd*meanRad
////    len = meanLen + gnoise(pd*meanLen)
//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
//      Execute "DoFFT()"
//      Duplicate/O iBin iBin_2g_256_sep40
//      iBin_2g_256_sep40 *= x1
//      
//      FastOp mat=0
//      rad = meanRad
////    len = meanLen + gnoise(pd*meanLen)
//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
//      Execute "DoFFT()"
//      iBin_2g_256_sep40 += x2*iBin
//      
//      FastOp mat=0
//      rad = meanRad+pd*meanRad
////    len = meanLen + gnoise(pd*meanLen)
//      FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
//      FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
//      Execute "DoFFT()"
//      iBin_2g_256_sep40 += x3*iBin
//      
//      iBin_2g_256_sep40 /= (x1+x2+x3)         //renormalize
        
/////////////////////////       
//
//
////
//      ii = 0  
//      do
//              FastOp mat=0
//              rad = meanRad + gnoise(pd*meanRad)
//              print ii,rad
//              FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)
//              Execute "DoFFT()"
//              if(ii==0)       
//                      Duplicate/O iBin iBin_1g_256
//              else
//                      iBin_1g_256 += iBin
//              endif
//              ii+=1
//      while(ii<nPass)
//      iBin_1g_256 /= nPass
//      
/////   
//      ii=0
//      nPass=1
//      do
//              FastOp mat=0
//              rad = meanRad + gnoise(pd*meanRad)
//              print ii,rad
//              FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)              //cylinder 1
//              rad = meanRad + gnoise(pd*meanRad)
//              print ii,rad
//              FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)   //cylinder 2
//              Execute "DoFFT()"
//              if(ii==0)       
//                      Duplicate/O iBin iBin_2g_256
//              else
//                      iBin_2g_256 += iBin
//              endif
//              ii+=1
//      while(ii<nPass)
//      iBin_2g_256 /= nPass
//      
        ii=0
        npass=1
        do
                FastOp mat=0
                rad = meanRad + gnoise(pd*meanRad)
                print ii,rad
                FillZCylinder(mat,grid,rad,xc,yc,zc,len,1)              //cylinder 1
                rad = meanRad + gnoise(pd*meanRad)
                print ii,rad
                FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)           //cylinder 2
                rad = meanRad + gnoise(pd*meanRad)
                print ii,rad
                FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,1)            //cylinder 3
                Execute "DoFFT()"
                if(ii==0)       
                        Duplicate/O iBin iBin_3g_256
                else
                        iBin_3g_256 += iBin
                endif
                ii+=1
        while(ii<nPass)
        iBin_3g_256 /= nPass
//      
        rad = meanRad + gnoise(pd*meanRad)
        print rad
        FillZCylinder(mat,grid,rad,xc,yc+aa,zc,len,1)
        Execute "DoFFT()"
        Duplicate/O iBin iBin_2
///     
        rad = meanRad + gnoise(pd*meanRad)
        print rad
        FillZCylinder(mat,grid,rad,xc+sina,yc+cosa,zc,len,1)
        Execute "DoFFT()"
        Duplicate/O iBin iBin_3
///     
        rad = meanRad + gnoise(pd*meanRad)
        print rad
        FillZCylinder(mat,grid,rad,xc+sina,yc-cosa,zc,len,1)
        Execute "DoFFT()"
        Duplicate/O iBin iBin_4
///     
        rad = meanRad + gnoise(pd*meanRad)
        print rad
        FillZCylinder(mat,grid,rad,xc,yc-aa,zc,len,1)
        Execute "DoFFT()"
        Duplicate/O iBin iBin_5
///     
        rad = meanRad + gnoise(pd*meanRad)
        print rad
        FillZCylinder(mat,grid,rad,xc-sina,yc-cosa,zc,len,1)
        Execute "DoFFT()"
        Duplicate/O iBin iBin_6
///     
        rad = meanRad + gnoise(pd*meanRad)
        print rad
        FillZCylinder(mat,grid,rad,xc-sina,yc+cosa,zc,len,1)
        Execute "DoFFT()"
        Duplicate/O iBin iBin_7

end


Proc Vary_N_In_Direction()

        FFTEraseMatrixButtonProc("")
        ConnectedRodFill(mat,2,30000*4,1)
        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
        Duplicate/O ival_XOP ival_XOP_2
        Duplicate/O qval_XOP qval_XOP_2
        
        FFTEraseMatrixButtonProc("")
        ConnectedRodFill(mat,3,30000*4,1)
        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
        Duplicate/O ival_XOP ival_XOP_3
        Duplicate/O qval_XOP qval_XOP_3
        
        FFTEraseMatrixButtonProc("")
        ConnectedRodFill(mat,4,30000*4,1)
        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
        Duplicate/O ival_XOP ival_XOP_4
        Duplicate/O qval_XOP qval_XOP_4
        
        FFTEraseMatrixButtonProc("")
        ConnectedRodFill(mat,5,30000*4,1)
        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
        Duplicate/O ival_XOP ival_XOP_5
        Duplicate/O qval_XOP qval_XOP_5
        
        FFTEraseMatrixButtonProc("")
        ConnectedRodFill(mat,6,30000*4,1)
        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
        Duplicate/O ival_XOP ival_XOP_6
        Duplicate/O qval_XOP qval_XOP_6
        
        FFTEraseMatrixButtonProc("")
        ConnectedRodFill(mat,7,30000*4,1)
        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
        Duplicate/O ival_XOP ival_XOP_7
        Duplicate/O qval_XOP qval_XOP_7

        FFTEraseMatrixButtonProc("")
        ConnectedRodFill(mat,8,30000*4,1)
        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
        Duplicate/O ival_XOP ival_XOP_8
        Duplicate/O qval_XOP qval_XOP_8
        
        FFTEraseMatrixButtonProc("")
        ConnectedRodFill(mat,9,30000*4,1)
        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
        Duplicate/O ival_XOP ival_XOP_9
        Duplicate/O qval_XOP qval_XOP_9
        
        FFTEraseMatrixButtonProc("")
        ConnectedRodFill(mat,10,30000*4,1)
        DoBinnedSpheresCalcFFTPanel(300,0.0001,0.5)
        Duplicate/O ival_XOP ival_XOP_10
        Duplicate/O qval_XOP qval_XOP_10
                
End