c******PROGRAM TR
c******CONTOURS RANDOM DATA IN X,Y,Z FORMAT WITH DELAUNAY TRIANGLES
       parameter(npoint=7000,nktetr=200)
       common/pnt/ pntx(npoint),pnty(npoint),pntz(npoint),trx(npoint*2),
     &      try(npoint*2),trr(npoint*2),xpnt(3,3),det(2,3)
       integer*2 itetr,istack,ktetr,itemp
       common/stack/itetr(npoint*2,3),istack(npoint*2),
     & ktetr(nktetr,2),itemp(3,2)
       character*1 ifmt,resp,resp2,resp3
       character*24 sta
       integer jsta(6)
       character*40 ifile,ofile
       character*8 buff,passwrd,contur
       dimension isl(1),dxp(2),dyp(2),xp(4),yp(4),idim(1)
       equivalence(sta,jsta(1))
       data passwrd/'sB8yr0hF'/,contur/'contour '/
       data xmin/1.0e38/,xmax/-1.0e38/,ymin/1.0e38/,ymax/-1.0e38/
     &     ,zmin/1.0e38/,zmax/-1.0e38/,isl(1)/1/,dval/1.0e38/
       data xp(2)/0./,yp(2)/0./,idim(1)/1/
      call getarg(1,buff,length)
      if(buff.ne.passwrd) then
        print *,'program can only be run from program GEOCON'
        stop
      endif
      call getarg(2,buff,length)
1       write(*,'(a\)') ' enter input random data file name:'
        read(*,'(a)') ifile
        call opnsta(10,ifile,ifmt,ierr,*1)
        if(ierr.ne.0) then
          pause 'press ENTER to return to menu'
          stop
        endif
      if(buff.ne.contur) then
        write(*,'(a\)') ' enter output grid file name:'
        read(*,'(a)') ofile
        write(*,'(a\)') ' enter grid increment:'
        read *,del
      else
        write(*,'(2a\)') ' want to input contour interval; if no, ',
     &  'the program determines 20 levels ? [y]'
        read(*,'(a)') resp3
        if(resp3.ne.'n'.and.resp3.ne.'N') then
          write(*,'(a\)') ' enter contour interval:'
          read *,dcval
        endif
        write(*,'(a\)') ' want to plot data locations ? [y]'
        read(*,'(a)') resp
        if(resp.ne.'n'.and.resp.ne.'N') then
          write(*,'(a\)') ' want to plot data values also ? [y]'
          read(*,'(a)') resp2
        endif
      endif

      print *,'checking input data'
      i = 4
10    if(ifmt.eq.'*') then
        read(10,*,err=90,end=12,iostat=icheck) pntx(i),pnty(i),pntz(i)
      else
        read(10,err=90,end=12,iostat=icheck) pntx(i),pnty(i),pntz(i)
      endif
      if(pntx(i).lt.xmin) xmin = pntx(i)
      if(pntx(i).gt.xmax) xmax = pntx(i)
      if(pnty(i).lt.ymin) ymin = pnty(i)
      if(pnty(i).gt.ymax) ymax = pnty(i)
      if(pntz(i).lt.zmin) zmin = pntz(i)
      if(pntz(i).gt.zmax) zmax = pntz(i)
      i = i + 1
      if(i.gt.npoint) then
        print *,'number of input points > ',npoint
        pause 'press ENTER to return to menu'
        stop
      endif
      go to 10
12    icnt = i -1
c*****ELIMINATE DUPLICATE POINTS
      do 22 i = 4,icnt-1
      if(pntx(i).eq.dval) go to 22
      do 21 j = i+1,icnt
      if(pntx(i).eq.pntx(j).and.pnty(i).eq.pnty(j)) pntx(j) = dval
21    continue
22    continue
      m = 5
      do 23 i=4,icnt -1
        if(pntx(i+1).eq.dval) go to 23
        pntx(m) = pntx(i+1) 
        pnty(m) = pnty(i+1) 
        pntz(m) = pntz(i+1) 
        m = m + 1
