c generate LWT distributions
c divide particles by single (SH, L/T < SHcut, nonSH L/T < nSHcut)
c  and non-spherical  (SH L/T > SHcut, nonSH L/T > nSHcut)
c NOTE: non-spherical consists of mostly double and multi-particles,
c but also misshapen and broken single particles
c later maybe will try dividing the non-spherical particle into
c double, triple, and so-on, using the knowledge I am getting from
c the non-spherical sphere simulations

      real psd(1000,2),trI(1000000),volume(1000000)
      real volsin(1000000),volmul(1000000)
      real xw(1000000),xl(1000000),xt(1000000),sa2d(1000000)
      real sieve(1000),sa3d(1000000),vd2d(1000000),xcmin(1000000)
      real vd3d(1000000),xll(26),xww(26),xtt(26)
      real diam(1000000),xn(1000000),vesdavg(26),curv(1000000)
      real diamsin(1000000),diammul(1000000)
      real size(30),volsize(30),ratio(1000000)
      real nonSHavgl,nonSHavgw,nonSHavgt,nonSHavgV
      real nonSHavgl2,nonSHavgw2,nonSHavgt2,nonSHavgV2
      real nSHvxlwavg,numLT
      real nSHvxlavg,nSHcut
      real nSHvxwavg
      real nSHvxlavg2
      real nSHvxwavg2
      real nSHvxlwavg2

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

c cutoff values for SH and nonSH for single vs. non-spherical
      SHcut=1.17
      nSHcut=1.10
c # SH
      NS=14581
c # non-SH
      N1=2390

      do 5544 mmm=1,7

      if(mmm.eq.1) then
      open(unit=7,file='ASTM-AMPM2-VESD-sin-mult-all.txt')
      open(unit=10,file='ASTM-AMPM2-VESD-sin-hist-all.txt')
      open(unit=11,file='ASTM-AMPM2-VESD-mult-hist-all.txt')
      open(unit=12,file='ASTM-AMPM2-VESD-hist-all.txt')
      write(7,*) ' VESD '
      letter=' VESD'
       end if

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

      if(mmm.eq.3) then
      open(unit=7,file='ASTM-AMPM2-W-sin-mult-all.txt')
      open(unit=10,file='ASTM-AMPM2-W-sin-hist-all.txt')
      open(unit=11,file='ASTM-AMPM2-W-mult-hist-all.txt')
      open(unit=12,file='ASTM-AMPM2-W-hist-all.txt')
      write(7,*) ' W '
      letter='  W  '
      end if

      if(mmm.eq.4) then
      open(unit=7,file='ASTM-AMPM2-T-sin-mult-all.txt')
      open(unit=10,file='ASTM-AMPM2-T-sin-hist-all.txt')
      open(unit=11,file='ASTM-AMPM2-T-mult-hist-all.txt')
      open(unit=12,file='ASTM-AMPM2-T-hist-all.txt')
      write(7,*) ' T '
      letter='  T  '
      end if

      if(mmm.eq.5) then
      open(unit=7,file='ASTM-AMPM2-LW-sin-mult-all.txt')
      open(unit=10,file='ASTM-AMPM2-LW-sin-hist-all.txt')
      open(unit=11,file='ASTM-AMPM2-LW-mult-hist-all.txt')
      open(unit=12,file='ASTM-AMPM2-LW-hist-all.txt')
      write(7,*) ' L/W '
      letter=' L/W '
      end if

      if(mmm.eq.6) then
      open(unit=7,file='ASTM-AMPM2-WT-sin-mult-all.txt')
      open(unit=10,file='ASTM-AMPM2-WT-sin-hist-all.txt')
      open(unit=11,file='ASTM-AMPM2-WT-mult-hist-all.txt')
      open(unit=12,file='ASTM-AMPM2-WT-hist-all.txt')
      write(7,*) ' W/T '
      letter=' W/T '
      end if

      if(mmm.eq.7) then
      open(unit=7,file='ASTM-AMPM2-LT-sin-mult-all.txt')
      open(unit=10,file='ASTM-AMPM2-LT-sin-hist-all.txt')
      open(unit=11,file='ASTM-AMPM2-LT-mult-hist-all.txt')
      open(unit=12,file='ASTM-AMPM2-LT-hist-all.txt')
      write(7,*) ' L/T '
      letter=' L/T '
      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-AMPM2-nonSH-len.dat')
      open(unit=9,file='ASTM-AMPM2-geom-len.dat')

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

