c
c     closurestats - calculate statistics (area, volume, etc) for
c                    fundamental closures identified by sortclosures
c                    program.
c
c     Saltus - 14 Sep 2001 - for NPRA oil & gas assessment project
c
      program closurestats
c
      parameter (MAXPTS=300,MAXROWS=500,MAXCOLS=600)
      parameter (MAXCLOSE=5000)
      parameter (dval1=1e30,dval2=-1e31)
      dimension xpts(MAXCLOSE,MAXPTS),ypts(MAXCLOSE,MAXPTS)
      dimension grid(MAXROWS,MAXCOLS)
      dimension hgrid(MAXROWS,MAXCOLS)
      dimension areaka(MAXCLOSE),volkm3(MAXCLOSE)
      dimension topdepth(MAXCLOSE),ttopdepth(MAXCLOSE)
      dimension ntopdepth(MAXCLOSE)
      dimension heightft(MAXCLOSE),narea(MAXCLOSE),nheight(MAXCLOSE)
      dimension nareahist(MAXCLOSE/10),nheighthist(MAXCLOSE/10)
      dimension ntophist(MAXCLOSE/10)
      dimension tareaka(MAXCLOSE),theightft(MAXCLOSE)
      dimension thmin(MAXCLOSE),thmax(MAXCLOSE),tzval(MAXCLOSE)
      dimension txcenter(MAXCLOSE),tycenter(MAXCLOSE),tnpts(MAXCLOSE)
      dimension itcont(MAXCLOSE)
      character id*56,pgm*8
      character buffer*80
      integer tnpts
c
      print *,'open chosen.xyz, grid.sg, and stats.xyz'
c
      open(10,file="chosen.xyz",form='formatted',status='old')
      open(11,file="grid.sg",form='unformatted',status='old')
      open(12,file="stats.xyz",form='formatted',status='unknown')
      open(13,file="centers.xyz",form='formatted',status='unknown')
      open(14,file="ccontours.xyz",form='formatted',status='unknown')
c
c     read grid into memory
c
      print *,'reading grid into memory...'
      read(11)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      print *, id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      cellarea=dx*dy
      ndatacells=0
      do 3 i=1,nrow
        read(11)dummy,(grid(i,j),j=1,ncol)
        do 2 k=1,ncol
          if ((grid(i,k).gt.-1e30).and.(grid(i,k).lt.1e30)) then
            ndatacells=ndatacells+1
          end if
    2   continue
    3 continue
c
      print *, 'grid input to memory'
      percentdata=real(ndatacells)/real(ncol*nrow)
      dataareakm=ndatacells*cellarea/(1000*1000)
c
c     1 km^2 = 0.391 sq mile
c     1 sq mile = 640 acres
c     1 km^2 = 250 acres
c     (km^2)/4 = k acres = ka
c
      dataareaka=dataareakm/4
      print *, 'data area (% of grid) =',percentdata*100.
      print *, 'data area (kacres) =',dataareaka
c
      buffer='"Automated statistics on closed highs"'
      write(12,*)buffer
      buffer='"from program closurestats v1.0 (Saltus, USGS, 10/2001)"'
      write(12,*)buffer
      buffer='"CONTOUR INTERVAL:"'
      write(12,*)buffer
      buffer='"GRID INTERVAL:"'
      write(12,*)buffer
      buffer='"PLAY NAME:"'
      write(12,*)buffer
      buffer='"SURFACE NAME:"'
      write(12,*)buffer
      write(12,*)
      buffer='"Total data area of grid (kacres)"'
      write(12,*)buffer
      write(12,*)dataareaka
      write(12,*)
c
c     initial hgrid (output of closure heights)
c
      do 4 i=1,nrow
      do 44 j=1,ncol
        hgrid(i,j)=dval2
   44 continue
    4 continue
