c Makes 2D projections for nonSH files, from the voxel representation
c does x and y and z projections
c sx is in the *-part-*.dat file along with all the voxels
c save(i,4) is just whether voxel is on surface or not, don't need
c that for this algorithm
c perimeter, FOurier coefficients, min/max Feret diameters computed

         real pi,sx,sy,lenx,leny,area
         integer save(10000000,4)

c needs to be changed with filename length changes
         character*33 filename
         character*39 filename2

         integer ierr,npes,Ntot,first,last,ios,first1
         integer nnn1
         integer n2,ind,ss,ss2
         integer rem
         integer*2 pix(500,500),burned

c USER
c set the number of nonSH files         
        Ntot=2390
c USER
c voxel length scale (micrometers) is sx

         pi=4.0*atan(1.0)
         nx=500
         ny=500

c USER
c replace characters before "-nonSH"
         open(unit=14,file='ASTM-AMPM-2-nonSH-proj.dat')
         open(unit=12,file='nonSH.lis')
         open(unit=7,file='ASTM-AMPM-2-nonSH-proj.out')

       do 3800 lmn=1,Ntot
        read(12,2999) filename
c USER
2999    format(a33)

         filename2(1:6)='nonSH/'
         filename2(7:39)=filename
         open(unit=19,file=filename2)
c  Read in voxel coordinates (middle of voxel)
c first read in number of voxels and scale (um)
         read(19,1003) isave,sx
1003     format(i10,f9.6)
         sy=sx
         vol=isave*sx*sx*sx
         do 390 n=1,isave
         read(19,*) save(n,1),save(n,2),save(n,3),save(n,4)
390      continue
         close(19)

c now loop over x, y, and z projections, write the results
c out separately so we have 3x the number of data points
c as there are particles

      do 8540 ijk=1,3
c area that proj computes has any holes still in it
      write(21,*) ' before call proj'
      call flush(21)
      call proj(pix,lmn,sx,isave,save,area,ijk,
     &          lenx,leny,xcm,ycm)
      write(21,*) ' after call proj'
      call flush(21)

c proj returns pix(i,j), which is the digital pixel projection
c of the particle. This is what anm needs.
c also have to first compute xcm,ycm,lenx,leny before 
c calling anm.
c anm uses small to fill in any holes
c in these multi-particle projections, then computes the 2D
c Fourier coefficients, then perimeter and min and max Feret diameters
c use area computed BEFORE holes are filled in
      burned=2
      write(21,*) ' before call anm'
      call anm(pix,filename,area,
     &         sx,sy,lenx,leny,nx,ny,xcm,ycm,burned,
     &          volpix,volume,numtot)
      write(21,*) ' after call anm'
      call flush(21)
c      write(14,399) filename,area
c399   format(a33,f18.10)
      call flush(12)
8540  continue

3800    continue
        
        end

      subroutine proj(pix,lmn,sx,isave,save,area,ijk,
     &      lenx,leny,xcm,ycm)

c Computes area of a projection in the x, y, or z direction
c controlled by value of ijk (1=x,2=y,3=z)

         real area,diam,sx,lenx,leny
         integer save(10000000,4)
         integer*2 pix(500,500)

         pi=4.0*atan(1.0)

c initialize projection image to matrix=1, projection will be 2
        do 100 j=1,500
        do 100 i=1,500
        pixii,j)=1
100     continue

c determine center of mass of particle, will move it to center of image
        icm=0
        jcm=0
        kcm=0
        do 1000 m=1,isave
        icm=icm+save(m,1)
        jcm=jcm+save(m,2)
        kcm=kcm+save(m,3)
1000    continue
        icm=icm/isave
        jcm=jcm/isave
        kcm=kcm/isave
        write(8,*) ' icm jcm kcm'
        write(8,*) icm,jcm,kcm
        ileft=10000
        jleft=10000
        iright=0
        jright=0
        do 1100 m=1,isave
c move image over to center of projected image (500x500)
c this will put center of mass at 250,250
        save(m,1)=save(m,1)-icm+250
        save(m,2)=save(m,2)-jcm+250
        save(m,3)=save(m,3)-kcm+250
        if(ileft.gt.save(m,1)) ileft=save(m,1)
        if(jleft.gt.save(m,2)) jleft=save(m,2)
        if(iright.le.save(m,1)) iright=save(m,1)
        if(jright.le.save(m,2)) jright=save(m,2)