c volume-based ratio averages for all particles
c L/T, W/T, L/W
      vxlwavg=0.0
      vxlavg=0.0
      vxwavg=0.0
      vxlavg2=0.0
      vxwavg2=0.0
      vxlwavg2=0.0
c volume-based ratio averages for single SH+nonSH particles
c L/T, W/T, L/W
      SHvxlwavg=0.0
      SHvxlavg=0.0
      SHvxwavg=0.0
      SHvxlavg2=0.0
      SHvxwavg2=0.0
      SHvxlwavg2=0.0
c volume-based ratio averages for non-spherical particles 
c (SH L/T>SHcutoff, nonSH L/T > nonSHcutoff)
      nSHvxlwavg=0.0
      nSHvxlavg=0.0
      nSHvxwavg=0.0
      nSHvxlavg2=0.0
      nSHvxwavg2=0.0
      nSHvxlwavg2=0.0

c #-based for all particles, L/T, W/T, L/W
      xlavg=0.0
      xlwavg=0.0
      xwavg=0.0
      xtavg=0.0
      voltot=0.0
c vol-based for all particles, LWT
      avgl=0.0
      avgw=0.0
      avgt=0.0
      avgV=0.0
      avgl2=0.0
      avgw2=0.0
      avgt2=0.0
      avgV2=0.0
c single particles (SH L/T < SHcutoff, nonSH L/T < nonSHcutoff)
c V = VESD
      SHavgl=0.0
      SHavgw=0.0
      SHavgt=0.0
      SHavgV=0.0
      SHavgl2=0.0
      SHavgw2=0.0
      SHavgt2=0.0
      SHavgV2=0.0
c non-spherical particles (SH L/T > SHcutoff, nonSH > nonSHcutoff)
c V = VESD
      nonSHavgl=0.0
      nonSHavgw=0.0
      nonSHavgt=0.0
      nonSHavgV=0.0
      nonSHavgl2=0.0
      nonSHavgw2=0.0
      nonSHavgt2=0.0
      nonSHavgV2=0.0

c find min and max diameter for all particles
      diammax=0.0
      diammin=10000.
c find min and max diameter for all non-spherical (double or more) particles
      dmaxmul=0.0
      dminmul=10000.
c find min and max diameter for all single particles
      dmaxsin=0.0
      dminsin=10000.

c all particles count
      ic=0
c single particle count
      icsin=0
c non-spherical particle count
      icmult=0
c number of nonSH particles that are considered single
      NONsin=0

c do non-SH particles first
      do 109 i=1,N1
      read(19,*) filename,volumei,xli,xwi,xti
      vd3di=(6.*volumei/pi)**(1./3.)
c adjust L>W.T (sometimes small difference the wrong way)
       if(xwi.lt.xti) then
        write(8,*) filename,' W < T'
        write(8,111) xli,xwi,xti,(xwi-xti)/0.5/(xti+xwi)
        yyy=max(xti,xwi)
       xti=yyy
       xwi=yyy
       write(8,111) xwi,xti
       end if
       if(xli.lt.xwi) then
        write(8,119) filename
119     format(a30,' L < W')
        write(8,111) xli,xwi,xti,(xli-xwi)/0.5/(xli+xwi)
