c  Program PLOTDIP
c
c  Plotting of depth, dip, etc. estimates from LOCDEP
c
c  written by J. Phillips
c
c      dimension xyz(1203,3)
c      dimension isl(5), dxp(2), dyp(2), xp(4), yp(4), x(100), y(100),
c     +          ie(10), ipn(12), post(6), flx(100,2)
      dimension isl(5), dxp(2), dyp(2), xp(4), yp(4),
     +          ie(10), ipn(12), post(6)
      character ifile*80, label*21,sta*8
      data ddval/1.e30/
      data ipn/100,900,300,1100,200,1000,1400,600,400,1200,500,1300/
      rad=3.1415927/180.
      print *,'input post file: '
      read(5,100) ifile
  100 format(a)
c      open(10,file='edep3d.pst',form='unformatted',status='old')
      open(10,file=ifile,form='unformatted',status='old')
1     print*,'1=depth, 2=error, 3=strike, 4=dip_azimuth, 5=dip, 6=susc'
      print *,'enter z-field for coloring (1 to 6):'
      read(5,*) iz
      xmin=1.e+37
      xmax=-1.e+37
      ymin=1.e+37
      ymax=-1.e+37
      zmin=1.e+37
      zmax=-1.e+37
      ifile(1:1)=' '
      nl=0
c      print *,'enter lxfile (or return):'
c      read(5,100) ifile
cc      print *,ifile
c      if(ifile(1:1).eq.' ') go to 9
c      open(11,file=ifile,form='formatted',status='old')
c      do 10 i=1,100
c      read(11,*,end=8) flx(i,1),flx(i,2)
cc      print *, flx(i,1), flx(i,2)
c      if(flx(i,1).ge.1.e+37) go to 10
c      if(flx(i,1).lt.xmin) xmin=flx(i,1)
c      if(flx(i,1).gt.xmax) xmax=flx(i,1)
c      if(flx(i,2).lt.ymin) ymin=flx(i,2)
c      if(flx(i,2).gt.ymax) ymax=flx(i,2)
c   10 continue
    8 nl=i-1
c    9 print *,'enter iplotr:'
c      read(5,*) iplotr
c  read the data
c      read(10) lab(4),xyz(4,1),xyz(4,2),post
c      read(10) xyz(4,1),xyz(4,2),xyz(4,3)
c      xmin=xyz(4,1)
c      xmax=xyz(4,1)
c      ymin=xyz(4,2)
c      ymax=xyz(4,2)
c      zmin=xyz(4,3)
c      zmax=xyz(4,3)
c      do 3 i=4,1203
c      read(10,end=4) lab(i),xyz(i,1),xyz(i,2),post
      print *,'maximum percent error in depth:'
      read(5,*) pcterr
    3 read(10,end=4) sta,x,y,post
      if(e.gt.pcterr) go to 3
      if(x.lt.xmin) xmin=x
      if(x.gt.xmax) xmax=x
      if(y.lt.ymin) ymin=y
      if(y.gt.ymax) ymax=y
      if(post(iz).lt.zmin) zmin=post(iz)
      if(post(iz).gt.zmax) zmax=post(iz)
      go to 3
c    3 continue
c      print *,'warning - input file may contain more than 1200 points'
c    4 n=i-1
    4 continue
      xmax = xmax+1.
      xmin = xmin-1.
      ymax = ymax+1.
      ymin = ymin-1.
      print *,'xmin = ',xmin,' xmax = ',xmax
      print *,'ymin = ',ymin,' ymax = ',ymax
      print *,'zmin = ',zmin,' zmax = ',zmax
c      print *, 'redefine z-range?:'
c      read(5,100) tri
      zmin0=zmin
      zmax0=zmax
c      if(tri.ne.'y'.and.tri.ne.'Y') go to 6
c      print *, 'enter new zmin, zmax:'
c      read(5,*) zmin,zmax
c      do 7 i=4,n
c      if(xyz(i,3).lt.zmin) xyz(i,3)=zmin
c      if(xyz(i,3).gt.zmax) xyz(i,3)=zmax
c    7 continue
c    6 print *, 'enter color interval (negative to reverse color scale):'
c      read(5,*) conti
      conti=(zmax-zmin)/12.
   50 print '(a,1p,g15.6,a\)', ' Enter contour interval [',conti,']:'
      call getfloat(conti,ierr)
      if(ierr.ne.0) go to 50
      aplotr=10
      print '(a\)',' Enter iplotr (5=HP,8=CGA,9=EGA,10=VGA) [10]:'