c
c     input contours and compute stats
c
      ncount=1
      ncontour=1
    5 continue
      read(10,*,end=999)xpts(ncontour,1),ypts(ncontour,1),
     & zval,ndir,num,tnpts(ncontour)
      itcont(ncount)=0
      xmin=xpts(ncontour,1)
      xmax=xpts(ncontour,1)
      ymin=ypts(ncontour,1)
      ymax=ypts(ncontour,1)
      do 10 i=2,tnpts(ncontour)
      read(10,*)xpts(ncontour,i),ypts(ncontour,i)
      if (xpts(ncontour,i).lt.xmin) xmin=xpts(ncontour,i)
      if (xpts(ncontour,i).gt.xmax) xmax=xpts(ncontour,i)
      if (ypts(ncontour,i).lt.ymin) ymin=ypts(ncontour,i)
      if (ypts(ncontour,i).gt.ymax) ymax=ypts(ncontour,i)
   10 continue
c     print *,'read',num,' contour with',tnpts(ncontour),' points.'
c     print *,xmin,xmax,ymin,ymax,' = xmin,xmax,ymin,ymax'
c
c     calculate statistics
c
      istart=1+(ymin-yo)/dy
      istart=max(1,istart)
      istart=min(nrow,istart)
      iend=1+(ymax-yo)/dy
      iend=max(1,iend)
      iend=min(nrow,iend)
      jstart=1+(xmin-xo)/dx
      jstart=max(1,jstart)
      jstart=min(ncol,jstart)
      jend=1+(xmax-xo)/dx
      jend=max(1,jend)
      jend=min(ncol,jend)
c     if ((istart.eq.ncol).or.(jstart.eq.nrow)) then
c       print *,'Outside grid area',istart,iend,jstart,jend
c       goto 5
c     end if

c
c     check that grid coords are within spec
c
c
c     print *,istart,iend,jstart,jend,
c    &    ' = istart,iend,jstart,jend'
c
      xcenter=xmin+(xmax-xmin)/2
      ycenter=ymin+(ymax-ymin)/2
c
      area=0.
      volpos=0.
      volneg=0.
      heightmax=0.
      heightmin=0.
      ncellp=0.
      ncelln=0.
c
      do 20 i=istart,iend
        do 15 j=jstart,jend
          xtest=xo+(j-1)*dx
          ytest=yo+(i-1)*dy
c     print *,'calling inside'
          intest=inside(xtest,ytest,xpts,ypts,tnpts(ncontour),
     &                  ncontour,MAXCLOSE,MAXPTS)
c     print *,'back from inside'
          if ((intest.eq.0).and.(grid(i,j).lt.1e20)
     &         .and.(grid(i,j).gt.-1e20)) then 
            cellheight=grid(i,j)-zval
            if (cellheight.gt.10.) then
              print *,ncontour,i,j,grid(i,j),zval
            end if
            if (cellheight.gt.0) then
              hgrid(i,j)=cellheight
              ncellp=ncellp+1
              volpos=volpos+cellheight*cellarea
              if (cellheight.gt.heightmax) heightmax=cellheight
            else
              ncelln=ncelln+1
              volneg=volneg+cellheight*cellarea
              if (cellheight.lt.heightmin) heightmin=cellheight
            end if
          end if
   15   continue
   20 continue
      volposkm3=volpos/(1000*1000*3.28)
      volnegkm3=volneg/(1000*1000*3.28)
      volume=volposkm3+volnegkm3
      areakm2=ncellp*cellarea/(1000*1000)
c
c     units - currently a mix of kft (zval) & meters (x,y)
c      
      if (volume.gt.0.0) then
        areaka(ncount)=areakm2/4.
        tareaka(ncount)=areakm2/4.
        volkm3(ncount)=volposkm3
        heightft(ncount)=3280.*(volposkm3/areakm2)
        theightft(ncount)=3280.*(volposkm3/areakm2)
        topdepth(ncount)=zval+heightmax
        ttopdepth(ncount)=topdepth(ncount)
c       print *,ncount,xcenter/1000,ycenter/1000,zval,areaka(ncount)
c       print *,ncontour
c       write(12,100)ncount,xcenter/1000,ycenter/1000,zval,areakm2/4,
c    &      volposkm3,heightft(ncount),heightmin*1000,heightmax*1000
        txcenter(ncount)=xcenter
        tycenter(ncount)=ycenter
        tzval(ncount)=zval
        thmin(ncount)=heightmin
        thmax(ncount)=heightmax
        itcont(ncount)=ncontour
        ncount=ncount+1
      end if
      ncontour=ncontour+1
      goto 5