111    format(4f12.6)
        yyy=max(xli,xwi)
       xli=yyy
       xwi=yyy
       write(8,111) xli,xwi
       end if

       ic=ic+1
      volume(ic)=volumei
      xl(ic)=xli
      xw(ic)=xwi
      xt(ic)=xti
      vd3d(ic)=vd3di

      voltot=voltot+volume(ic)
      avgl=avgl+xl(ic)*volume(ic)
      avgw=avgw+xw(ic)*volume(ic)
      avgt=avgt+xt(ic)*volume(ic)
      avgV=avgV+vd3d(ic)*volume(ic)
      avgl2=avgl2+(xl(ic)**2)*volume(ic)
      avgw2=avgw2+(xw(ic)**2)*volume(ic)
      avgt2=avgt2+(xt(ic)**2)*volume(ic)
      avgV2=avgV2+(vd3d(ic)**2)*volume(ic)
      xlavg=xlavg+xl(ic)/xt(ic)
      xwavg=xwavg+xw(ic)/xt(ic)
      xlwavg=xlwavg+xl(ic)/xw(ic)

      vxlavg=vxlavg+volume(ic)*xl(ic)/xt(ic)
      vxwavg=vxwavg+volume(ic)*xw(ic)/xt(ic)
      vxlwavg=vxlwavg+volume(ic)*xl(ic)/xw(ic)
      vxlavg2=vxlavg2+volume(ic)*(xl(ic)/xt(ic))**2
      vxwavg2=vxwavg2+volume(ic)*(xw(ic)/xt(ic))**2
      vxlwavg2=vxlwavg2+volume(ic)*(xl(ic)/xw(ic))**2

      if(mmm.eq.1)  diam(ic)=vd3d(ic)
      if(mmm.eq.2)  diam(ic)=xl(ic)
      if(mmm.eq.3)  diam(ic)=xw(ic)
      if(mmm.eq.4)  diam(ic)=xt(ic)
      if(mmm.eq.5)  diam(ic)=xl(ic)/xw(ic)
      if(mmm.eq.6)  diam(ic)=xw(ic)/xt(ic)
      if(mmm.eq.7)  diam(ic)=xl(ic)/xt(ic)

c non-spherical particles
      if((xl(ic)/xt(ic)).ge.nSHcut) then
       icmult=icmult+1
       diammul(icmult)=diam(ic)
       volmul(icmult)=volume(ic)
       vtmul=vtmul+volmul(icmult)
      nonSHavgl=nonSHavgl+xl(ic)*volmul(icmult)
      nonSHavgw=nonSHavgw+xw(ic)*volmul(icmult)
      nonSHavgt=nonSHavgt+xt(ic)*volmul(icmult)
      nonSHavgV=nonSHavgV+vd3d(ic)*volmul(icmult)
      nonSHavgl2=nonSHavgl2+(xl(ic)**2)*volmul(icmult)
      nonSHavgw2=nonSHavgw2+(xw(ic)**2)*volmul(icmult)
      nonSHavgt2=nonSHavgt2+(xt(ic)**2)*volmul(icmult)
      nonSHavgV2=nonSHavgV2+(vd3d(ic)**2)*volmul(icmult)
      nSHvxlavg=nSHvxlavg+volmul(icmult)*xl(ic)/xt(ic)
      nSHvxwavg=nSHvxwavg+volmul(icmult)*xw(ic)/xt(ic)
      nSHvxlwavg=nSHvxlwavg+volmul(icmult)*xl(ic)/xw(ic)
      nSHvxlavg2=nSHvxlavg2+volmul(icmult)*(xl(ic)/xt(ic))**2
      nSHvxwavg2=nSHvxwavg2+volmul(icmult)*(xw(ic)/xt(ic))**2
      nSHvxlwavg2=nSHvxlwavg2+volmul(icmult)*(xl(ic)/xw(ic))**2
      if(dmaxmul.lt.diammul(icmult)) dmaxmul=diammul(icmult)
      if(dminmul.gt.diammul(icmult)) dminmul=diammul(icmult)
      end if
