c generate volume-based distribution based on area in 2D projections
c do all particles together
c various 2D parameters

      real psd(1000,2),volume(1000000)
      real sa2d(1000000)
      real sieve(1000),sa3d(1000000),vd2d(1000000)
      real vd3d(1000000),xll(26),xww(26),xtt(26)
      real diam(1000000),xn(1000000),vesdavg(26)
      real diamsin(1000000),diammul(1000000)
      real size(30),volsize(30),ratio(1000000)
      real nSHvxlavg,xlavg

      integer nn(26),gaussI
      character*50 filename
      character*9 letter

      do 6600 mmm=1,6
      if(mmm.eq.1) then
      open(unit=7,file='ASTM-AMPM-2-AECD-sin-mult-all.txt')
      open(unit=12,file='ASTM-AMPM-2-AECD-hist-all.txt')
      write(7,*) ' AECD '
      letter='   AECD  '
       end if

      if(mmm.eq.2) then
      open(unit=7,file='ASTM-AMPM-2-Xcmin-sin-mult-all.txt')
      open(unit=12,file='ASTM-AMPM-2-Xcmin-hist-all.txt')
      write(7,*) ' Xcmin '
      letter='  Xcmin  '
      end if

      if(mmm.eq.3) then
      open(unit=7,file='ASTM-AMPM-2-Xcmax-sin-mult-all.txt')
      open(unit=12,file='ASTM-AMPM-2-Xcmax-hist-all.txt')
      write(7,*) ' Xcmax '
      letter='  Xcmax  '
      end if

      if(mmm.eq.4) then
      open(unit=7,file='ASTM-AMPM-2-AR-sin-mult-all.txt')
      open(unit=12,file='ASTM-AMPM-2-AR-hist-all.txt')
      write(7,*) ' AR=Fermax/Fermin '
      letter='    AR   '
      end if

      if(mmm.eq.5) then
      open(unit=7,file='ASTM-AMPM-2-C-sin-mult-all.txt')
      open(unit=12,file='ASTM-AMPM-2-C-hist-all.txt')
      write(7,*) ' C=4piA/P/P '
      letter='    C    '
      end if

      pi=4.0*atan(1.0)
      numtot=0
      voltot=0.0
      vtsin=0.0
      vtmul=0.0
      voltotSH=0.0

      open(unit=19,file='ASTM-AMPM-2-nonSH-proj.dat')
      open(unit=9,file='ASTM-AMPM-2-proj.dat')

c SH
      N=43743
c non-SH
      N1=7170

c number-based volume average for only the SH particles
c      volnumavg=0.0

      avg=0.0
      voltot=0.0
      avgl=0.0
      avgl2=0.0

c find min and max diameter for all particles
      diammax=0.0
      diammin=10000.

c all particles count
      ic=0

c do non-SH particles first
      do 109 i=1,N1
      read(19,*) filename,areai,perim,
     &              fermin,fermax,xcmin,xcmax,nnnF
      vd3di=(4.*areai/pi)**(1./2.)

       ic=ic+1
      volume(ic)=areai
      aecd=sqrt(4.*areai/pi)
      vd3d(ic)=aecd
      voltot=voltot+volume(ic)
      xlavg=xlavg+areai*vd3di

      if(mmm.eq.1)  diam(ic)=vd3d(ic)
      if(mmm.eq.2)  diam(ic)=xcmin
      if(mmm.eq.3)  diam(ic)=xcmax
      if(mmm.eq.4)  diam(ic)=fermax/fermin
      if(mmm.eq.5)  diam(ic)=4.*pi*areai/perim/perim

      avg=avg+areai*diam(ic)
      if(diammax.lt.diam(ic)) diammax=diam(ic)
      if(diammin.gt.diam(ic)) diammin=diam(ic)

109   continue
      close(19)
      voltotr=voltot
      d1=diammin
      d2=diammax

      do 100 i=1,N
      read(9,*) filename,areai,perim,
     &              fermin,fermax,xcmin,xcmax,nnnF

      ic=ic+1
      volume(ic)=areai
      vd3di=(4.*areai/pi)**(1./2.)
      vd3d(ic)=vd3di
      xlavg=xlavg+areai*vd3di

      if(mmm.eq.1)  diam(ic)=vd3di
      if(mmm.eq.2)  diam(ic)=xcmin
      if(mmm.eq.3)  diam(ic)=xcmax
      if(mmm.eq.4)  diam(ic)=fermax/fermin
      if(mmm.eq.5)  diam(ic)=4.*pi*areai/perim/perim
      avg=avg+areai*diam(ic)

      if(diammax.lt.diam(ic)) diammax=diam(ic)
      if(diammin.gt.diam(ic)) diammin=diam(ic)

      voltot=voltot+volume(ic)
100   continue

      write(7,*) ' number of SH particles = ',N
      write(7,*) ' number of nonSH particles = ',N1
      write(7,*) ' total # of particles ',ic
      write(7,*)
      close(9)
      N=ic
      xlavg=xlavg/voltot
      avg=avg/voltot

      write(7,*) ' min max ',letter,' for all particles'
      write(7,*) diammin,diammax
      write(7,*) ' average area-wtd AECD = ',xlavg
      write(7,*) ' average area-wtd ',letter,avg
      write(7,*) 

c generate diff. prob. based on volume fraction
c      nbins=100
      nbins=40

c do all particles together
      write(7,*) ic,' all particles'
c set diammin and diammax to be the same for all 5 particle types
c for L and L/W
c      if(mmm.eq.1) then
c      diammax=155.
c      diammin=9.0
c      end if
      dx=(diammax-diammin)/float(nbins)
      write(7,*) nbins,dx
      do 356 i=1,1000
      psd(i,1)=0.0
      psd(i,2)=0.0
356   continue

c do differential probability distribution for diam choice

      do 519 i=1,ic
        m=1+(diam(i)-diammin)/dx
        if(m.gt.nbins) m=nbins
        psd(m,1)=psd(m,1)+volume(i)/voltot
        psd(m,2)=psd(m,2)+1.0
519   continue

      sum=0.0
      do 366 i=1,nbins
      sum=sum+psd(i,1)
      x=diammin+(i-1)*dx
      write(7,113) x,x+dx,psd(i,1),psd(i,2)
113   format(3f11.6,f8.0)
366   continue
      write(7,*) sum

      write(12,129) diammin,0.0
129   format(2f11.6)
      do 466 i=1,nbins
      x=diammin+(i-1)*dx
      write(12,129) x,psd(i,1)
      write(12,129) x+dx,psd(i,1)
      write(12,129) x+dx,0.0
466   continue
      close(7)
      close(12)
6600  continue
      end