c
  999 continue
      ncount=ncount-1
      print *,ncount,' closed contours with positive volume'
      print *,'creating hgrid.sg output grid...'
c
      open(15,file="hgrid.sg",form='unformatted',status='unknown')
      id="Closure height grid"
      pgm="clsrstat"
      write(15)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      do 50 i=1,nrow
        write(15)dummy,(hgrid(i,j),j=1,ncol)
   50 continue
c
c     find 5%, 50%, and 95% values for area & height
c
      i5=nint(real(ncount)*.05)
      i50=nint(real(ncount)*.5)
      i95=nint(real(ncount)*.95)
      if (i5.le.0) i5=1
      if (i95.gt.ncount) i95=ncount
      print *,'i5, i50, i95 =',i5,i50,i95
c
      call rshsort(areaka,narea,ncount)
      call rshsort(heightft,nheight,ncount)
      call rshsort(topdepth,ntopdepth,ncount)
      write(12,*)
      buffer='"AREA distribution (kacre)"'
      write(12,*)buffer
      write(12,150)'M5%','M50%','M95%'
      write(12,155)areaka(i5),areaka(i50),areaka(i95)
      print *,'Area (kacres) distribution:'
      print *,'M5%=',areaka(i5),' M50%=',areaka(i50),
     & ' M95%=',areaka(i95)
      write(12,*)
      buffer='"HEIGHT distribution (ft)"'
      write(12,*)buffer
      write(12,150)'M5%','M50%','M95%'
      write(12,155)heightft(i5),heightft(i50),heightft(i95)
      print *,'Height (ft) distribution:'
      print *,'M5%=',heightft(i5),' M50%=',heightft(i50),
     & ' M95%=',heightft(i95)
      write(12,*)
      buffer='"TOPDEPTH distribution (kft)"'
      write(12,*)buffer
      write(12,150)'M5%','M50%','M95%'
      write(12,155)topdepth(i5),topdepth(i50),topdepth(i95)
      print *,'Topdepth (kft) distribution:'
      print *,'M5%=',topdepth(i5),' M50%=',topdepth(i50),
     & ' M95%=',topdepth(i95)
  150 format(3a12)
  155 format(3f12.5) 
c
c     construct area and height histograms
c
      print *,'Histograms...'
c
      nbins=ncount/10.
      areaint=areaka(ncount)/nbins
      heightint=heightft(ncount)/nbins
      topint=(topdepth(ncount)-topdepth(1))/nbins
c
      print *,'area interval =',areaint
      print *,'height interval =',heightint
      print *,'top interval =',topint
c
      do 60 i=1,nbins
        nareahist(i)=0
        nheighthist(i)=0 
        ntophist(i)=0 
   60 continue
c
      do 65 i=1,ncount
        iabin=ifix(areaka(i)/areaint)+1
        if (iabin.gt.nbins) iabin=nbins
        ihbin=ifix(heightft(i)/heightint)+1
        if (ihbin.gt.nbins) ihbin=nbins
        itbin=ifix((topdepth(i)-topdepth(1))/topint)+1
        if (itbin.gt.nbins) itbin=nbins
        nareahist(iabin)=nareahist(iabin)+1
        nheighthist(ihbin)=nheighthist(ihbin)+1
        ntophist(itbin)=ntophist(itbin)+1
c       print *,areaka(i),iabin,heightft(i),ihbin
   65 continue
c
      write(12,*)
      buffer='"Histograms"'
      write(12,*)buffer
      write(12,160)'Bin','Area_ka','#pts','Height_ft','#pts',
     & 'Top_kft','#pts'
      do 70 i=1,nbins
        write(12,165)i,areaint*i,nareahist(i),
     &                 heightint*i,nheighthist(i),
     &                 topdepth(1)+topint*(i-1),ntophist(i)
  160   format(a5,a12,a5,a12,a5,a12,a5)
  165   format(i5,f12.2,i5,f12.0,i5,f12.3,i5)
   70 continue  
c
c     write stats file header
c
      write(12,*)
      buffer='"Summary table for closed highs"'
      write(12,*)buffer
      write(12,105)'num',' xcenter_km',' ycenter_km',
     & ' bot_kft','top_kft',' area_ka',' vol_km3',' vol_kaft',
     & ' have_ft',' hmin_ft',' hmax_ft'
  100   format(i5,2f12.3,2f8.1,f12.2,2f12.2,f8.1,2f10.1)
  105   format(a5,2a12,2a8,3a12,a8,2a10)