c single particles
      if((xl(ic)/xt(ic)).lt.nSHcut) then
      NONsin=NONsin+1
       icsin=icsin+1
       diamsin(icsin)=diam(ic)
       volsin(icsin)=volume(ic)
       vtsin=vtsin+volsin(icsin)
      SHavgl=SHavgl+xl(ic)*volsin(icsin)
      SHavgw=SHavgw+xw(ic)*volsin(icsin)
      SHavgt=SHavgt+xt(ic)*volsin(icsin)
      SHavgV=SHavgV+vd3d(ic)*volsin(icsin)
      SHavgl2=SHavgl2+(xl(ic)**2)*volsin(icsin)
      SHavgw2=SHavgw2+(xw(ic)**2)*volsin(icsin)
      SHavgt2=SHavgt2+(xt(ic)**2)*volsin(icsin)
      SHavgV2=SHavgV2+(vd3d(ic)**2)*volsin(icsin)
      SHvxlavg=SHvxlavg+volsin(icsin)*xl(ic)/xt(ic)
      SHvxwavg=SHvxwavg+volsin(icsin)*xw(ic)/xt(ic)
      SHvxlwavg=SHvxlwavg+volsin(icsin)*xl(ic)/xw(ic)
      SHvxlavg2=SHvxlavg2+volsin(icsin)*(xl(ic)/xt(ic))**2
      SHvxwavg2=SHvxwavg2+volsin(icsin)*(xw(ic)/xt(ic))**2
      SHvxlwavg2=SHvxlwavg2+volsin(icsin)*(xl(ic)/xw(ic))**2
      if(dmaxsin.lt.diamsin(icsin)) dmaxsin=diamsin(icsin)
      if(dminsin.gt.diamsin(icsin)) dminsin=diamsin(icsin)
      end if

      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

c now do particles that passed into geom file
      numSH=0
c num of SH particles with L/T > SHcut
      numLT=0
      do 100 i=1,NS
      read(9,*) filename,x1,x2,y1,y2,z1,z2,volumei,sa3di,sar,
     &           vd3di,trIr,nnn,gaussI,curvi,xli,xwi,xti

c adjust L>W>T (sometimes small difference the wrong way)
       if(xwi.lt.xti) then
        write(8,*) filename,' W < T'
        write(8,111) xli,xwi,xti,(xwi-xti)/0.5/(xti+xwi)
        yyy=max(xti,xwi)
       xti=yyy
       xwi=yyy
       write(8,111) xwi,xti
       end if
       if(xli.lt.xwi) then
        write(8,119) filename
        write(8,111) xli,xwi,xti,(xli-xwi)/0.5/(xli+xwi)
        yyy=max(xli,xwi)
       xli=yyy
       xwi=yyy
       write(8,111) xli,xwi
       end if

       ic=ic+1
       numSH=numSH+1
      volume(ic)=volumei
      volnumavg=volnumavg+volume(ic)
      xl(ic)=xli
      xw(ic)=xwi
      xt(ic)=xti
      vd3d(ic)=vd3di

      if(mmm.eq.1)  diam(ic)=vd3d(ic)
      if(mmm.eq.2)  diam(ic)=xl(ic)
      if(mmm.eq.3)  diam(ic)=xw(ic)
      if(mmm.eq.4)  diam(ic)=xt(ic)
      if(mmm.eq.5)  diam(ic)=xl(ic)/xw(ic)
      if(mmm.eq.6)  diam(ic)=xw(ic)/xt(ic)
      if(mmm.eq.7)  diam(ic)=xl(ic)/xt(ic)
c record all single and non-spherical SH particles