1100    continue
c add lenx, leny computation, just like in part*label.f 
      xr=iright-250
      xl=250-ileft
      yr=jright-250
      yl=250-jleft
      lenx=2*max(xl,xr)
      leny=2*max(yl,yr)
      write(8,*) 'ij left right'
      write(8,*) ileft,iright,jleft,jright
      write(8,*) xl,xr,yl,yr
      write(8,*) 'lenx leny'
      write(8,*) lenx,leny
c check that projection will fit in digital image, centered at image
        if((iright-ileft).ge.500.or.(jright-jleft).ge.500.or.
     &    (kright-kleft).ge.500) then
        write(7,*) ileft,jleft,kleft,iright,jright,kright
        end if 

c now do projection
        area=0.0
        do 1300 m=1,isave
        if(ijk.eq.1) then
        pix(save(m,2),save(m,3))=2
        end if
        if(ijk.eq.2) then
        pix(save(m,1),save(m,3))=2
        end if
        if(ijk.eq.3) then
        pix(save(m,1),save(m,2))=2
        end if
1300    continue
c compute area
        xcm=0
        ycm=0
        do 1400 j=1,500
        do 1400 i=1,500
        if(pix(i,j).eq.2) then
        area=area+1.0
        xcm=xcm+i
        ycm=ycm+j
        end if
1400    continue
        xcm=xcm/area
        ycm=ycm/area
        write(8,*) ' proj cm = ',xcm,ycm
c convert area to real units, compute diam
        area=area*sx*sx

        return
        end
      subroutine anm(pix,filename,area,
     &         sx,sy,lenx,leny,nx,ny,xcm,ycm,burned,
     &          volpix,volume,numtot)

c  Interpolate surface function (rr(phi)) of shape 
c  (0 = matrix, ipart = particle). Take regular increments of phi
c  according to Gaussian quadrature points,
c  get rr(phi) corresponding to these points.
c  Go along (phi) vector (r unit vector), approximately 
c  find surface points by taking pairs inside and outside object.
c  Use this numerical characterization of r(phi) to compute
c  the Fourier coefficients
c  output into new file with labelling scheme

         real rf(1000),len
         real lenx,leny,area
         real x1,y1,x2,y2,z1,z2
c aa, bb are cos and sin Fourier coefficienti
         real xg(1000),wg(1000),pi,xx1,yy1,xx,yy,aa(0:1000),bb(0:1000)
         real hh,kk,rr1,rr2,meanrc,lineeq(1000000,6)
         real xtry,ytry,rtry,phi,phipt,xproj
         real endpts(10000,8),perim,kkurv
         real r,rphi,rphiphi,kurv(1000000),tan(1000000,3)
c lineeq: 1 = slope, 2 = intercept. 3,4 - x,y of one end of segment, 5,6 = other end
         real slope,bint,sx,sy
c added for xcmin computation
      real xc(1000,1000,4),xcmax(1000),yk(1000),yF(100000)
      real yB(100000),phiB(100000),phiF(100000),rBB(100000),rFF(100000)
      real rpp(1000)

         integer*2 pix(500,500)
         integer*2 burned
c USER
         character*34 filename
         character*31 filename2

c sx, sy are the scale factors

         pi=4.0*atan(1.0)
c USER
         ng=120
         open(unit=30,file='gauss120.dat')

         do 1100 i=1,ng
         read(30,*) xg(i),wg(i)
1100     continue
         close(30)

c  Interpolate r(phi) for a regular arrangement of (phi)
c  cm coords have been computed in subroutine fire
         
         xcm=xcm*sx
         ycm=ycm*sy
         write(8,*) ' xcm ycm scaled in anm =',xcm,ycm

         write(8,*) ' lenx leny in anm = ',lenx,leny
         len=min(lenx*sx,leny*sy)
         if(len.eq.0.0) then
         len=0.5*max(lenx*sx,leny*sy)
         write(8,*) ' len = 0 '
         end if
         write(8,*) ' len = ',len
         cutoff=len/100000.0
         write(8,*) ' cutoff in anm = ',cutoff
         write(8,*) ' start interpolation numtot = ',numtot
         call flush(8)

