c  Program DETOURG
c
c  Contouring of gridded data using Delaunay tessellation
c
c  modified from  D.F. Watson by J. Phillips
c
c  reference: Watson, D.F., 1982, Acord: Automatic contouring of raw data:
c             Computers and Geosciences, v.8, no. 1, p.97-101.
c
      dimension xyz(1203,3), tetr(2401,3), xpnt(3,3), det(2,3)
      dimension isl(5), dxp(2), dyp(2), xp(4), yp(4), x(2), y(2),
     +          ie(10), ipn(12)
c      dimension title(14),pgm(2),row(1000,2)
      dimension row(1000,2)
      integer itetr(2401,3), istack(2401), ktetr(50,2), itemp(3,2)
      character ifile*80, tri*1, dat*1, label*9, title*56, pgm*8
      data itemp,xpnt/1,1,2,2,3,3,-1.,5.,-1.,-1.,-1.,5.,2.,2.,18./
c      data ipn/100,300,200,600,400,500/
      data ipn/100,900,300,1100,200,1000,1400,600,400,1200,500,1300/
      data ddval/1.0e30/
      print *,'input grid file: '
      read(5,100) ifile
  100 format(a)
      open(10,file=ifile,form='unformatted',status='old')
      read(10) title,pgm,nc,nr,nz,xo,dx,yo,dy
      xmin=xo
      ymin=yo
      xmax=xo+(nc-1)*dx
      ymax=yo+(nr-1)*dy
      zmin=1.e+38
      zmax=-1.e+38
      do j=1,nr
      read(10) dlt,(row(i,1),i=1,nc)
      do i=1,nc
      if(row(i,1).gt.ddval) go to 23
      if(row(i,1).gt.zmax) zmax=row(i,1)
      if(row(i,1).lt.zmin) zmin=row(i,1)
   23 continue
      end do
      end do
      rewind 10
      read(10) title
      print *,'enter iplotr (5=HP,8=CGA,9=EGA,10=VGA):'
      read(5,*) iplotr
      size=1.0
      print *,'xmin = ',xmin,' xmax = ',xmax
      print *,'ymin = ',ymin,' ymax = ',ymax
      print *,'zmin = ',zmin,' zmax = ',zmax
      print *, 'enter contour interval:'
      read(5,*) conti
c      contm=zmin-amod(zmin,conti)+conti
      contm=zmin-amod(zmin,conti)
      print *,'minimum contour level = ',contm
c      lc=int((zmax-contm)/conti)
      lc=int((zmax-contm)/conti)+1
      print *,'number of contours = ',lc
      print *
      print *,'Press ENTER to continue'
      lc=lc+1
      if(contm.eq.zmin) contm=contm+.01*conti
      read(5,100) tri
      xr=xmax-xmin
      datax=xr
      yr=ymax-ymin
      if(yr.gt.datax) datax=yr
      call pltset(iplotr,xboard,yboard,isl)
      dxp(1)=0.0
      dxp(2)=xr/datax
      dyp(1)=0.0
      dyp(2)=yr/datax
      yp(1)=yboard-0.7
      xp(1)=yp(1)*xr/yr
      if(xp(1).gt.xboard-1.2) then
        xp(1)=xboard-1.2
        yp(1)=xp(1)*yr/xr
      endif
      xp(2)=0.
      xp(3)=0.1
      xp(4)=xboard
      yp(2)=0.
      yp(3)=0.6
      yp(4)=yboard
      nopts=4
      call scale(dxp,dyp,xp,yp,nopts,ier)
      if(ier.ne.0) stop 'could not scale the plot'
      call neatl
      call vchar(.15,.35,ifile,80,3,.1,0.,0.,0.)
      lt=len_trim(title)
      x=xp(3)+xp(1)-lt*0.1+.05
      if(lt.gt.0) call vchar(x,.35,title,lt,3,.1,0.,0.,0.)
      read(10) dlt,(row(i,1),i=1,nc)
      do 200 j=2,nr
        do 1 i=1,3
          itetr(1,i)=i
          tetr(1,i)=xpnt(i,3)
          do 1 jj=1,3
    1   xyz(i,jj)=xpnt(i,jj)
        do 2 i=2,2401
    2   istack(i)=i
        read(10) dlt,(row(i,2),i=1,nc)
        ix=3
        ilast=0
        do 290 i=1,nc-1
          boxmin=min(row(i,1),row(i+1,1),row(i,2),row(i+1,2))
          boxmax=max(row(i,1),row(i+1,1),row(i,2),row(i+1,2))
          cont=contm
          do ic=1,lc
            if(boxmin.le.cont.and.boxmax.ge.cont) go to 299
            cont=cont+conti
          end do
          go to 290