c non-spherical particles
      if((xl(ic)/xt(ic)).ge.SHcut) then
       numLT=numLT+1
       icmult=icmult+1
       diammul(icmult)=diam(ic)
       volmul(icmult)=volume(ic)
       vtmul=vtmul+volmul(icmult)
      nonSHavgl=nonSHavgl+xl(ic)*volmul(icmult)
      nonSHavgw=nonSHavgw+xw(ic)*volmul(icmult)
      nonSHavgt=nonSHavgt+xt(ic)*volmul(icmult)
      nonSHavgV=nonSHavgV+vd3d(ic)*volmul(icmult)
      nonSHavgl2=nonSHavgl2+(xl(ic)**2)*volmul(icmult)
      nonSHavgw2=nonSHavgw2+(xw(ic)**2)*volmul(icmult)
      nonSHavgt2=nonSHavgt2+(xt(ic)**2)*volmul(icmult)
      nonSHavgV2=nonSHavgV2+(vd3d(ic)**2)*volmul(icmult)
      nSHvxlavg=nSHvxlavg+volmul(icmult)*xl(ic)/xt(ic)
      nSHvxwavg=nSHvxwavg+volmul(icmult)*xw(ic)/xt(ic)
      nSHvxlwavg=nSHvxlwavg+volmul(icmult)*xl(ic)/xw(ic)
      nSHvxlavg2=nSHvxlavg2+volmul(icmult)*(xl(ic)/xt(ic))**2
      nSHvxwavg2=nSHvxwavg2+volmul(icmult)*(xw(ic)/xt(ic))**2
      nSHvxlwavg2=nSHvxlwavg2+volmul(icmult)*(xl(ic)/xw(ic))**2
      if(dmaxmul.lt.diammul(icmult)) dmaxmul=diammul(icmult)
      if(dminmul.gt.diammul(icmult)) dminmul=diammul(icmult)
      end if
c single particles
      if((xl(ic)/xt(ic)).lt.SHcut) then
       icsin=icsin+1
       diamsin(icsin)=diam(ic)
       volsin(icsin)=volume(ic)
       vtsin=vtsin+volsin(icsin)
      SHavgl=SHavgl+xl(ic)*volsin(icsin)
      SHavgw=SHavgw+xw(ic)*volsin(icsin)
      SHavgt=SHavgt+xt(ic)*volsin(icsin)
      SHavgV=SHavgV+vd3d(ic)*volsin(icsin)
      SHavgl2=SHavgl2+(xl(ic)**2)*volsin(icsin)
      SHavgw2=SHavgw2+(xw(ic)**2)*volsin(icsin)
      SHavgt2=SHavgt2+(xt(ic)**2)*volsin(icsin)
      SHavgV2=SHavgV2+(vd3d(ic)**2)*volsin(icsin)
      SHvxlavg=SHvxlavg+volsin(icsin)*xl(ic)/xt(ic)
      SHvxwavg=SHvxwavg+volsin(icsin)*xw(ic)/xt(ic)
      SHvxlwavg=SHvxlwavg+volsin(icsin)*xl(ic)/xw(ic)
      SHvxlavg2=SHvxlavg2+volsin(icsin)*(xl(ic)/xt(ic))**2
      SHvxwavg2=SHvxwavg2+volsin(icsin)*(xw(ic)/xt(ic))**2
      SHvxlwavg2=SHvxlwavg2+volsin(icsin)*(xl(ic)/xw(ic))**2
      if(dmaxsin.lt.diamsin(icsin)) dmaxsin=diamsin(icsin)
      if(dminsin.gt.diamsin(icsin)) dminsin=diamsin(icsin)
      end if

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

      voltot=voltot+volume(ic)
      avgl=avgl+xl(ic)*volume(ic)
      avgw=avgw+xw(ic)*volume(ic)
      avgt=avgt+xt(ic)*volume(ic)
      avgV=avgV+vd3d(ic)*volume(ic)
      avgl2=avgl2+(xl(ic)**2)*volume(ic)
      avgw2=avgw2+(xw(ic)**2)*volume(ic)
      avgt2=avgt2+(xt(ic)**2)*volume(ic)
      avgV2=avgV2+(vd3d(ic)**2)*volume(ic)
      xlavg=xlavg+xl(ic)/xt(ic)
      xwavg=xwavg+xw(ic)/xt(ic)
      xlwavg=xlwavg+xl(ic)/xw(ic)
      vxlavg=vxlavg+volume(ic)*xl(ic)/xt(ic)
      vxwavg=vxwavg+volume(ic)*xw(ic)/xt(ic)
      vxlwavg=vxlwavg+volume(ic)*xl(ic)/xw(ic)
      vxlavg2=vxlavg2+volume(ic)*(xl(ic)/xt(ic))**2
      vxwavg2=vxwavg2+volume(ic)*(xw(ic)/xt(ic))**2
      vxlwavg2=vxlwavg2+volume(ic)*(xl(ic)/xw(ic))**2