c      read(5,*) iplotr
      call getfloat(aplotr,ierr)
      iplotr=aplotr
      if(conti.lt.0.0) then
        do 7 i=1,6
        isav=ipn(i)
        ipn(i)=ipn(13-i)
        ipn(13-i)=isav
    7   continue
        conti=-conti
      endif
      if(conti.eq.0.) go to 13
      zmin=int(zmin0/conti)*conti
      zmax=int(zmax0/conti)*conti
   10 if(zmin.ge.zmin0) go to 11
      zmin=zmin+conti
      go to 10
   11 if(zmax.le.zmax0) go to 12
      zmax=zmax-conti
      go to 11
   12 if(zmax.le.zmin) stop 'unacceptable range'
c      contm=zmin-amod(zmin,conti)+conti
c      contm=amod(zmin,conti)
c      print*,'zmin=',zmin,' amod(zmin,conti)=',contm
c      stop
c      contm=zmin-amod(zmin,conti)
      contm=zmin
      lc=int(.000001+(zmax-contm)/conti)
      if(zmin.ne.zmin0) lc=lc+1
c      if(zmax.ne.zmax0) lc=lc+1
13    xr=xmax-xmin
      datax=xr
      yr=ymax-ymin
      if(yr.gt.datax) datax=yr
c  normalize the data
c      do 5 i=4,n
c      xyz(i,1)=(xyz(i,1)-xmin)/datax
c    5 xyz(i,2)=(xyz(i,2)-ymin)/datax
      isp=1
      id=2
      call pltset(iplotr,xboard,yboard,isl)
      dxp(1)=-0.1
      dxp(2)=xr/datax+0.1
      dyp(1)=-0.1
      dyp(2)=yr/datax+0.1
      yp(1)=yboard-0.2
      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.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
      rewind(10)
c      do 20 i=4,n
      if(conti.eq.0.) go to 21
   20 read(10,end=21) sta,x,y,post
      x=(x-xmin)/datax
      y=(y-ymin)/datax
      e=post(2)
      s=post(3)
      s=-rad*s
c      t=post(4)
c      t=rad*t
c      d=post(5)
c      d=cos(rad*d)
      if(post(iz).ge.ddval) go to 20
c      print*,z,contm,conti
      jc=int(1.0+(post(iz)-contm)/conti)
      if(jc.gt.lc) jc=lc
cc      if(z.lt.zmin) jc=0
cc      if(z.gt.zmax) z=zmax
      ip=1+jc*12/(lc+1)
cc      print*,ip,jc,lc
cc      if(lab(i).eq.'contact ') then
cc        call vchar(xyz(i,1),xyz(i,2),'C',1,ipn(ip),.1,0.,0.,0.)
cc      else if(lab(i).eq.'sill    ') then
cc        call vchar(xyz(i,1),xyz(i,2),'S',1,ipn(ip),.1,0.,0.,0.)
cc      else if(lab(i).eq.'pipe    ') then
cc        call vchar(xyz(i,1),xyz(i,2),'P',1,ipn(ip),.1,0.,0.,0.)
cc      else if(lab(i).eq.'dipole  ') then
cc        call vchar(xyz(i,1),xyz(i,2),'D',1,ipn(ip),.1,0.,0.,0.)
cc      endif
cc        print*,ip,ipn(ip),jc,lc
        if(e.le.pcterr) then
          size=0.05+0.25*(post(iz)-zmin)/(zmax-zmin)
          call vchar(x,y,3,1,ipn(ip),
     1               size,0.,0.,0.)
19      call vchar(x,y,108,1,1500,
     1             .175,s,0.,0.)
        endif
c        ipend=1500
c        if(conti.eq.0.) ipend=1200
c        call line(x,y,1,0,ipend)
c        x=x+.0175*d*sin(t)
c        y=y+.0175*d*cos(t)
c        call line(x,y,1,1,ipend)
      go to 20
c   20 continue
   21 rewind(10)
   30 read(10,end=31) sta,x,y,post
      x=(x-xmin)/datax
      y=(y-ymin)/datax
      e=post(2)
      s=post(3)
      s=-rad*s
      t=post(4)
      t=rad*t
      d=post(5)
      d=cos(rad*d)