c
c     write stats file table
c
      do 80 i=1,ncount
        j=narea(ncount+1-i)
        write(12,100)i,txcenter(j)/1000,tycenter(j)/1000,
     &              tzval(j),ttopdepth(j),tareaka(j),
     &      volkm3(j),volkm3(j)*820,
     &      theightft(j),thmin(j)*1000,thmax(j)*1000
        write(13,*)txcenter(j),tycenter(j),i,tzval(j)
        write(14,'("LINE ",i5)')i
        do 30 ix=1,tnpts(itcont(j))
          write(14,*)xpts(itcont(j),ix),ypts(itcont(j),ix)
   30   continue
   80 continue
c
      close(11)
      close(12)
      close(13)
      close(14)
      close(15)
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine rshsort (array,neword,num)
c     Modified from Kernighan and Pflauger, Elements of Programming Style,
c       page 134.

c      array = list of reals; List is returned sorted...
c      neword = integer array
c         (If i is index in sorted list,
c          then neword(i) is position in original list.)
c      num  = number of items in list

      real array(num)
      integer neword(num)

      do 2 i=1,num
  2   neword(i)=i

      igap=num
  5   if(igap.le.1) return
      igap=igap/2
      imax=num-igap

 10   iex=0

      do 20 i=1,imax
      iplusg=i+igap
      if(array(i).le.array(iplusg)) go to 20
        save=array(i)
        array(i)=array(iplusg)
        array(iplusg)=save
        indtmp=neword(i)
        neword(i)=neword(iplusg)
        neword(iplusg)=indtmp
      iex=1
 20   continue

      if(iex.ne.0) goto 10
      goto 5

      end
C+
C      INSIDE - tests whether point is inside polygon
C
C    Call INSIDE (x,y,xbod,ybod,npt,nbd,mbod,mcorn)
C
C    x,y= point to test
C    xbod,ybod= two-dimensional arrays with corners of body
C    npt= number of points in body
C    nbd= body number (first subscript of xbod, ybod)
C    mbod= 1st dimension of xbod,ybod arrays
C    mcorn= 2nd dimension of xbod,ybod arrays
C
C    RETURNS:  0= x,y inside
C              1= x,y outside
C-
C  *********************************************
      Integer Function inside(x,y,xbod,ybod,npt,nbd,mbod,mcorn)
      Dimension xbod(mbod,mcorn),ybod(mbod,mcorn)
c     Double Precision 
C
C     bail for trivial cases npt=1 or 2
C
      if (npt.le.2) then
        inside=1
        return
      end if
C
C     find xmin,xmax,ymin,ymax for polygon
C
      xmin=xbod(nbd,1)
      xmax=xbod(nbd,1)
      ymin=ybod(nbd,1)
      ymax=ybod(nbd,1)
      do 10 i=1,npt
        if (xbod(nbd,i).lt.xmin) xmin=xbod(nbd,i)
        if (xbod(nbd,i).gt.xmax) xmax=xbod(nbd,i)
        if (ybod(nbd,i).lt.ymin) ymin=ybod(nbd,i)
        if (ybod(nbd,i).gt.ymax) ymax=ybod(nbd,i)
   10 continue
c     print *,'xmin,xmax,ymin,ymax=',xmin,xmax,ymin,ymax
C
C     check if test point completely outside polygon box
C
      if ((x.gt.xmax).or.(x.lt.xmin).or.(y.gt.ymax)
     &               .or.(y.lt.ymin)) then
        inside=1
c       print *,'test point outside polygon bounding box'
        return
      end if
C
C     construct radial line from test point to a point
C      outside the polygon min, max box
C
      x2=xmax+2.*abs(xmax)
      y2=ymax+3.*abs(ymax)