100   continue
      close(9)
      write(7,*) ' # SH parts = ',NS
      write(7,*) ' # SH parts with L/T =>',SHcut,numLT
      write(7,*) ' # SH parts with L/T <',SHcut,NS-numLT
      write(7,*) ' # nonSH parts = ',N1
      write(7,*) ' # nonSH parts with L/T =>',nSHcut,' ',N1-NONsin
      write(7,*) ' # nonSH parts w/ L/T <',nSHcut,' ',NONsin
      write(7,*) ' total # of particles ',ic
      write(7,*) ' # single particles ',icsin
      write(7,*) ' # non-spherical particles ',icmult
      write(7,*)
      write(7,*) ' volume of single particles = ',vtsin
      write(7,*) ' volume of non-spherical particles = ',vtmul
      write(7,*) ' total powder volume = ',voltot/1.0e+9,' mm**3'
      write(7,*) ' total powder volume = ',(vtsin+vtmul)/1.0e+9,' mm**3'
      write(7,*) ' % single-particle volume = ',100.*vtsin/voltot
      write(7,*) ' % non-spherical-particle volume = ',100.*vtmul/voltot
      write(7,*) 
      N=ic
      xlwavg=xlwavg/float(N)
      xlavg=xlavg/float(N)
      xwavg=xwavg/float(N)
      vxlwavg=vxlwavg/voltot
      vxlavg=vxlavg/voltot
      vxwavg=vxwavg/voltot
      vxlwavg2=vxlwavg2/voltot
      vxlavg2=vxlavg2/voltot
      vxwavg2=vxwavg2/voltot
      vxlwavg2=sqrt(vxlwavg2-vxlwavg**2)
      vxlavg2=sqrt(vxlavg2-vxlavg**2)
      vxwavg2=sqrt(vxwavg2-vxwavg**2)
c ratio volume-based averages for single particles
      SHvxlavg=SHvxlavg/vtsin
      SHvxwavg=SHvxwavg/vtsin
      SHvxlwavg=SHvxlwavg/vtsin
      SHvxlavg2=SHvxlavg2/vtsin
      SHvxwavg2=SHvxwavg2/vtsin
      SHvxlwavg2=SHvxlwavg2/vtsin
      SHvxlwavg2=sqrt(SHvxlwavg2-SHvxlwavg**2)
      SHvxlavg2=sqrt(SHvxlavg2-SHvxlavg**2)
      SHvxwavg2=sqrt(SHvxwavg2-SHvxwavg**2)