c      if(conti.eq.0.) go to 19
c      if(post(iz).ge.ddval) go to 19
cc      print*,z,contm,conti
c      jc=int(1.0+(post(iz)-contm)/conti)
c      if(jc.gt.lc) jc=lc
cc      if(z.lt.zmin) jc=0
cc      if(z.gt.zmax) z=zmax
c      ip=1+jc*12/(lc+1)
cc      print*,ip,jc,lc
cc      if(lab(i).eq.'contact ') then
cc        call vchar(xyz(i,1),xyz(i,2),'C',1,ipn(ip),.1,0.,0.,0.)
cc      else if(lab(i).eq.'sill    ') then
cc        call vchar(xyz(i,1),xyz(i,2),'S',1,ipn(ip),.1,0.,0.,0.)
cc      else if(lab(i).eq.'pipe    ') then
cc        call vchar(xyz(i,1),xyz(i,2),'P',1,ipn(ip),.1,0.,0.,0.)
cc      else if(lab(i).eq.'dipole  ') then
cc        call vchar(xyz(i,1),xyz(i,2),'D',1,ipn(ip),.1,0.,0.,0.)
cc      endif
cc        print*,ip,ipn(ip),jc,lc
c        if(e.le.pcterr) then
c          size=0.05+0.25*(post(iz)-zmin)/(zmax-zmin)
c          call vchar(x,y,3,1,ipn(ip),
c     1               size,0.,0.,0.)
c        endif
c29      call vchar(x,y,108,1,1500,
c     1             .175,s,0.,0.)
29      continue
c        ipend=1500
c        if(conti.eq.0.) ipend=1200
c        call line(x,y,1,0,ipend)
c        x=x+.0175*d*sin(t)
c        y=y+.0175*d*cos(t)
c        call line(x,y,1,1,ipend)
      go to 30
c
   31 ipold=1
      y=0.
      cont=zmin0
      write(label(1:21),'(g9.3)') cont
      call vchar(xp(1)+0.75,y+.2,label,9,3,.1,0.,0.,0.)
c      print*,contm
      do 22 i=1,lc+1
      ip=1+i*12/(lc+1)
      if(ip.ne.ipold) then
        cont=contm+(i-1)*conti
        if(zmin.eq.zmin0) cont=cont+conti
        if(cont.gt.zmax) cont=zmax0
c        cont1=cont+conti
c       cont=contm+(float(i)-.5)*conti
c        write(label(1:21),'(g9.3,3h - ,g9.3)') cont,cont1
       write(label(1:21),'(g9.3)') cont
        y=y+.4
        size=0.05+0.25*(cont-zmin)/(zmax-zmin)
        call vchar(xp(1)+0.5,y,3,1,ipn(ipold)+3,
     1             size,0.,0.,0.)
c        call vchar(xp(1)+0.75,y,label,21,ipn(ipold)+3,.1,0.,0.,0.)
        call vchar(xp(1)+0.75,y+.2,label,9,3,.1,0.,0.,0.)
c        cont=cont1
        ipold=ip
      endif
c      ipold=ip
   22 continue
c      if(nl.ne.0) then
c        j=0
c   13   j=j+1
c        if(j.ge.nl-1) go to 23
c        do 12 i=1,nl-1
c        if(flx(j,1).ge.1.e+37) go to 11
c        x(i)=(flx(j,1)-xmin)/datax
c        y(i)=(flx(j,2)-ymin)/datax
c        j=j+1
c   12   continue
c   11   call line(x,y,i-1,0,0)
c        go to 13
c      endif
   23 call endpt(ie)
      rewind(10)
      go to 1
      stop
      end
c******************************************************************************
      subroutine getfloat(default,ierr)
      character string*80, fmt*40
      ierr=0
      read(*,'(a)') string
      leng=len_trim(string)
      if(leng.ne.0) then
        lenm=leng-index(string(1:leng),'.')
        if(lenm.eq.leng) lenm=0
        write(fmt,'(a,i3,a,i3,a)') '(g',leng,'.',lenm,')'
        read(string(1:leng),fmt,err=990) default
      endif
      go to 999
  990 ierr=1
  999 return
      end