c
c
c     print *,'x2,y2=',x2,y2
C
C     check for intersections of radial line with all polygon sides
C     int is number of intersections
C     if int=odd, then pt is inside
C     if int=even, then pt is outside
C
      int=0
      ipar=0
      Do 50 i=1,npt-1
      xp1=xbod(nbd,i)
      yp1=ybod(nbd,i)
      xp2=xbod(nbd,i+1)
      yp2=ybod(nbd,i+1)
      intest=insect2(x,y,x2,y2,xp1,yp1,xp2,yp2,int,ipar)
c     print *,'i,int=',i,int
      if (intest.ne.0) then
        x2t=xmin-3.*abs(xmin)
        y2t=ymin-2.*abs(ymin)
        intest=insect2(x,y,x2t,y2t,xp1,yp1,xp2,yp2,int,ipar)
        if (intest.ne.0) then
          print *,'Possible missed intersection'
        end if
      end if
   50 Continue
      xp1=xbod(nbd,npt)
      yp1=ybod(nbd,npt)
      xp2=xbod(nbd,1)
      yp2=ybod(nbd,1)
      intest=insect2(x,y,x2,y2,xp1,yp1,xp2,yp2,int,ipar)
c     print *,'i,int=',i+1,int
      if (intest.ne.0) then
        x2t=xmin-3.*abs(xmin)
        y2t=ymin-2.*abs(ymin)
        intest=insect2(x,y,x2t,y2t,xp1,yp1,xp2,yp2,int,ipar)
        if (intest.ne.0) then
          print *,'Possible missed intersection'
        end if
      end if
C
C     examine int and report results
C
c     print *,'int=',int
      If (mod(int,2).EQ.0) then
        inside=1
c       int1=int
c
c     if first answer is outside, try again (just in case)
c
c       x2=xmin-2.*abs(xmin)
c       y2=ymin-3.*abs(ymin)
c       int=0
c       ipar=0
c       Do 150 i=1,npt-1
c       xp1=xbod(nbd,i)
c       yp1=ybod(nbd,i)
c       xp2=xbod(nbd,i+1)
c       yp2=ybod(nbd,i+1)
c       intest=insect2(x,y,x2,y2,xp1,yp1,xp2,yp2,int,ipar)
c       print *,'i,int=',i,int
c       if (intest.ne.0) then
c         x2t=xmin-3.*abs(xmin)
c         y2t=ymin-2.*abs(ymin)
c         intest=insect2(x,y,x2t,y2t,xp1,yp1,xp2,yp2,int,ipar)
c         if (intest.ne.0) then
c           print *,'Possible missed intersection'
c         end if
c       end if
c 150   Continue
c       xp1=xbod(nbd,npt)
c       yp1=ybod(nbd,npt)
c       xp2=xbod(nbd,1)
c       yp2=ybod(nbd,1)
c       intest=insect2(x,y,x2,y2,xp1,yp1,xp2,yp2,int,ipar)
c       print *,'i,int=',i+1,int
c       if (intest.ne.0) then
c         x2t=xmin-3.*abs(xmin)
c         y2t=ymin-2.*abs(ymin)
c         intest=insect2(x,y,x2t,y2t,xp1,yp1,xp2,yp2,int,ipar)
c         if (intest.ne.0) then
c           print *,'Possible missed intersection'
c         end if
c       end if
c       if (mod(int,2).ne.mod(int1,2)) then
c         print *,'Changed result with 2nd insect2 test...'
c       end if
c       if (mod(int,2).EQ.0) then
c         inside=1
c       else
c         inside=0
c       end if
      else
        inside=0
      end if
      Return
      End
C
C     INSECT2  finds intersection of two lines 
C
C  *********************************************
      Integer function insect2(x1,y1,x2,y2,a1,b1,a2,b2,int,ipar)
C
C     if an intersection is found, then increment int
C     if lines are parallel, increment ipar
C
C     insect2=0 normal return
C            =1 intersection falls on a1,b1
C            =2 intersection falls on a2,b2
C
C     conditions 1 & 2 could mean that line is tangent to
C       polygon (no intersection) or that it crosses at 
C       vertex (intersection).  Further testing is required
C       to solve this condition.
C
C     xsect(x1,y1,x2,y2,yloc)=(yloc-y1)*((x2-x1)/(y2-y1))+x1
C 
C     calculate line parameters
C
C     slopes (s) & intercepts (yint)
C
C     ltype = 0 for finite slope
C           = 1 for vertical line
C
      insect2=0