c ratio volume-based averages for non-spherical particles
      nSHvxlavg=nSHvxlavg/vtmul
      nSHvxwavg=nSHvxwavg/vtmul
      nSHvxlwavg=nSHvxlwavg/vtmul
      nSHvxlavg2=nSHvxlavg2/vtmul
      nSHvxwavg2=nSHvxwavg2/vtmul
      nSHvxlwavg2=nSHvxlwavg2/vtmul
      nSHvxlwavg2=sqrt(nSHvxlwavg2-nSHvxlwavg**2)
      nSHvxlavg2=sqrt(nSHvxlavg2-nSHvxlavg**2)
      nSHvxwavg2=sqrt(nSHvxwavg2-nSHvxwavg**2)

      avgl=avgl/voltot
      avgw=avgw/voltot
      avgt=avgt/voltot
      avgV=avgV/voltot
      avgl2=avgl2/voltot
      avgw2=avgw2/voltot
      avgt2=avgt2/voltot
      avgV2=avgV2/voltot
      avgl2=sqrt(avgl2-avgl**2)
      avgw2=sqrt(avgw2-avgw**2)
      avgt2=sqrt(avgt2-avgt**2)
      avgV2=sqrt(avgV2-avgV**2)

      volnumavg=volnumavg/float(numSH)

c LWT volume-based averages for single particles
      SHavgl=SHavgl/vtsin
      SHavgw=SHavgw/vtsin
      SHavgt=SHavgt/vtsin
      SHavgV=SHavgV/vtsin
      SHavgl2=SHavgl2/vtsin
      SHavgw2=SHavgw2/vtsin
      SHavgt2=SHavgt2/vtsin
      SHavgV2=SHavgV2/vtsin
      SHavgl2=sqrt(SHavgl2-SHavgl**2)
      SHavgw2=sqrt(SHavgw2-SHavgw**2)
      SHavgt2=sqrt(SHavgt2-SHavgt**2)
      SHavgV2=sqrt(SHavgV2-SHavgV**2)
c LWT volume-based averages for non-spherical particles
      nonSHavgl=nonSHavgl/vtmul
      nonSHavgw=nonSHavgw/vtmul
      nonSHavgt=nonSHavgt/vtmul
      nonSHavgV=nonSHavgV/vtmul
      nonSHavgl2=nonSHavgl2/vtmul
      nonSHavgw2=nonSHavgw2/vtmul
      nonSHavgt2=nonSHavgt2/vtmul
      nonSHavgV2=nonSHavgV2/vtmul
      nonSHavgl2=sqrt(nonSHavgl2-nonSHavgl**2)
      nonSHavgw2=sqrt(nonSHavgw2-nonSHavgw**2)
      nonSHavgt2=sqrt(nonSHavgt2-nonSHavgt**2)
      nonSHavgV2=sqrt(nonSHavgV2-nonSHavgV**2)

      write(7,*) ' min max ',letter,' for all particles'
      write(7,*) diammin,diammax
      write(7,*) ' min max ',letter,' for single particles'
      write(7,*) dminsin,dmaxsin
      write(7,*) ' min max ',letter,' for non-spherical particles'
      write(7,*) dminmul,dmaxmul
      write(7,*) 
      write(7,*) ' vol averages of L/T W/T L/W for all particles'
      write(7,*) vxlavg,vxwavg,vxlwavg
      write(7,*) ' vol std of L/T W/T L/W for all particles'
      write(7,*) vxlavg2,vxwavg2,vxlwavg2
      write(7,*) 
      write(7,*) ' vol averages of L/T W/T L/W for single particles'
      write(7,*) SHvxlavg,SHvxwavg,SHvxlwavg
      write(7,*) ' vol std of L/T W/T L/W for single particles'
      write(7,*) SHvxlavg2,SHvxwavg2,SHvxlwavg2
      write(7,*) 
      write(7,*) ' vol avgs L/T W/T L/W for non-spherical particles'
      write(7,*) nSHvxlavg,nSHvxwavg,nSHvxlwavg
      write(7,*) ' vol std of L/T W/T L/W for non-spherical particles'
      write(7,*) nSHvxlavg2,nSHvxwavg2,nSHvxlwavg2
      write(7,*) 
      write(7,*) ' vol averages of LWT for all particles'
      write(7,*) avgl,avgw,avgt
      write(7,*) ' vol std of LWT for all particles'
      write(7,*) avgl2,avgw2,avgt2
      write(7,*) ' vol std of LWT / avg for all particles'
      write(7,*) avgl2/avgl,avgw2/avgw,avgt2/avgt
      write(7,*) 
      write(7,*) ' vol averages of LWT for single particles'
      write(7,*) SHavgl,SHavgw,SHavgt
      write(7,*) ' vol std of LWT for single particles'
      write(7,*) SHavgl2,SHavgw2,SHavgt2
      write(7,*) ' vol std of LWT / avg for single particles'
      write(7,*) SHavgl2/SHavgl,SHavgw2/SHavgw,SHavgt2/SHavgt
      write(7,*) 
      write(7,*) ' vol averages of LWT for non-spherical particles'
      write(7,*) nonSHavgl,nonSHavgw,nonSHavgt
      write(7,*) ' vol std of LWT for non-spherical particles'
      write(7,*) nonSHavgl2,nonSHavgw2,nonSHavgt2
      write(7,*) ' vol std of LWT / avg for non-spherical particles'
      write(7,*) nonSHavgl2/nonSHavgl,nonSHavgw2/nonSHavgw,
     $               nonSHavgt2/nonSHavgt

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