23    continue
      m = m - 1
      icnt = m
      if(buff.ne.contur) then
        if(xmax.le.xmin.or.ymax.le.ymin) then
          print *,' something wrong with input data'
          go to 99
        endif
        nc = nint((xmax-xmin)/abs(del)) + 1
        nr = nint((ymax-ymin)/abs(del)) + 1
        open(12,file='tr.dat',status='unknown')
        write(12,'(2a)') ' output file= ',ofile
        write(12,*) xmin,xmax,ymin,ymax
        write(12,*) del,nc,nr
        do 25 i=4, icnt
        write(12,*) pntx(i),pnty(i),pntz(i)
25      continue
      endif
      if(buff.eq.contur) then
c*****SET UP PLOTTER
        call pltset(9,xbrd,ybrd,isl)
        dxp(1) = xmin
        dxp(2) = xmax
        dyp(1) = ymin
        dyp(2) = ymax
        xp(1) = xbrd - 1.
        xp(3) = .1
        xp(4) = xp(1) + xp(3) + .9
        yp(1) = ybrd - 1.
        yp(3) = .1
        yp(4) = yp(1) + yp(3) + .9
        call scale(dxp,dyp,xp,yp,4,ier)
        if(ier.ne.0) then
          print *,'error in subroutine scale'
          pause 'press ENTER to return to menu'
          stop
        endif
        call neatl
        if(resp.ne.'n'.and.resp.ne.'N') then
c*****PLOT DATA
          do 14 i = 4,icnt
          call vchar(pntx(i),pnty(i),4,1,2,.06,0.,0.,0.)
          if(resp2.ne.'n'.and.resp2.ne.'N') then
            write(sta,'(f9.3)') pntz(i)
            n = leftj(sta)
            call vchar(pntx(i),pnty(i),jsta,n,2,.10,.3,.1,0.)
          endif
14      continue
        endif
      endif
c*****NORMALIZE THE DATA
      datax = xmax - xmin
      datay = ymax - ymin
      do 15 i = 4,icnt
      pntx(i) = (pntx(i)-xmin)/datax
      pnty(i) = (pnty(i)-ymin)/datay
15    continue
	ITEMP(1,1)=1
	ITEMP(2,1)=1
	ITEMP(3,1)=2
	ITEMP(1,2)=2
	ITEMP(2,2)=3
	ITEMP(3,2)=3

	XPNT(1,1)=-1
	XPNT(1,2)=-1
	XPNT(1,3)=2
	XPNT(2,1)=5
	XPNT(2,2)=-1
	XPNT(2,3)=2
	XPNT(3,1)=-1
	XPNT(3,2)=5
	XPNT(3,3)=18
C******FORM TRIANGLES
       if(buff.ne.contur) then
         print *,'forming triangles'
       else
         write(*,'(a)')'  forming triangles'
       endif
       call formtr(isp,icnt)
       if(buff.ne.contur) then
         print *,'augmenting data'
       else
         write(*,'(a)')'+                  '
       endif
       if(resp3.eq.'n'.or.resp3.eq.'N'.or.buff.ne.contur) then
c******FIND REASONABLE CONTOUR INTERVAL
         call set(zmin,zmax,dcval,20)
       endif
       temin=zmin/dcval
       icont=int(temin)
       if(aint(temin).ne.temin.and.temin.gt.0) icont=icont+1
c******CONTOUR LOOP
16     cont=icont*dcval
       if(cont.gt.zmax) go to 17
       call trcont(isp,cont,xmin,datax,ymin,datay,buff)
       icont = icont + 1
       go to 16
17     continue
       call endpt(idim)
      stop
90    print *,'error no. ',icheck,' reading input file'
99    pause 'press ENTER to return to menu'
      stop
      end
c***********************************************************************
      subroutine opnsta(ista,ifile2,ifmtv,ierr,*)
      character*(*) ifile2,ifmtv
      character*1 resp
      logical exists,opn
      ierr=0
      inquire(file=ifile2,exist=exists)
        if(.not.exists) then
          print *,'can''t find file ',ifile2
          ifile2=' '
          write(*,'(a\)') ' TRY AGAIN ? [y]: '
          read(*,'(a)') resp
          if(resp.eq.'n'.or.resp.eq.'N') then
            ierr=1
            return 
          else
            return 1
          endif
        endif
        inquire(unit=ista,opened=opn)
          if(opn) close(ista)
        open(unit=ista,file=ifile2,
     &  status='old',err=442,iostat=icheck)
        read(ista,*,err=442,end=442,iostat=icheck) x,y,z
        rewind ista
        ifmtv(1:1) = '*'
        return