299       continue
c          print *,boxmin,cont,boxmax
          if(ilast.ne.i) then
            ix=ix+1
            xyz(ix,1)=xo+float(i-1)*dx
            xyz(ix,2)=yo+float(j-2)*dy
            xyz(ix,3)=row(i,1)
            ix=ix+1
            xyz(ix,1)=xo+float(i-1)*dx
            xyz(ix,2)=yo+float(j-1)*dy
            xyz(ix,3)=row(i,2)
          endif
          ix=ix+1
          xyz(ix,1)=xo+float(i)*dx
          xyz(ix,2)=yo+float(j-2)*dy
          xyz(ix,3)=row(i+1,1)
          ix=ix+1
          xyz(ix,1)=xo+float(i)*dx
          xyz(ix,2)=yo+float(j-1)*dy
          xyz(ix,3)=row(i+1,2)
c            if(xyz(i,3).gt.ddval) xyz(i,3)=0.001
          ilast=i+1
290     continue
c        print *,ix
c  normalize the data
        do ic=4,ix
          xyz(ic,1)=(xyz(ic,1)-xmin)/datax
          xyz(ic,2)=(xyz(ic,2)-ymin)/datax
        end do

        isp=1
        id=2
        do 13 nuc=4,ix
          km=0
c  loop thru the established 3-tuples
          do 10 jt=1,isp
c  test if new data point is within the jt circumcircle
            dsq=tetr(jt,3)-(xyz(nuc,1)-tetr(jt,1))**2
            if(dsq.lt.0.0) go to 10
            dsq=dsq-(xyz(nuc,2)-tetr(jt,2))**2
            if(dsq.lt.0.0) go to 10
c  delete this 2-tuple but save its edges
            id=id-1
            istack(id)=jt
c  add edges to ktetr but delete if already listed
            do 9 ii=1,3
              l1=itemp(ii,1)
              l2=itemp(ii,2)
              if(km.le.0) go to 8
              kmt=km
              do 7 jj=1,kmt
                if(itetr(jt,l1).ne.ktetr(jj,1)) go to 7
                if(itetr(jt,l2).ne.ktetr(jj,2)) go to 7
                km=km-1
                if(jj.gt.km) go to 9
                do 6 k=jj,km
                  k1=k+1
                  do 6 l=1,2
    6             ktetr(k,l)=ktetr(k1,l)
                go to 9
    7         continue
    8         km=km+1
              ktetr(km,1)=itetr(jt,l1)
              ktetr(km,2)=itetr(jt,l2)
    9       continue
   10     continue
c  form new 3-tuples
          do 12 ii=1,km
            kt = istack(id)
            id = id + 1
c  calculate the circumcircle center and radius
c  squared of points ktetr(i,*) and place in tetr(kt,*)
            do 11 jz = 1,2
              i2 = ktetr(ii,jz)
              det(jz,1) = xyz(i2,1) - xyz(nuc,1)
              det(jz,2) = xyz(i2,2) - xyz(nuc,2)
   11         det(jz,3) = det(jz,1) * (xyz(i2,1) + xyz(nuc,1))/2.0
     +          + det(jz,2) * (xyz(i2,2) + xyz(nuc,2))/2.0
            dd = det(1,1) * det(2,2) - det(1,2) * det(2,1)
            if(dd.eq.0.0) stop 'dd'
            tetr(kt,1) = (det(1,3) * det(2,2) - det(2,3) * det(1,2))/dd
            tetr(kt,2) = (det(1,1) * det(2,3) - det(2,1) * det(1,3))/dd
            tetr(kt,3) = (xyz(nuc,1) - tetr(kt,1))**2
     +           + (xyz(nuc,2) - tetr(kt,2))**2
            itetr(kt,1) = ktetr(ii,1)
            itetr(kt,2) = ktetr(ii,2)
   12       itetr(kt,3) = nuc
   13     isp = isp + 2
        do 19 jt = 1,isp
          if(itetr(jt,1) .lt. 4 .or. tetr(jt,3) .gt. 1.) go to 19