c need to call small before interpolation and computation of 
c Fourier coefficients and r(phi)
         write(21,*) ' before call small'
         call flush(21)
         call small(pix,nx,ny,burned,
     &            bcount,ipore,sx,sy)
         write(21,*) ' after call small'
         call flush(21)

         write(21,*) ' before interpolation'
         call flush(21)
         do 200 j=1,ng
         phi=pi*xg(j)+pi
         phi=pi-phi
         dd=len/10.
         rstart=0.0

45       continue
         do 40 ii=1,10
         rr1=rstart+(ii-1)*dd
         rr2=rr1+dd
c  position with respect to center of mass
         x1=rr1*cos(phi)
         y1=rr1*sin(phi)
         x2=rr2*cos(phi)
         y2=rr2*sin(phi)

c  actual coordinate
         x1=x1+xcm
         x2=x2+xcm
         y1=y1+ycm
         y2=y2+ycm

c  nearest pixel center
         i1=int(x1/sx+0.5)
         i2=int(x2/sx+0.5)
         j1=int(y1/sy+0.5)
         j2=int(y2/sy+0.5)

c  check if 2nd point goes beyond image box, on the first try
           if(i1.gt.nx.or.i1.lt.1.or.j1.gt.ny.or.j1.lt.1) then
             write(8,*) ' i,j beyond image bd for 1',i1,j1
c             dd=dd/2.0
c             goto 45
           end if
           if(i2.gt.nx.or.i2.lt.1.or.j2.gt.ny.or.j2.lt.1) then
             dd=dd/2.0
             write(8,*) ' i,j beyond image bd for 2 ',i2,j2
             write(8,*) ' rr1 rr2 ',rr1,rr2
             goto 45
           end if

c  if pair of points are both within particle, go to next pair
      if(pix(i1,j1).eq.burned.and.pix(i2,j2).eq.burned) then
           if(ii.lt.10) goto 40
             if(ii.eq.10) then
             rstart=rr2
             goto 45
             end if
           end if

c  if pair crossed boundary, reduce step size, rename rstart at first point
c  of pair, and do over.  If step size is small enough, then end.
         if(pix(i1,j1).eq.burned.and.pix(i2,j2).ne.burned) then
          if(abs(dd).lt.cutoff) then
          goto 49
          end if
          rstart=rr1
          dd=dd/10.0
          goto 45
         end if

40       continue

49       rf(j)=0.5*(rr1+rr2)

200      continue
         write(21,*) ' after interpolation'
         call flush(21)

c  compute Fourier expansion of surface

c rf(j) is the 2D set of points on the projection surface
c now fit Fourier series to these values
         
c  nnnF is how many sin and cos terms to use in series (n limit)
c bb(0)=0.0
         nnnF=26
         do 391 n=0,nnnF
         aa(n)=0.0
         bb(n)=0.0
391      continue

         do 831 j=1,ng
         phi=pi*xg(j)+pi
         aa(0)=aa(0)+rf(j)*wg(j)
         do 400 n=1,nnnF
         aa(n)=aa(n)+rf(j)*wg(j)*cos(n*phi)
         bb(n)=bb(n)+rf(j)*wg(j)*sin(n*phi)
400      continue

831      continue

         do 832 j=1,ng
         phi=pi*xg(j)+pi
         r=aa(0)*0.5
         do 402 n=1,nnnF
         r=r+aa(n)*cos(n*phi)+bb(n)*sin(n*phi)
402      continue
         xx1=r*cos(phi)
         yy1=r*sin(phi)
         write(8,973) j,rf(j),r,xx1,yy1
973      format(i4,4f14.8)
         
832      continue

         perim=0.0

         factor=pi
c redo x and y extents, since projection is different from a slice
         x1=10000.
         y1=10000.
         x2=0.0
         y2=0.0
         do 430 j=1,ng
         phi=pi*xg(j)+pi

         r=aa(0)*0.5
         rphi=0.0
         rphiphi=0.0