442     close(ista)
        open(unit=ista,file=ifile2,form='unformatted',
     &  status='old',err=443,iostat=icheck)
        read(ista,err=444,end=444,iostat=icheck) x,y,z
        nword = lenrec(ista,'u')
        if(nword.ne.3.and.nword.ne.-1) go to 444
        ifmtv = ' '
        return
443    print *,'error no. ',icheck,' with file: ',ifile2
       go to 445
444    print *,'file is not an xyz file'
445    ifile2=' '
       close(ista)
       ierr = 1
       return 
       end
c***********************************************************************
	function lenrec( ista,rtype )
c  brute force recovery of record length
       parameter(nwork=50)
	real irec(nwork)
       character*1 rtype

	istep = 10
	rewind ista

	do 5 max = istep, nwork, istep
         if(rtype.eq.'u') then
           read( ista, err=10 ) ( irec(i), i=1, max )
         else
	    read( ista,*, err=10 ) ( irec(i), i=1, max )
         endif
	  rewind ista
5       continue
	lenrec = -1
	go to 99

10	rewind ista
	do 30 j = max, max - istep, -1
         if(rtype.eq.'u') then
	    read( ista, err=20 ) ( irec(i), i=1, j )
         else
	    read( ista,*, err=20 ) ( irec(i), i=1, j )
         endif
	  lenrec=j
	  go to 99
20	  rewind ista
30      continue

99	rewind ista
	return
	end
c***********************************************************************
       subroutine formtr(isp,icnt)
       parameter(npoint=7000,nktetr=200)
       common/pnt/ pntx(npoint),pnty(npoint),pntz(npoint),trx(npoint*2),
     &      try(npoint*2),trr(npoint*2),xpnt(3,3),det(2,3)
       integer*2 itetr,istack,ktetr,itemp
       common/stack/itetr(npoint*2,3),istack(npoint*2),
     & ktetr(nktetr,2),itemp(3,2)
	ISP=1
	ID=2

       trx(1) = xpnt(1,3)
       try(1) = xpnt(2,3)
       trr(1) = xpnt(3,3)
	do 2 I=1 , 3
		ITETR(1,I)=I
c		TRR(I)=XPNT(I,3)
		PNTX(I)=XPNT(I,1)
		PNTY(I)=XPNT(I,2)
2      continue
	do 6 I=2 , npoint*2
              ISTACK(I)=I
6      continue
	do 50 NUC=4 , ICNT
		KM=0
		do 30 JT=1 , ISP	
c******TEST IF NEW DATA POINT IS WITHIN THE JT CIRCUMCIRCLE
 			DX=TRX(JT)-PNTX(NUC)
 			DY=TRY(JT)-PNTY(NUC)
c			DSQ=TRR(JT)-(DX*DX + DY*DY)
                     dsq=trr(jt) -(dx*dx)
 			IF (DSQ.lt.0.0) GO TO 30
 			dsq = dsq -(dy*dy)
 			IF (DSQ.lt.0.0) GO TO 30
c******DELETE THIS 3-TUPLE BUT SAVE ITS EDGES
			ID=ID-1
			ISTACK(ID)=JT
c******ADD EDGES TO KTETR BUT DELETE IF ALREADY LISTED
			do 28 I=1 , 3		
				L1=ITEMP(I,1)
				L2=ITEMP(I,2)
				IF (KM.le.0)  GO TO 26
				KMT=KM
				do 24 J=1 , KMT		
				   IF (ITETR(JT,L1).ne.KTETR(J,1))GO TO 24
        			   IF (ITETR(JT,L2).ne.KTETR(J,2))GO TO 24
					KM=KM-1
					IF (J.gt.KM) GO TO 28
					do 20 K=J , KM		
						K1=K+1
						KTETR(K,1)=KTETR(K1,1)
						KTETR(K,2)=KTETR(K1,2)
20					continue
					GO TO 28