c     print *,'x1,y1,x2,y2=',x1,y1,x2,y2
c     print *,'a1,b1,a2,b2=',a1,b1,a2,b2
      if (x2.eq.x1) then
        ltype1=1
      else
        ltype1=0
        s1=(y2-y1)/(x2-x1)
        yint1=y1-s1*x1
c       print *,'s1,yint1=',s1,yint1
      end if
      if (a2.eq.a1) then
        ltype2=1
      else
        ltype2=0
        s2=(b2-b1)/(a2-a1)
        yint2=b1-s2*a1
c       print *,'s2,yint2=',s2,yint2
      end if
C
C     Set min, max for line segments
C
      if (x1.lt.x2) then
        xmin=x1
        xmax=x2
      else
        xmin=x2
        xmax=x1
      end if
      if (y1.lt.y2) then
        ymin=y1
        ymax=y2
      else
        ymin=y2
        ymax=y1
      end if
      if (a1.lt.a2) then
        amin=a1
        amax=a2
      else
        amin=a2
        amax=a1
      end if
      if (b1.lt.b2) then
        bmin=b1
        bmax=b2
      else
        bmin=b2
        bmax=b1
      end if
C
C     find intersection and test if it is on both segments
C
      if ((ltype1.eq.0).and.(ltype2.eq.0)) then
C
C     both lines not vertical
C
        if (s1.eq.s2) then
C
C     equal slopes = parallel
C
          ipar=ipar+1
        else
C
C     not parallel
C
          xi=(yint2-yint1)/(s1-s2)
          yi=s1*xi+yint1
          if (s1.eq.0) then
C
C     line 1 horizontal
C
            if ((xi.gt.xmin).and.(xi.lt.xmax).and.
     &          (xi.gt.amin).and.(xi.lt.amax).and.
     &          (yi.gt.bmin).and.(yi.lt.bmax)) then
              int=int+1
            end if
          else if (s2.eq.0) then
C
C     line 2 horizontal
C
            if ((xi.gt.xmin).and.(xi.lt.xmax).and.
     &          (xi.gt.amin).and.(xi.lt.amax).and.
     &          (yi.gt.ymin).and.(yi.lt.ymax)) then
              int=int+1
            end if
          else
C
C     neither line horizontal
C
            if ((xi.gt.xmin).and.(xi.lt.xmax).and.
     &          (xi.gt.amin).and.(xi.lt.amax).and.
     &          (yi.gt.ymin).and.(yi.lt.ymax).and.
     &          (yi.gt.bmin).and.(yi.lt.bmax)) then
              int=int+1
            end if
          end if
        end if       
      else if ((ltype1.eq.1).and.(ltype2.eq.0)) then
C
C     line 1 vertical, line 2 not
C
        xi=x1
        yi=s2*xi+yint2
        if (s2.eq.0) then
C
C     line 2 horizontal
C
          if ((xi.gt.amin).and.(xi.lt.amax).and.
     &        (yi.gt.ymin).and.(yi.lt.ymax)) then
            int=int+1
          end if
        else
C
C     line 2 not horizontal
C
          if ((xi.gt.amin).and.(xi.lt.amax).and.
     &        (yi.gt.ymin).and.(yi.lt.ymax).and.
     &        (yi.gt.bmin).and.(yi.lt.bmax)) then
            int=int+1
          end if
        end if
      else if ((ltype1.eq.0).and.(ltype2.eq.1)) then
C
C     line 2 vertical, line 1 not
C
        xi=x2
        yi=s1*xi+yint1
        if (s1.eq.0) then
C
C     line 1 horizontal
C
          if ((xi.gt.xmin).and.(xi.lt.xmax).and.
     &        (yi.gt.bmin).and.(yi.lt.bmax)) then
            int=int+1
          end if
        else
C
C     line 1 not horizontal
C
          if ((xi.gt.xmin).and.(xi.lt.xmax).and.
     &        (yi.gt.ymin).and.(yi.lt.ymax).and.
     &        (yi.gt.bmin).and.(yi.lt.bmax)) then
            int=int+1
          end if
        end if
      else
C
C     both lines vertical = parallel
C
        ipar=ipar+1
      end if
C
      Return
      End