c  ASDEP
c
c  Crude magnetic depth estimates from the analytic signal amplitude (ASA)
c  grid.  Applying this method to the ASA of magnetic or pseudomagnetic data
c  should yield reasonably accurate depths to the tops of very thick bodies.
c  Applying this method to the ASA of gravity or pseudogravity data should
c  yield reasonably accurate depths to very thin bodies.
c
c  A 5x5 data window is passed over the ASA grid (normally created by AS or
C  AS2).  Within each window, an attempt is made to find several maxima
c  and determine a linear trend for the maxima.  If the attempt is successful,
c  a depth and standard error are estimated from data within a belt
c  perpendicular to the linear trend.  Normally a user should start by plotting
c  data windows in order to get a feel for the largest acceptable standard
c  error.  Once the user is confident of his parameters, he/she can restart the
c  program by entering the letter 'q'.

c  Output is sent to the post file ASDEP.PST, which contains the fields:
c  post(1)=depth, post(2)=%error_in_depth, and post(3)=strike_of_contact.
c  This data can be plotted with the programs PLOTPST or DETOUR.
c
c  Written by Jeff Phillips.
c
      dimension as(2000,5),at(2000)
      dimension id(14),pgm(2),pk(3,100)
      character ifile*80, pans*1, nans*1, dum*1, sta*8
      data zero/0.0/,sta/'        '/
      data ddval/1.e+20/
      rad=3.1415927/180.
      pi2=3.1415927/2.
  100 format(a)
      print *,'Enter analytic signal magnitude file:'
      read(5,100) ifile
   10 open(14,file=ifile,form='unformatted',status='old')
      open(15,file='asdep.pst',form='unformatted',status='unknown')
c      print *,'Enter the number of maxima used to get trend'
c      read(5,*) maxn
      print *,'Enter the skip factor (normally 0):'
      read(5,*) iskip
      is=iskip+1
      print *,'Enter largest acceptable percent standard error'
      print *,'for the depth estimates:'
      read(5,*) dlim
      print *,'Enter height above datum in same units as dx,dy'
      print *,'(all depths will be relative to the datum):'
      read(5,*) dmin
      print *,'Keep negative depths?:'
      read(5,100) nans
      if(nans.eq.'N') nans='n'
      print *,'plot data windows?:'
      read(5,100) pans
      read(14) id,pgm,nc,nr,nz,xo,dx,yo,dy
      dxs=dx*float(is)
      dys=dy*float(is)
      do 60 kk=1,is
      do 20 j=1,4
c   20 read(14) dlt,(as(i,j),i=1,nc)
      read(14) dlt,(at(i),i=1,nc)
      do 19 i=1,nc
      if(at(i).ge.ddval) then
        as(i,j)=at(i)
      else
        as(i,j)=at(i)**2
      endif
   19 continue
      if(iskip.gt.0) then
        do i=1,iskip
        read(14) dlt
        enddo
      endif
   20 continue
c      del=sqrt(dx*dx+dy*dy)
      del=sqrt(dxs*dxs+dys*dys)
      numdep=0
c  Begin main loop through rows.
c      do 30 k=3,nr-2
c      read(14) dlt,(as(i,5),i=1,nc)
      do 30 k=kk+2*is,nr-2*is-kk+1,is
      read(14) dlt,(at(i),i=1,nc)
      do 35 i=1,nc
      if(at(i).ge.ddval) then
        as(i,5)=at(i)
      else
        as(i,5)=at(i)**2
      endif
   35 continue
c  Begin main loop through columns.
c      do 40 i=3,nc-2
      do 40 i=3+2*iskip,nc-2-2*iskip
      do jjj=1,5
c      do iii=i-2,i+2
      do iii=i-2*is,i+2*is,is
      if(as(iii,jjj).ge.ddval) go to 40
      enddo
      enddo
c  Try to find maxima along 12 3-point transects through the window.
      npeak=0