c note that psd(i,1) is vol-based, psd(i,2) is num-based
c do single particles first
      if(mmm.eq.3) then
      dminsin=diammin
      dmaxsin=diammax
      end if
      write(7,*) icsin,' single particles'
      dx=(dmaxsin-dminsin)/float(nbins)
      write(7,*) nbins,dx
      do 350 i=1,1000
      psd(i,1)=0.0
      psd(i,2)=0.0
350   continue

c do differential probability distribution for diam choice

      do 500 i=1,icsin
        m=1+(diamsin(i)-dminsin)/dx
        if(m.gt.nbins) m=nbins
        psd(m,1)=psd(m,1)+volsin(i)/vtsin
        psd(m,2)=psd(m,2)+1.0
500   continue

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

      write(10,129) dminsin,0.0
129   format(2f11.6)
      do 460 i=1,nbins
      x=dminsin+(i-1)*dx
      write(10,129) x,psd(i,1)
      write(10,129) x+dx,psd(i,1)
      write(10,129) x+dx,0.0
460   continue

c do non-spherical particles next
      if(mmm.eq.3) then
      dminmul=diammin
      dmaxmul=diammax
      end if
      write(7,*) icmult,' non-spherical particles'
      dx=(dmaxmul-dminmul)/float(nbins)
      write(7,*) nbins,dx
      do 355 i=1,1000
      psd(i,1)=0.0
      psd(i,2)=0.0
355   continue

c do differential probability distribution for diam choice

      do 509 i=1,icmult
        m=1+(diammul(i)-dminmul)/dx
        if(m.gt.nbins) m=nbins
        psd(m,1)=psd(m,1)+volmul(i)/vtmul
        psd(m,2)=psd(m,2)+1.0
509   continue

      sum=0.0
      do 365 i=1,nbins
      sum=sum+psd(i,1)
      x=dminmul+(i-1)*dx
      write(7,113) x,x+dx,psd(i,1),psd(i,2),sum
365   continue
      write(7,*) sum

      write(11,129) dminmul,0.0
c129   format(2f11.6)
      do 465 i=1,nbins
      x=dminmul+(i-1)*dx
      write(11,129) x,psd(i,1)
      write(11,129) x+dx,psd(i,1)
      write(11,129) x+dx,0.0
465   continue

c do single+ non-spherical 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=135.
c      diammin=7.0
c      end if
c      if(mmm.eq.5) then
c      diammax=3.0
c      diammin=1.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),sum
366   continue
      write(7,*) sum

      write(12,129) diammin,0.0
      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(10)
      close(11)
      close(12)
5544  continue
      end