c  derivatives of(0,0) term is zero, since it is a constant
         do 62 n=1,nnnF
         r=r+aa(n)*cos(n*phi)+bb(n)*sin(n*phi)
         rphi=rphi+n*(-aa(n)*sin(n*phi)+bb(n)*cos(n*phi))
         rphiphi=rphiphi-n*n*(aa(n)*cos(n*phi)+bb(n)*sin(n*phi))
62       continue
         xx1=r*cos(phi)
         yy1=r*sin(phi)
         if(x1.gt.xx1) x1=xx1*1.00001
         if(y1.gt.yy1) y1=yy1*1.00001
         if(x2.lt.xx1) x2=xx1*1.00001
         if(y2.lt.yy1) y2=yy1*1.00001
         
         perim=perim+wg(j)*sqrt(rphi**2+r**2)
430      continue
         dmin=0.25*min(x2-x1,y2-y1)
         perim=perim*factor
         write(7,*) 'perim ',filename2,perim
         write(7,*) 'perim ',perim/pi


c go around object, use Nj angles, remember those with positive curvature
         Nj=5000
         jcurv=0
c initialize lineeq line segment end points so I know which have been used
         do 489 i=1,1000000
         lineeq(i,3)=-10000.
         lineeq(i,4)=-10000.
         lineeq(i,5)=-10000.
         lineeq(i,6)=-10000.
489      continue

         do 530 j=1,Nj
         phi=(j-1)*2.*pi/float(Nj)

         r=aa(0)*0.5
         rphi=0.0
         rphiphi=0.
c  derivatives of (0,0) term is zero, since it is a constant
         do 32 n=1,nnnF
         r=r+aa(n)*cos(n*phi)+bb(n)*sin(n*phi)
         rphi=rphi+n*(-aa(n)*sin(n*phi)+bb(n)*cos(n*phi))
         rphiphi=rphiphi-n*n*(aa(n)*cos(n*phi)+bb(n)*sin(n*phi))
32       continue
         xx1=r*cos(phi)
         yy1=r*sin(phi)
         if(xx1.lt.x1) x1=1.00001*xx1
         if(xx1.gt.x2) x2=1.00001*xx1
         if(yy1.lt.y1) y1=1.00001*yy1
         if(yy1.gt.y2) y2=1.00001*yy1
         
         kkurv=(r*r+2.*(rphi)**2-r*rphiphi)
         kkurv=kkurv/(r*r+rphi**2)**(1.5)

         if(kkurv.ge.0.0) then
          jcurv=jcurv+1
          kurv(jcurv)=kkurv
c work out tangent unit vector
          tan(jcurv,1)=(rphi*cos(phi)-r*sin(phi))/sqrt(rphi**2+r**2)
          tan(jcurv,2)=(rphi*sin(phi)+r*cos(phi))/sqrt(rphi**2+r**2)
          tan(jcurv,3)=tan(jcurv,2)/tan(jcurv,1)
c work out slope and y-intercept of tangent line
          slope=tan(jcurv,3)
          bint=r*sin(phi)-tan(jcurv,3)*r*cos(phi)
          lineeq(jcurv,1)=slope
          lineeq(jcurv,2)=bint

c check intersection at top or bottom, then left or right
c this assumes that slope is never vertical or zero

c this is the first if statement, no end points have yet been found, so if one 
c is found here, it goes in (3,4)
          xx=(y1-bint)/slope
          if(xx.ge.x1.and.xx.le.x2) then
           lineeq(jcurv,3)=xx
           lineeq(jcurv,4)=y1
          end if

          xx=(y2-bint)/slope
          if(xx.ge.x1.and.xx.le.x2) then
           if(lineeq(jcurv,3).eq.-10000.) then
            lineeq(jcurv,3)=xx
            lineeq(jcurv,4)=y2
            goto 1200
           end if
            lineeq(jcurv,5)=xx
            lineeq(jcurv,6)=y2
          end if

