c pick out VRML files where L/T is in segments. First do all of them
c then pick random ones, maybe 10 per interval, to show that
c L/T = 1.4 is a good cutoff for single vs. multiple particles
c for the SH particles. All non-SH particles are assumed to be 
c multi-particles, so maybe pick some of those as well

c generate LWT distributions
c divide particles by single (SH, L/T < 1.4) and multiple 
c (SH, L/T > 1.4 and non-SH)
c later maybe will try dividing the multiple particle into
c double, triple, and so-on, using the knowledge I am getting from
c the multiple 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)

      integer nn(26),gaussI
      character*50 filename

      mmm=7

      if(mmm.eq.7) then
      open(unit=7,file='ASTM-AMPS2-LT-SH-VRML-list.txt')
      open(unit=8,file='SH-VRML-list.txt')
      open(unit=17,file='ASTM-AMPS2-LT-nonSH-VRML-list.txt')
      open(unit=18,file='nonSH-VRML-list.txt')
      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 SH
      N=14581
c non-SH
      N1=2390

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

c volume-based ratio averages for all particles
      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 particles
      SHvxlwavg=0.0
      SHvxlavg=0.0
      SHvxwavg=0.0
      SHvxlavg2=0.0
      SHvxwavg2=0.0
      SHvxlwavg2=0.0
c volume-based ratio averages for multiple particles (L/T<1.4 SH,nonSH)
      nSHvxlwavg=0.0
      nSHvxlavg=0.0
      nSHvxwavg=0.0
      nSHvxlavg2=0.0
      nSHvxwavg2=0.0
      nSHvxlwavg2=0.0

      xlavg=0.0
      xlwavg=0.0
      xwavg=0.0
      xtavg=0.0
      voltot=0.0
      avgl=0.0
      avgw=0.0
      avgt=0.0
      avgl2=0.0
      avgw2=0.0
      avgt2=0.0
c single particles (SH L/T < 1.4)
      SHavgl=0.0
      SHavgw=0.0
      SHavgt=0.0
      SHavgl2=0.0
      SHavgw2=0.0
      SHavgt2=0.0
c multiple particles (SH L/T > 1.4, all non-SH)
      nonSHavgl=0.0
      nonSHavgw=0.0
      nonSHavgt=0.0
      nonSHavgl2=0.0
      nonSHavgw2=0.0
      nonSHavgt2=0.0

c find min and max diameter for all particles
      diammax=0.0
      diammin=10000.
c find min and max diameter for all multiple (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 multiple particle count
      icmult=0

c do non-SH particles first
      do 1009 ijk=1,10
      cut1=1.0+(ijk-1)*0.1
      cut2=1.0+ijk*0.1
      ncut=0
      do 109 i=1,N1
      read(19,*) filename,volumei,xli,xwi,xti
      vd3di=(6.*volumei/pi)**(1./3.)
       if(xwi.lt.xti) then
        yyy=max(xti,xwi)
       xti=yyy
       xwi=yyy
       end if
       if(xli.lt.xwi) then
119     format(a30,' L < W')
111    format(4f12.6)
        yyy=max(xli,xwi)
       xli=yyy
       xwi=yyy
       end if

       ic=ic+1
      xl(ic)=xli
      xw(ic)=xwi
      xt(ic)=xti

      if(mmm.eq.7)  diam(ic)=xl(ic)/xt(ic)
      if(mmm.eq.7) then
      if(diam(ic).gt.cut1.and.diam(ic).le.cut2) then
       ncut=ncut+1
       if(ncut.le.10) write(17,*) filename,diam(ic)
       if(ncut.le.10) write(18,*) filename
      end if
      end if

109   continue
      rewind(19)
1009  continue

      do 1000 ijk=1,10
      cut1=1.0+(ijk-1)*0.1
      cut2=1.0+ijk*0.1
      ncut=0
      do 100 i=1,N
      read(9,*) filename,x1,x2,y1,y2,z1,z2,volumei,sa3di,sar,
     &           vd3di,trIr,nnn,gaussI,curvi,xli,xwi,xti

       if(xwi.lt.xti) then
        yyy=max(xti,xwi)
       xti=yyy
       xwi=yyy
       end if
       if(xli.lt.xwi) then
        yyy=max(xli,xwi)
       xli=yyy
       xwi=yyy
       end if

       ic=ic+1
      xl(ic)=xli
      xw(ic)=xwi
      xt(ic)=xti

      if(mmm.eq.7)  diam(ic)=xl(ic)/xt(ic)
      if(mmm.eq.7) then
      if(diam(ic).gt.cut1.and.diam(ic).le.cut2) then
       ncut=ncut+1
       if(ncut.le.10) write(7,*) filename,diam(ic)
       if(ncut.le.10) write(8,*) filename
      end if
      end if

100   continue
      rewind(9)
1000  continue
      end