c      call interp(as(i-1,3),as(i,3),as(i+1,3),dx,xmax,ymax,ierr)
      call interp(as(i-is,3),as(i,3),as(i+is,3),dxs,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=xmax
        pk(2,npeak)=0.0
        pk(3,npeak)=ymax
      endif
      call interp(as(i,2),as(i,3),as(i,4),dys,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=0.0
        pk(2,npeak)=xmax
        pk(3,npeak)=ymax
      endif
c      call interp(as(i-1,2),as(i,3),as(i+1,4),del,xmax,ymax,ierr)
      call interp(as(i-is,2),as(i,3),as(i+is,4),del,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=dxs*xmax/del
        pk(2,npeak)=dys*xmax/del
        pk(3,npeak)=ymax
      endif
c      call interp(as(i+1,2),as(i,3),as(i-1,4),del,xmax,ymax,ierr)
      call interp(as(i+is,2),as(i,3),as(i-is,4),del,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=-dxs*xmax/del
        pk(2,npeak)=dys*xmax/del
        pk(3,npeak)=ymax
      endif
c  If less than 2 maxima in the central area - move the window.
      if(npeak.lt.2) go to 40
c      call interp(as(i-1,2),as(i,2),as(i+1,2),dx,xmax,ymax,ierr)
      call interp(as(i-is,2),as(i,2),as(i+is,2),dxs,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=xmax
        pk(2,npeak)=-dys
        pk(3,npeak)=ymax
      endif
c      call interp(as(i-1,4),as(i,4),as(i+1,4),dx,xmax,ymax,ierr)
      call interp(as(i-is,4),as(i,4),as(i+is,4),dxs,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=xmax
        pk(2,npeak)=dys
        pk(3,npeak)=ymax
      endif
c      call interp(as(i-1,2),as(i-1,3),as(i-1,4),dy,xmax,ymax,ierr)
      call interp(as(i-is,2),as(i-is,3),as(i-is,4),dys,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=-dxs
        pk(2,npeak)=xmax
        pk(3,npeak)=ymax
      endif
c      call interp(as(i+1,2),as(i+1,3),as(i+1,4),dy,xmax,ymax,ierr)
      call interp(as(i+is,2),as(i+is,3),as(i+is,4),dys,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=dxs
        pk(2,npeak)=xmax
        pk(3,npeak)=ymax
      endif
c      call interp(as(i,1),as(i-1,2),as(i-2,3),del,xmax,ymax,ierr)
      call interp(as(i,1),as(i-is,2),as(i-2*is,3),del,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=dxs*(-xmax/del-1)
        pk(2,npeak)=dys*(xmax/del-1)
        pk(3,npeak)=ymax
      endif
c      call interp(as(i,1),as(i+1,2),as(i+2,3),del,xmax,ymax,ierr)
      call interp(as(i,1),as(i+is,2),as(i+2*is,3),del,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=dxs*(xmax/del+1)
        pk(2,npeak)=dys*(xmax/del-1)
        pk(3,npeak)=ymax
      endif
c      call interp(as(i-2,3),as(i-1,4),as(i,5),del,xmax,ymax,ierr)
      call interp(as(i-2*is,3),as(i-is,4),as(i,5),del,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=dxs*(xmax/del-1)
        pk(2,npeak)=dys*(xmax/del+1)
        pk(3,npeak)=ymax
      endif
c      call interp(as(i+2,3),as(i+1,4),as(i,5),del,xmax,ymax,ierr)
      call interp(as(i+2*is,3),as(i+is,4),as(i,5),del,xmax,ymax,ierr)
      if(ierr.eq.0) then
        npeak=npeak+1
        pk(1,npeak)=dxs*(-xmax/del+1)
        pk(2,npeak)=dys*(xmax/del+1)
        pk(3,npeak)=ymax
      endif
c  Sort the maxima by magnitude in descending order.
         call bsort(pk,npeak,ierr)
c  Use weighted least squares to fit a line to the maxima.
         kkk=npeak
         xmin=pk(1,1)
         xmax=pk(1,1)
         ymin=pk(2,1)
         ymax=pk(2,1)
         do iii=2,kkk
         if(pk(1,iii).gt.xmax) xmax=pk(1,iii)
         if(pk(1,iii).lt.xmin) xmin=pk(1,iii)
         if(pk(2,iii).gt.ymax) ymax=pk(2,iii)
         if(pk(2,iii).lt.ymin) ymin=pk(2,iii)
         end do
         xrange=xmax-xmin
         yrange=ymax-ymin
         x=0.
         x2=0.
         y=0.
         y2=0.
         xy=0.
         sw=0.
         do iii=1,kkk
         w=(kkk-iii+1)**2
         x=x+w*pk(1,iii)
         x2=x2+w*pk(1,iii)**2
         y2=y2+w*pk(2,iii)**2
         y=y+w*pk(2,iii)
         xy=xy+w*pk(1,iii)*pk(2,iii)
         sw=sw+w
         end do
         if(xrange.ge.yrange) then
           den=sw*x2-x*x
           if(den.eq.0.) stop 'den'
           a=(sw*xy-x*y)/den
           b=(x2*y-x*xy)/den
         else
           den=sw*y2-y*y
           if(den.eq.0.) stop 'den'
           a=(sw*xy-y*x)/den
           b=(y2*x-y*xy)/den
           if(abs(a).lt.1.e-30) a=1.e-30
           a=1./a
           b=-b*a
         endif
c  Fill the rest of the array for contouring.
         do jjj=1,5
c         do iii=i-2,i+2
         do iii=i-2*is,i+2*is,is
         npeak=npeak+1
         pk(1,npeak)=float(iii-i)*dx
c         pk(2,npeak)=float(jjj-3)*dy
         pk(2,npeak)=float(jjj-3)*dys
         pk(3,npeak)=as(iii,jjj)
         end do
c  Find the largest maxima within the central area of the window.
c  The axis of the belt perpendicular to the linear trend will be made to
c  pass through this point.
       do iii=1,npeak
       if(abs(pk(1,iii)).lt.abs(dxs).and.abs(pk(2,iii)).lt.abs(dys))
     1 then
         xzero=pk(1,iii)
         yzero=pk(2,iii)
         go to 38
       endif
       enddo
   38 continue
         end do
         if(pans.eq.'y') then
c  Plot the data in the window.
           aa=-1./a
           bb=yzero-aa*xzero
           xd=(bb-b)/(a+1./a)
           bb1=(a+1./a)*(xd+del/(2.*sqrt(a*a+1.)))+b
           bb2=bb-(bb1-bb)
           call detour(kkk,npeak,pk,a,b,aa,bb,bb1,bb2)
         endif
  222 format(5f16.2)
      theta=-atan2(b,-b/a)
c  Extract data within the belt and solve for the squared depth.
c  (Ignoring interference terms and assuming that the xp-origin equals the
c   trend line).
      f=0.
      f2=0.
      x2f=0.
      x2f2=0.
      num=0
      sw=0.
      do 50 jj=-2,2
      y=float(jj)*dys
c      do 50 ii=-2,2
      do 50 ii=-2*is,2*is,is
      x=float(ii)*dx
      axp=as(i+ii,jj+3)
      xp=(y-yzero)*cos(theta)-(x-xzero)*sin(theta)
      yp=-(x-xzero)*cos(theta)-(y-yzero)*sin(theta)
      if(abs(yp).gt.del/2) go to 50
      w=del/(abs(xp)+del)
c      if(abs(xp).gt.del/2) go to 50
c      if(pans.eq.'y') print *, x, y, xp, yp
      num=num+1
      f=f+w*axp
      f2=f2+w*axp*axp
      x2f=x2f+w*xp*xp*axp
      x2f2=x2f2+w*xp*xp*axp*axp
      sw=sw+w
   50 continue
      if(num.lt.3) go to 40
      anum=float(num)
c      den1=anum*f2-f*f
      den1=sw*f2-f*f
c      d=(f*x2f-anum*x2f2)/den1
      d=(f*x2f-sw*x2f2)/den1
c      derr=anum/den1
      derr=sw/den1
      a1=(f2*x2f-f*x2f2)/den1
c  If the squared depth is negative, move the window.
      if(d.lt.0.) then
        derr=1.e38
        d=1.e38
        if(pans.eq.'y') then
          call setcur(0,0)
          print *,'negative depth'
          print*,'(Enter q to exit)'
c          ie(2)=3
c          call endpt(ie)
          read(5,100) dum
          if(dum.eq.'q'.or.dum.eq.'Q') then
            call setmod(2)
            close(14)
            close(15)
            go to 10
          endif
        endif
        if(nans.eq.'n') go to 40
        go to 52
      endif
c  Compute the standard error of the squared depth
      err=0
      do 51 jj=-2,2
      y=float(jj)*dys
c      do 51 ii=-2,2
      do 51 ii=-2*is,2*is,is
      x=float(ii)*dx
      axp=as(i+ii,jj+3)
      xp=(y-yzero)*cos(theta)-(x-xzero)*sin(theta)
      yp=-(x-xzero)*cos(theta)-(y-yzero)*sin(theta)
      if(abs(yp).gt.del/2) go to 51
c      if(abs(xp).gt.del/2) go to 51
      w=del/(abs(xp)+del)
c      err=err+(xp**2*axp-(a1-axp*d))**2
      err=err+w*(xp**2*axp-(a1-axp*d))**2
   51 continue
c      err=err/(anum-2.)
      err=anum*(err/sw)/(anum-2.)
      derr=sqrt(err*derr)
c  Estimate the standard error of the depth as a percent
      derr=sqrt(d+derr)
      d=sqrt(d)
      derr=100.*(derr-d)/d
      d=d-dmin
   52  xzero=xo+float(i-1)*dx+xzero
       yzero=yo+float(k-1)*dy+yzero
       dtheta=90.-theta/rad
       if(dtheta.gt.180.) dtheta=dtheta-180.
c       if(pans.eq.'y') print *, dtheta
       if(pans.eq.'y'.and.d.lt.1.e38) then
         call setcur(0,0)
c         print*,'     x:',xzero
c         call setcur(2,0)
c         print*,'     y:',yzero
c         call setcur(3,0)
c         print*,' depth:',d
c         call setcur(4,0)
c         print*,'stderr:',derr
         print *,'       x               y             depth          %s
     1tderr         strike'
c         print *,'1234567890123456789012345678901234567890123456789012'
         write(6,101) xzero, yzero, d, derr, dtheta
         print *,'(Enter q to exit)'
c          ie(2)=3
c         call endpt(ie)
         read(5,100) dum
          if(dum.eq.'q'.or.dum.eq.'Q') then
            call setmod(2)
            close(14)
            close(15)
            go to 10
          endif
       endif
       if(d.lt.0.0.and.nans.eq.'n') derr=1.e38
       if(derr.le.dlim) then
c         if(pans.ne.'y') write(6,101) xzero, yzero, d, derr, dtheta
         write(6,101) xzero, yzero, d, derr, dtheta
         write(15) sta, xzero,yzero,d,derr,dtheta,zero,zero,zero
         numdep=numdep+1
       endif
   40 continue
      do 28 j=1,4
      do 28 i=1,nc
      as(i,j)=as(i,j+1)
   28 continue
      if(iskip.gt.0) then
        do 29 j=1,iskip
        read(14,end=30) dlt
   29   continue
      endif
   30 continue
      rewind(14)
      read(14) id
      do i=1,kk
      read(14) dlt
      enddo
   60 continue
      if(pans.eq.'y') call setmod(2)
      print*,numdep,' solutions written to asdep.pst'
  101 format(5g16.6)
      stop
      end
c---------------------------------------------------------------------------
      subroutine bsort(d,kk,ierr)
      dimension d(3,kk)
      ierr=1
      if(kk.lt.2) return
      do 21 i=1,kk-1
      do 21 j=i+1,kk
      if(d(3,i)-d(3,j)) 20,21,21
   20 ax=d(1,i)
      ay=d(2,i)
      ak=d(3,i)
      d(1,i)=d(1,j)
      d(2,i)=d(2,j)
      d(3,i)=d(3,j)
      d(1,j)=ax
      d(2,j)=ay
      d(3,j)=ak
   21 continue
      ierr=0
      return
      end
c---------------------------------------------------------------------------
      subroutine interp(y1,y2,y3,d,xmax,ymax,ierr)
c
c     Subroutine interp fits a parabola through ordinants y1,y2 and y3, assumed
c  to be evenly spaced, and reports the x and y coordinate of the parabola's
c  maximum.  Note...y1 must be .lt. y2 and y3 must be .lt. y2.
c
      ierr=1
      if(y1.lt.y2.and.y3.lt.y2) then
        ierr=0
      else
        return
      endif
      if(d.eq.0.) stop 'd'
      a=0.5*(y1-2.*y2+y3)/(d*d)
      b=0.5*(y3-y1)/d
      c=y2
      if(a.eq.0.)then
        ierr=1
        return
      end if
      if(a.eq.0.) stop 'a'
      xmax=-0.5*b/a
      ymax=a*xmax**2+b*xmax+c
      return
      end
c---------------------------------------------------------------------------
      subroutine detour(n1,n,pk,a,b,aa,bb,bb1,bb2)
c
c  Contouring of scattered 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 pk(3,100),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),
     +          ipn(12)
c     +          ie(10), ipn(12)
      integer itetr(2401,3), istack(2401), ktetr(50,2), itemp(3,2)
      character lab*8
      data itemp,xpnt/1,1,2,2,3,3,-1.,5.,-1.,-1.,-1.,5.,2.,2.,18./
      data ipn/100,900,300,1100,200,1000,1400,600,400,1200,500,1300/
  100 format(a)
      iplotr=9
      do 1 i=1,3
      itetr(1,i)=i
      tetr(1,i)=xpnt(i,3)
      do 1 j=1,3
    1 xyz(i,j)=xpnt(i,j)
      do 2 i=2,65
    2 istack(i)=i
      size=1.0
      xyz(4,1)=pk(1,1)
      xyz(4,2)=pk(2,1)
      xyz(4,3)=pk(3,1)
      xmin=xyz(4,1)
      xmax=xyz(4,1)
      ymin=xyz(4,2)
      ymax=xyz(4,2)
      zmin=xyz(4,3)
      zmax=xyz(4,3)
c      print*,n
      do 3 i=5,n+3
c      print*,i
      xyz(i,1)=pk(1,i-3)
      xyz(i,2)=pk(2,i-3)
      xyz(i,3)=pk(3,i-3)
      if(xyz(i,1).lt.xmin) xmin=xyz(i,1)
      if(xyz(i,1).gt.xmax) xmax=xyz(i,1)
      if(xyz(i,2).lt.ymin) ymin=xyz(i,2)
      if(xyz(i,2).gt.ymax) ymax=xyz(i,2)
      if(xyz(i,3).lt.zmin) zmin=xyz(i,3)
      if(xyz(i,3).gt.zmax) zmax=xyz(i,3)
    3 continue
    4 n=n+3
      n1=n1+3
       conti=(zmax-zmin)/10.
       if(conti.eq.0.) stop 'conti'
      contm=zmin-amod(zmin,conti)+conti
      lc=int((zmax-contm)/conti)
      xr=xmax-xmin
      datax=xr
      yr=ymax-ymin
      if(yr.gt.datax) datax=yr
c  normalize the data
      if(datax.eq.0.) stop 'datax'
      do 5 i=4,n
      xyz(i,1)=(xyz(i,1)-xmin)/datax
    5 xyz(i,2)=(xyz(i,2)-ymin)/datax
      isp=1
      id=2
      do 13 nuc=4,n
      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 i=1,3
      l1=itemp(i,1)
      l2=itemp(i,2)
      if(km.le.0) go to 8
      kmt=km
      do 7 j=1,kmt
      if(itetr(jt,l1).ne.ktetr(j,1)) go to 7
      if(itetr(jt,l2).ne.ktetr(j,2)) go to 7
      km=km-1
      if(j.gt.km) go to 9
      do 6 k=j,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 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 tetr(kt,*)
      do 11 jz = 1,2
      i2 = ktetr(i,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)
c      if(dd.eq.0) stop 'dd'
      if(abs(dd).lt.1.e-10) go to 56
      if(kt.lt.1) kt=1
      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(i,1)
      itetr(kt,2) = ktetr(i,2)
   12 itetr(kt,3) = nuc
   13 isp = isp + 2
   56 call pltset(iplotr,xboard,yboard,isl)
      dxp(1)=0.0
      dxp(2)=xr/datax
      dyp(1)=0.0
      dyp(2)=yr/datax
c      yp(1)=yboard-0.2
      yp(1)=yboard-1.0
      if(yr.eq.0.) stop 'yr'
      xp(1)=yp(1)*xr/yr
      if(xp(1).gt.xboard-1.2) then
        xp(1)=xboard-1.2
        if(xr.eq.0.) stop 'xr'
        yp(1)=xp(1)*yr/xr
      endif
      xp(2)=0.
      xp(3)=0.1
      xp(4)=xboard
      yp(2)=0.
      yp(3)=0.1
      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
      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
      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
      do 20 i=4,n
      call vchar(xyz(i,1),xyz(i,2),3,1,1500,.1,0.,0.,0.)
   20 continue
      do 21 i=4,n1
      call vchar(xyz(i,1),xyz(i,2),42,1,1500,.1,0.,0.,0.)
   21 continue
      ipold=0
      y=0.
      do 22 i=1,lc
      ip=1+i*12/(lc+1)
      cont=contm+(i-1)*conti
      if(ip.ne.ipold) then
      write(lab(1:8),'(g8.2)') cont
      y=y+.4
      call vchar(xp(1)+0.5,y,lab,8,ipn(ip)+3,.1,0.,0.,0.)
      endif
      ipold=ip
   22 continue
      b=(a*xmin+b-ymin)/datax
      x(1)=0.
      y(1)=b
      x(2)=1.
      y(2)=a+b
      call line(x,y,2,0,0)
      bb=(aa*xmin+bb-ymin)/datax
      y(1)=bb
      y(2)=aa+bb
      call line(x,y,2,0,0)
      bb1=(aa*xmin+bb1-ymin)/datax
      y(1)=bb1
      y(2)=aa+bb1
      call line(x,y,2,0,1)
      bb2=(aa*xmin+bb2-ymin)/datax
      y(1)=bb2
      y(2)=aa+bb2
      call line(x,y,2,0,1)
c   21 call endpt(ie)
c   21 continue
      return
      end