1200      yy=slope*x1+bint
          if(yy.ge.y1.and.yy.le.y2) then
           if(lineeq(jcurv,3).eq.-10000.) then
            lineeq(jcurv,3)=x1
            lineeq(jcurv,4)=yy
            goto 1300
           end if
           if(lineeq(jcurv,5).eq.-10000.) then
            lineeq(jcurv,5)=x1
            lineeq(jcurv,6)=yy
           end if
          end if

1300      yy=slope*x2+bint
          if(yy.ge.y1.and.yy.le.y2) then
           if(lineeq(jcurv,3).eq.-10000.) then
            lineeq(jcurv,3)=x2
            lineeq(jcurv,4)=yy
           end if
           if(lineeq(jcurv,5).eq.-10000.) then
            lineeq(jcurv,5)=x2
            lineeq(jcurv,6)=yy
           end if
          end if

c now endpoints of line segment are established, determine if points on line
c segment are inside particle or not. If not, then keep this point. If yes,
c then throw away this point by letting jcurv=jcurv-1
c does not check actual endpoints, only points in-between
c assumes origin is at center of particle (which is how the extents, anm, etc. are defined

         nover=100
         do 475 it=1,nover-1
         t=float(it)/float(nover)
         xtry=lineeq(jcurv,3)*(1.-t)+t*lineeq(jcurv,5)
         ytry=lineeq(jcurv,4)*(1.-t)+t*lineeq(jcurv,6)
         rtry=sqrt(xtry**2+ytry**2)
       phipt=atan( ytry/xtry )
       if(ytry.lt.0.0.and.xtry.lt.0.0) phipt=phipt+pi
       if(ytry.gt.0.0.and.xtry.lt.0.0) phipt=phipt+pi
       if(ytry.lt.0.0.and.xtry.gt.0.0) phipt=phipt+2.*pi
       r=0.5*aa(0)
       do 20 n=1,nnnF
       r=r+aa(n)*cos(n*phipt)+bb(n)*sin(n*phipt)
20     continue
        if(rtry.le.r) then
         lineeq(jcurv,3)=-10000.
         lineeq(jcurv,4)=-10000.
         lineeq(jcurv,5)=-10000.
         lineeq(jcurv,6)=-10000.
         jcurv=jcurv-1
        goto 530
        end if
475      continue

c this is the end of curvature > 0 if then statement
         end if

530      continue

c if program has made it to here, then jcurv is the total number of slopes at positive
c curvature points, and line segment to bounding box does not overlap the particle
c now compare pairs of these to find calipers, take maximum of pairs for each slope considered
c to get the caliper distance for that slope of tangent lines
         write(7,*) ' number of good points = ',jcurv
         npp=0
         distmax=-10.0
         distmin=10000.0
         do 610 j=1,jcurv-1
         np=0
c comparing slopes, if difference is more than 1 %, go to next pair
         do 600 jj=j+1,jcurv
         diff=abs( (lineeq(j,1)-lineeq(jj,1))/lineeq(j,1))
         if(diff.gt.0.01) goto 600
         dist=abs(lineeq(j,2)-lineeq(jj,2))/sqrt(1.+lineeq(j,1)**2)
         if(dist.lt.dmin) goto 600
         if(distmax.lt.dist) then
         np=np+1
         distmax=dist
         j1=j
         j2=jj
         end if
         if(distmin.gt.dist) then
         np=np+1
         distmin=dist
         jj1=j
         jj2=jj
         end if
600      continue

610      continue

         Fermax=distmax
         Fermin=distmin

c now compute xcmin
c xc(m,k,4) m=phi0s, k=y's, 1=phiF, 2=phiB (x - coordinates)
c xcmax(m) = max over k of abs(xc(m,k,1)-xc(m,k,2))
c xcmin = min over m of xcmax(m)

c loop in rotation angle phi0 (index m)
c only need to do 0-pi, since particle up is the same as particle down
c this is different orientations of the projection
c Note - think of this as: rprime(phi)=r(phi-phi0)
c so when I compute the corresponding y value, I use phi, not phi-phi0
       Mphi0=200
       do 4180 m=1,Mphi0
c rotation angle
       phi0=(m-1)*pi/float(Mphi0)

c compute ymin and ymax and angles at these values for given phi0
       ymin=10000.
       ymax=0.0
c Nex = number of angles used to calculate ymin and ymax
c phi is the true angle, phi-phi0 gives the rotated version
       Nex=1000
       do 4130 j=1,Nex
       phi=(j-1)*2.*pi/float(Nex)
       r=aa(0)*0.5
       ppp=phi-phi0
       do 612 n=1,nnnF
       r=r+aa(n)*cos(n*ppp)+bb(n)*sin(n*ppp)
612    continue
       yy1=r*sin(phi)
       if(ymin.gt.yy1) then
       ymin=yy1
       phimin=phi
       end if
       if(ymax.lt.yy1) then
       ymax=yy1
       phimax=phi
       end if
4130   continue

c  ymin,phimin and ymax,phimax are now computed for m (phi0) orientation
c lay out yk(k) array for KK-1 points
       KK=200
       do 4131 k=1,KK-1
       yk(k)=ymin+k*(ymax-ymin)/float(KK)
       if(yk(k).eq.0.0) yk(k)=yk(k)+0.1*(ymax-ymin)/float(KK)
4131   continue

c now determine angles needed to match set of yk(k)'s
c start with angle phimin, and go both ways until reaching phimax
      NNphi=10000
c F = forward (counter clockwise) B = backwards (clockwise)
c increment in angle for backward direction
      phiincB=abs(phimax-phimin)/float(NNphi)
c increment in angle for forward direction
      phiincF=(2.*pi-abs(phimax-phimin))/float(NNphi)
      kF=0
      kB=0
      do 4150 nphi=1,NNphi-1
      phiF(nphi)=phimin+nphi*phiincF
      phiB(nphi)=phimin-nphi*phiincB
      rFF(nphi)=aa(0)*0.5
      rBB(nphi)=aa(0)*0.5
c need to rotate angles when calculating rFF and rBB
      aaa=phiF(nphi)-phi0
      bbb=phiB(nphi)-phi0
      do 613 n=1,nnnF
      rFF(nphi)=rFF(nphi)+aa(n)*cos(n*aaa)+bb(n)*sin(n*aaa)
      rBB(nphi)=rBB(nphi)+aa(n)*cos(n*bbb)+bb(n)*sin(n*bbb)
613   continue
c use unrotated angle when computing y coordinate
      yF(nphi)=rFF(nphi)*sin(phiF(nphi))
      yB(nphi)=rBB(nphi)*sin(phiB(nphi))

      if(yF(nphi).lt.yk(kF+1)) then
       if(yB(nphi).lt.yk(kB+1)) then
       goto 4150
       end if
      end if

      if(kF.eq.(KK-1)) goto 4114
      if(yF(nphi).gt.yk(kF+1)) then
      xm=(yF(nphi)-yF(nphi-1))/(phiF(nphi)-phiF(nphi-1))
      phiFsave=phiF(nphi-1)+(yk(kF+1)-yF(nphi-1))/xm
c compute r(phiFsave)
      aaa=phiFsave-phi0
      rrF=aa(0)*0.5
      do 614 n=1,nnnF
      rrF=rrF+aa(n)*cos(n*aaa)+bb(n)*sin(n*aaa)
614   continue
      kF=kF+1
      xc(m,kF,1)=rrF*cos(phiFsave)
      xc(m,kF,3)=rrF*sin(phiFsave)
      end if

4114  if(kB.eq.(KK-1)) goto 4115
      if(yB(nphi).gt.yk(kB+1)) then
      xm=(yB(nphi)-yB(nphi-1))/(phiB(nphi)-phiB(nphi-1))
      phiBsave=phiB(nphi-1)+(yk(kB+1)-yB(nphi-1))/xm
c compute r(phiBsave)
      bbb=phiBsave-phi0
      rrB=aa(0)*0.5
      do 617 n=1,nnnF
      rrB=rrB+aa(n)*cos(n*bbb)+bb(n)*sin(n*bbb)
617   continue
      kB=kB+1
      xc(m,kB,2)=rrB*cos(phiBsave)
      xc(m,kB,4)=rrB*sin(phiBsave)
      end if
4115  if(kB.eq.(KK-1).and.kF.eq.(KK-1)) goto 5555

4150  continue
c5555  if(me.eq.0) write(33,*) 'kB kF after 4150 loop ',kB,kF
5555  continue

c now find maximum over kB and kF to determine xcmax(m)
      xcmax(m)=0.0
      do 4155 k=1,KK-1
      dd=abs(xc(m,k,1)-xc(m,k,2))
      if(xcmax(m).lt.dd) xcmax(m)=dd
c      write(32,*) filename,area,xc(m,k,1),xc(m,k,3),xc(m,k,2),
c     &      xc(m,k,4),xcmax(m) 
4155  continue
c go to another rotation angle
4180  continue

      xcmin=10000.0
c xcmm = max of xcmax
      xcmm=0.0
      do 4190 m=1,Mphi0
      if(xcmin.gt.xcmax(m)) xcmin=xcmax(m)
      if(xcmm.lt.xcmax(m)) xcmm=xcmax(m)
4190  continue

c now xcmin has been calculated, write out all results

      write(14,122) filename,area,perim,Fermin,Fermax,xcmin,xcmm,nnnF
122   format(a33,6f18.6,i4)
      call flush(14)
c end of x,y,z projection loopfor a given particle
3888     continue

3800     continue
c end of particle loop

2201     continue
         return
         end


      subroutine small(pix,nx,ny,burned,
     &            bcount,ipore,sx,sy)

      integer*2 pix(500,500)
      integer*2 old(1000000,2),new(1000000,2)
      integer*2 phase,burned,irestore,tempburn
      integer*2 in(4),jn(4)
      integer*4 inew,iold
c Start with a box containing the 2D projection. 
c Turn all outer lines into tempburn, then 
c build a "particle" or fire inward from each side.
c When "particle" or fire is stopped,
c then just need to turn inner real pores into burned particle phase

      pi=4.0*atan(1.0)
c  Direction labels to check for burning path (nearest neighbor information
c  in 2-D digital system).
      in(1)=-1
      in(2)=1
      in(3)=0
      in(4)=0
      jn(1)=0
      jn(2)=0
      jn(3)=-1
      jn(4)=1

      tempburn=-1500
c      write(33,*) ' in small ',nx,ny,burned,tempburn
c  find first matrix pixels around faces of particle system
      write(21,*) ' before setting up burned frame'
      do 990 j=1,ny
      i=1
      pix(i,j)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      i=nx
      pix(i,j)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
990   continue
      do 991 i=1,nx
      j=1
      pix(i,j)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
      j=ny
      pix(i,j)=tempburn
      iold=iold+1
      old(iold,1)=i
      old(iold,2)=j
991   continue

c  Now start building up new pixels from old set of pixels,
c  thus propagating the fire inward and identifying new particle 
c  pixels connnected to the ones already burned.

      write(21,*) ' before fire'
      icount=1
60    inew=0
      do 100 ijk=1,iold
      ii=old(ijk,1)
      jj=old(ijk,2)
c  check all six nearest neighbors of previously burned pixel
      do 90 n=1,4
      i1=ii+in(n)
      j1=jj+jn(n)

c  Hard boundary conditions (can't pass beyond boundary)

      if(j1.lt.1.or.j1.gt.ny) goto 90
      if(i1.lt.1.or.i1.gt.nx) goto 90

c  Store (i,j) labels of newly burned pixels in array new().
      if(pix(i1,j1).ne.burned.and.
     &      pix(i1,j1).ne.tempburn) then
      pix(i1,j1)=tempburn
      inew=inew+1
      new(inew,1)=i1
      new(inew,2)=j1
      end if
90    continue
100   continue

c  If new pixels were burned, then transfer labels to old() array, start
c  burning process over again.
      if(inew.gt.0) then
      iold=inew
      do 150 ijk=1,inew
      old(ijk,1)=new(ijk,1)
      old(ijk,2)=new(ijk,2)
150   continue
      icount=icount+1
      goto 60
      end if

6060    continue
      write(21,*) ' before set internal pores to solid'

c set internal pores to solid
c these were set to 1 in projection subroutine
       do 225 j=1,ny
       do 225 i=1,nx
       if(pix(i,j).eq.1) then
       pix(i,j)=burned
       end if
225    continue

      return
      end