c  find contour intersections
          abit = 0.0
          if(xyz(itetr(jt,1),3) .eq. xyz(itetr(jt,2),3) .or.
     +   xyz(itetr(jt,1),3) .eq. xyz(itetr(jt,3),3) .or.
     +   xyz(itetr(jt,2),3) .eq. xyz(itetr(jt,3),3)) abit = 1.e-10
          top = amax1(xyz(itetr(jt,1),3),xyz(itetr(jt,2),3),
     +            xyz(itetr(jt,3),3))
          bot = amin1(xyz(itetr(jt,1),3),xyz(itetr(jt,2),3),
     +            xyz(itetr(jt,3),3))
          do 18 jc =  1,lc
            if(abit.ne.0.0) go to 18
            cont=contm+(jc-1)*conti
            ip=1+jc*12/(lc+1)
            if(cont .gt. top .or.  cont .lt. bot) go to 18
            cz = (cont - xyz(itetr(jt,1),3))/(xyz(itetr(jt,2),3)
     +           - xyz(itetr(jt,1),3) + abit)
            if(cz .le. 0.0 .or. cz .ge. 1.0) go to 15
            x1 = (xyz(itetr(jt,1),1) + (xyz(itetr(jt,2),1)
     +   -  xyz(itetr(jt,1),1)) * cz) * size
            y1 = (xyz(itetr(jt,1),2) + (xyz(itetr(jt,2),2)
     +   -  xyz(itetr(jt,1),2)) * cz) * size
   14       cz = (cont - xyz(itetr(jt,1),3))/(xyz(itetr(jt,3),3)
     +           - xyz(itetr(jt,1),3) + abit)
            if(cz .lt. 0.0 .or. cz .gt. 1.0) go to 16
            x2 = (xyz(itetr(jt,1),1) + (xyz(itetr(jt,3),1)
     +   -  xyz(itetr(jt,1),1)) * cz) * size
            y2 = (xyz(itetr(jt,1),2) + (xyz(itetr(jt,3),2)
     +   -  xyz(itetr(jt,1),2)) * cz) * size
            go to 17
   15       cz = (cont - xyz(itetr(jt,1),3))/(xyz(itetr(jt,3),3)
     +           - xyz(itetr(jt,1),3) + abit)
            if(cz .lt. 0.0 .or. cz .gt. 1.0) go to 18
            x1 = (xyz(itetr(jt,1),1) + (xyz(itetr(jt,3),1)
     +   -  xyz(itetr(jt,1),1)) * cz) * size
            y1 = (xyz(itetr(jt,1),2) + (xyz(itetr(jt,3),2)
     +   -  xyz(itetr(jt,1),2)) * cz) * size
   16       cz = (cont - xyz(itetr(jt,2),3))/(xyz(itetr(jt,3),3)
     +           - xyz(itetr(jt,2),3) + abit)
            if(cz .lt. 0.0 .or. cz .gt. 1.0) go to 18
            x2 = (xyz(itetr(jt,2),1) + (xyz(itetr(jt,3),1)
     +   -  xyz(itetr(jt,2),1)) * cz) * size
            y2 = (xyz(itetr(jt,2),2) + (xyz(itetr(jt,3),2)
     +   -  xyz(itetr(jt,2),2)) * cz) * size
   17       continue
            x(1)=x1
            x(2)=x2
            y(1)=y1
            y(2)=y2
            call line(x,y,2,0,ipn(ip))
   18     continue
   19   continue
c        do 20 ii=4,ix
cc          print *,xyz(ii,1),xyz(ii,2)
c          call vchar(xyz(ii,1),xyz(ii,2),13,1,1500,.1,0.,0.,0.)
c   20   continue
        do 200 ic=1,nc
        row(ic,1)=row(ic,2)
  200 continue
c plot scale
      ipold=0
      y=.3
      do 22 i=1,lc
      ip=1+i*12/(lc+1)
      cont=contm+(i-1)*conti
      if(ip.ne.ipold) then
c      lab(1:1)="'"
c      lab(8:8)="'"
      write(label(1:9),'(g9.3)') cont
      y=y+.4
c      print *,lab
c      call vchar(6.0,y,3,1,ipn(ip)+3,.1,0.,0.,0.)
      call vchar(xp(1)+0.4,y,label,9,ipn(ip)+3,.1,0.,0.,0.)
c      call vchar(6.0,y,'   333',6,ipn(ip)+3,.1,0.,0.,0.)
      endif
      ipold=ip
   22 continue
   21 call endpt(ie)
      stop
      end