24				continue
26				KM=KM+1
			if(km.gt.nktetr) then
 			  print *,'data not suitable for this program'
                       pause 'press ENTER to return to menu'
                       stop
                     endif
				KTETR(KM,1)=ITETR(JT,L1)
				KTETR(KM,2)=ITETR(JT,L2)
28			continue
30			continue
c******FORM NEW 3-TUPLES
		do 48 I=1 , KM					
			KT=ISTACK(ID)
			ID=ID+1
c******CALCULATE THE CIRCUMCIRCLE CENTER AND RADIUS
c******SQUARED OF POINTS KTETR(I,*) AND PLACE IN TRX(KT),TRY(KT),TRR(KT)
			do 44 JZ=1 , 2				
				I2=KTETR(I,JZ)
				DET(JZ,1)=PNTX(I2)-PNTX(NUC)
				DET(JZ,2)=PNTY(I2)-PNTY(NUC)
				DET(JZ,3)=DET(JZ,1)*(PNTX(I2)+PNTX(NUC))/2.0
     &			       + DET(JZ,2)*(PNTY(I2)+PNTY(NUC))/2.0
44			continue
			DD=DET(1,1)*DET(2,2)-DET(1,2)*DET(2,1)
			TRX(KT)=(DET(1,3)*DET(2,2) - DET(2,3)*DET(1,2))/DD
			TRY(KT)=(DET(1,1)*DET(2,3) - DET(2,1)*DET(1,3))/DD
			TRR(KT)=(PNTX(NUC)-TRX(KT))**2 +
     &		       (PNTY(NUC)-TRY(KT))**2
			ITETR(KT,1)=KTETR(I,1)
			ITETR(KT,2)=KTETR(I,2)
			ITETR(KT,3)=NUC
48		continue
		if(isp.ge.(npoint*2-1)) then
		   print *,'number of triangles > ',npoint*2
                 pause 'press ENTER to return to menu'
		   stop
		endif
		ISP=ISP+2
50	continue
       return
       end
c***********************************************************************
      subroutine trcont(isp,cur,xmin,datax,ymin,datay,buff)
       parameter(npoint=7000,nktetr=200)
       common/pnt/ pntx(npoint),pnty(npoint),pntz(npoint),trx(npoint*2),
     &      try(npoint*2),trr(npoint*2),xpnt(3,3),det(2,3)
       integer*2 itetr,istack,ktetr,itemp
       common/stack/itetr(npoint*2,3),istack(npoint*2),
     & ktetr(nktetr,2),itemp(3,2)
       dimension x(3),y(3)
       character * 8 buff,contur
       character*24 sta
       integer jsta(6)
       equivalence(sta,jsta(1))
       data contur/'contour '/

		do 90 JT=1 , ISP			
			IF (ITETR(JT,1).lt.4.OR.TRR(JT).gt.1) GO TO 90
c******FIND CONTOUR CUTS
			ABIT=0.0
			IF (PNTZ(ITETR(JT,1)).eq.PNTZ(ITETR(JT,2)).OR.
     &			PNTZ(ITETR(JT,1)).eq.PNTZ(ITETR(JT,3)).OR.
     &			PNTZ(ITETR(JT,2)).eq.PNTZ(ITETR(JT,3)))
     &		       ABIT=1.0E-10
			TOP=amax1(PNTZ(ITETR(JT,1)),PNTZ(ITETR(JT,2)),
     &			PNTZ(ITETR(JT,3)))
			BOT=amin1(PNTZ(ITETR(JT,1)),PNTZ(ITETR(JT,2)),
     &			PNTZ(ITETR(JT,3)))
			IF (CUR.gt.TOP.OR.CUR.lt.BOT) GO TO 90
			CZ=(CUR-PNTZ(ITETR(JT,1)))/(PNTZ(ITETR(JT,2))
     &			-PNTZ(ITETR(JT,1)) +ABIT)
			IF (CZ.le.0.0.OR.CZ.ge.1.0) GO TO 83
			X1=(PNTX(ITETR(JT,1)) +(PNTX(ITETR(JT,2))
     &			-PNTX(ITETR(JT,1))) *CZ)
			Y1=(PNTY(ITETR(JT,1))+(PNTY(ITETR(JT,2))
     &			-PNTY(ITETR(JT,1))) *CZ)
