c analyzes particle porosities stored in *porosity.dat files

      real porosity(100000,2)
      character*10 ffile(100000)
      character*10 fSH,fSHmax
      character*10 fnonSH,fnonSHmax
      open(unit=8,file='ASTM-AMPM-2-SH-internal-porosity.dat')
      open(unit=9,file='ASTM-AMPM-2-nonSH-internal-porosity.dat')
      open(unit=7,file='ASTM-AMPM-2-intern-poros-analysis.txt')
      open(unit=10,file='ASTM-AMPM-2-intern-poros-list.txt')

      Nsh=14680
      Nnonsh=2390

c loop in both nonSH and SH particles, do separate statistics
c and together

c number of porous particles, average porosity,
c and total pore volume
c divided by total volume of porous particles

c total count for porous particles
      ncount=0
c porecount(n,3)=1 for SH, 0 for nonSH
c do for SH particles
      SH=0.0
      voltotSH=0.0
      avgSH=0.0
      avgtotSH=0.0
      xmaxporeSH=0.0
      xminporeSH=10000.0
      do 100 i=1,Nsh
      read(8,1) fSH,ivol,pore
1     format(13x,a10,7x,i10,23x,f12.8)
      if(pore.gt.0.000000) then
      ncount=ncount+1
      porosity(ncount,1)=ivol
      porosity(ncount,2)=pore
      ffile(ncount)=fSH
c      print *,ffile(ncount)
       SH=SH+1.0
       if(xmaxporeSH.lt.pore) then
        xmaxporeSH=pore
        fSHmax=fSH
       end if
       if(xminporeSH.gt.pore) xminporeSH=pore
       voltotSH=voltotSH+ivol
       avgtotSH=avgtotSH+pore*ivol
       avgSH=avgSH+pore
      end if
100   continue
      avgtotSH=avgtotSH/voltotSH
      avgSH=avgSH/SH
      write(7,*) ' SH '
      write(7,*) ' num particles with porosity = ',SH
      write(7,*) ' total SH particles  = ',Nsh
      write(7,*) ' avg porosity = ',avgSH,avgtotSH
      write(7,*) ' max porosity ',fSHmax,xmaxporeSH
      write(7,*) 

c now do for nonSH particles
      SHnon=0.0
      voltotnonSH=0.0
      avgnonSH=0.0
      avgtotnonSH=0.0
      xmaxporenonSH=0.0
      xminporenonSH=10000.0
      do 200 i=1,Nnonsh
      read(9,2) fnonSH,ivol,pore
2     format(15x,a10,11x,i10,23x,f12.8)
      if(pore.gt.0.000000) then
      ncount=ncount+1
      porosity(ncount,1)=ivol
      porosity(ncount,2)=pore
      ffile(ncount)=fnonSH
c      print *,ffile(ncount)
       SHnon=SHnon+1.0
       if(xmaxporenonSH.lt.pore) then
        xmaxporenonSH=pore
        fnonSHmax=fnonSH
       end if
       if(xminporenonSH.gt.pore) xminporenonSH=pore
       voltotnonSH=voltotnonSH+ivol
       avgtotnonSH=avgtotnonSH+pore*ivol
       avgnonSH=avgnonSH+pore
      end if
200   continue
      avgtotnonSH=avgtotnonSH/voltotnonSH
      avgnonSH=avgnonSH/SHnon
      write(7,*) ' nonSH '
      write(7,*) ' num particles with porosity = ',SHnon
      write(7,*) ' total SH particles  = ',Nnonsh
      write(7,*) ' avg porosity = ',avgnonSH,avgtotnonSH
      write(7,*) '  max porosity ',fnonSHmax,xmaxporenonSH
      write(7,*) 

c combine statistics
      totpore=SH+SHnon
      voltot=voltotSH+voltotnonSH
      avgtot=(avgtotSH*voltotSH+avgtotnonSH*voltotnonSH)/voltot
      avg=(avgSH*SH+avgnonSH*SHnon)/totpore
      total=Nsh+Nnonsh

      write(7,*) ' total nonSH+SH '
      write(7,*) ' num particles with porosity = ',totpore
      write(7,*) ' total num particles  = ',total
      write(7,*) ' avg porosity = ',avg,avgtot
      write(7,*) 

c list of all porous particles
      do 300 i=1,ncount
      write(10,5) ffile(i),porosity(i,1),porosity(i,2)
5     format(a11,f12.0,f12.8)
300   continue

      end