82			CZ=(CUR-PNTZ(ITETR(JT,1)))/(PNTZ(ITETR(JT,3))
     &			-PNTZ(ITETR(JT,1)) +ABIT)
			IF (CZ.lt.0.0.OR.CZ.gt.1.0) GO TO 84
			X2=(PNTX(ITETR(JT,1)) +(PNTX(ITETR(JT,3))
     &			-PNTX(ITETR(JT,1))) *CZ)
			Y2=(PNTY(ITETR(JT,1))+(PNTY(ITETR(JT,3))
     &			-PNTY(ITETR(JT,1))) *CZ)
			GOTO 85
83			CZ=(CUR-PNTZ(ITETR(JT,1)))/(PNTZ(ITETR(JT,3))
     &			-PNTZ(ITETR(JT,1)) +ABIT)
			IF (CZ.lt.0.0.OR.CZ.gt.1.0)  GO TO 90
			X1=(PNTX(ITETR(JT,1)) +(PNTX(ITETR(JT,3))
     &			-PNTX(ITETR(JT,1))) *CZ)
			Y1=(PNTY(ITETR(JT,1))+(PNTY(ITETR(JT,3))
     &			-PNTY(ITETR(JT,1))) *CZ)
84			CZ=(CUR-PNTZ(ITETR(JT,2)))/(PNTZ(ITETR(JT,3))
     &			-PNTZ(ITETR(JT,2)) +ABIT)
			IF (CZ.lt.0.0.OR.CZ.gt.1.0) GO TO 90
			X2=(PNTX(ITETR(JT,2)) +(PNTX(ITETR(JT,3))
     &			-PNTX(ITETR(JT,2))) *CZ)
			Y2=(PNTY(ITETR(JT,2))+(PNTY(ITETR(JT,3))
     &			-PNTY(ITETR(JT,2)))*CZ)
c******UNNORMALIZE FOR PLOTTING
85            x(1) = x1 * datax + xmin
              x(2) = x2 * datax + xmin
              y(1) = y1 * datay + ymin
              y(2) = y2 * datay + ymin
              xdif = (x(2)-x(1)) * .25
              ydif = (y(2)-y(1)) * .25
              xx = x(1) + xdif
              yy = y(1) + ydif
              xxx = x(2) - xdif
              yyy = y(2) - ydif
              if(buff.eq.contur) then
                call line(x,y,2,0,0)
              else
                write(12,*) xx,yy,cur
                write(12,*) xxx,yyy,cur
              endif
c              x(1) = pntx(itetr(jt,1)) * datax + xmin
c              x(2) = pntx(itetr(jt,2)) * datax + xmin
c              x(3) = pntx(itetr(jt,3)) * datax + xmin
c              y(1) = pnty(itetr(jt,1)) * datay + ymin
c              y(2) = pnty(itetr(jt,2)) * datay + ymin
c              y(3) = pnty(itetr(jt,3)) * datay + ymin
c              call line(x,y,3,0,1400)
90		continue
87	return
       end
c***********************************************************************
        function leftj(a)
c  left justifies a string and returns the position
c  of the last nonblank character
        character a*(*)
        n=len(a)
        if(a(1:1).ne.' ') go to 15
        do 1 m=2,n
1       if(a(m:m).ne.' ') go to 5
        leftj=0
        return
5       i2=1
        do 10 i=m,n
        a(i2:i2)=a(i:i)
        a(i:i)=' '
10      i2=i2+1
        n=n-m+1
15      do 20 leftj=n,1,-1
20      if(a(leftj:leftj).ne.' ') go to 25
25      return
        end
c***********************************************************************
      subroutine set(zmin,zmax,dx,maxint)
c  adjust interval 
c  input x-min max range and maxint-intervals
c  returns dx
      ixpn(r)=int(alog10(abs(r))+100.)-100
      t=(zmax-zmin)/float(maxint)
      if(t.lt.1.e-20) return
      p10=10.**ixpn(t)
      p10=sign(1.0,t)*p10
      t1=t/p10
      if(t1.le.1.0) dx=p10
      if(t1.gt.1.0) dx=2.*p10
      if(t1.gt.2.0) dx=5.*p10
      if(t1.gt.5.0) dx=10.*p10
      return
      end
