c es.for
c  solves for c only; assigns z=ar; fits data minus epsilon_2.
c  graphics added.
c  input from post (.pos) data files.
c  depth calculated from station h
c  input data converted to z,h positive down system
c  To compile: fl es.for graphics.lib
c  Lin Cordell, 13 Feb 91.
c      include 'fgraph.fi'
c      include 'fgraph.fd'
      character*50 ifile,ofile,ofile2
      character*1 ans
      dimension id(14),pgm(2)
      dimension g(10000),rad(2),ct(2),v(6)
      dimension x(10000),y(10000),h(10000),f(10000)
      dimension z(10000),c(10000)
      data dval/1.0e+38/
c  Read data.
      print*,'Enter input post file name:'
      print 1
1     format(' *',$)
      read(5,800)ifile
800   format(a)
      open(10,file=ifile,form='unformatted',status='old')
      print*,'POST FIELD   GRAVITY             MAGNETICS'
      print*,'    1        free-air anomaly    total field residual'
      print*,'    2        Bouguer  anomaly    total field'
      print*,'    3        elevation           terrain clearance'
      print*,'    4        inner zone TC       barometric altitude'
      print*,'    5        outer zone TC       fiducial number'
      print*,'    6        observed gravity    julien date'
      print*,'Enter post fields for data, elevation:'
      print 1
      read(5,*) ifld, ielv
79    print*,'Elevation in feet, meters, or km?'
      print 1
      read(5,800)ans
      efact=0
      if(ans.eq.'f'.or.ans.eq.'F')efact=.0003048
      if(ans.eq.'m'.or.ans.eq.'M')efact=.001
      if(ans.eq.'k'.or.ans.eq.'K')efact=1.
      if(efact.eq.0.)go to 79
      k=0
      xmin=1.e20
      xmax=-1.e20
      ymin=1.e20
      ymax=-1.e20
      hmin=1.e20
      hmax=-1.e20
      fmin=1.e20
      fmax=-1.e20
40    k=k+1
      if(k.gt.10000) stop 'to many data points'
      read(10,end=41)id(1),id(2),x(k),y(k),v
      f(k)=v(ifld)
      h(k)=v(ielv)
c      read(10,end=41)id(1),id(2),x(k),y(k),v,f(k),h(k),v,v,v
      if(x(k).gt.xmax)xmax=x(k)
      if(x(k).lt.xmin)xmin=x(k)
      if(y(k).gt.ymax)ymax=y(k)
      if(y(k).lt.ymin)ymin=y(k)
      h(k)=-h(k)*efact
c  Elevation converted to +down system, in km.
      if(h(k).gt.hmax)hmax=h(k)
      if(h(k).lt.hmin)hmin=h(k)
      if(f(k).gt.fmax)fmax=f(k)
      if(f(k).lt.fmin)fmin=f(k)
      goto 40
41    kmax=k-1
      close(10)
      write(6,42)kmax,xmin,xmax,ymin,ymax,hmin,hmax,fmin,fmax
42    format(i5,' Data points.'/' xmin,max = ',2e12.3,/' ymin,max'
     &'= ',2e12.3,/' hmin,max = ',2e12.3,' km'/' fmin,max = ',2e12.3)
      print*,'Enter reduction level (- if up,in km):'
      print 1
      read*,hcalc
c  Phony data.
      print*,' Want to include phony data?'
      read(5,800)ans
      if(ans.ne.'y'.and.ans.ne.'Y')go to 43
      print*,'Enter input phony-data post file name:'
      print 1
      read(5,800)ifile
      open(10,file=ifile,form='unformatted',status='old')
      k=kmax
44    k=k+1
      if(k.gt.10000) stop 'too many data points'
c      read(10,end=45)id(1),id(2),x(k),y(k),v,f(k),h(k),v,v,v
      read(10,end=45)id(1),id(2),x(k),y(k),v
      f(k)=v(ifld)
      h(k)=v(ielv)
      h(k)=-h(k)*efact
c  Elevation converted to +down system, in km.
      if(h(k).gt.hmax)hmax=h(k)
      if(h(k).lt.hmin)hmin=h(k)
      if(f(k).gt.fmax)fmax=f(k)
      if(f(k).lt.fmin)fmin=f(k)
      goto 44
45    print*, k-1-kmax,' phony data points'
      kmax=k-1
      print*, kmax,' total data points'
      write(6,42)kmax,xmin,xmax,ymin,ymax,hmin,hmax,fmin,fmax
      close(10)
43    continue
c  Begin computations.
46    print*
      print*,' a = ratio of source depth to nearest-neighbor distance'
      print*,'     (eg. 1.4).'
      print*,' nmax = maximum iterations (eg. 200).'
      print*,' eps = desired maximum error (eg. 1.0)'
      print*,' eps2 = desired fit to each data point (eg. 0.1)'
      print*
      print*,' Enter a,nmax,eps,eps2:'
      print 1
      read*, a,nmax,eps,eps2
      if(eps.lt.0.or.eps2.lt.0.or.eps2.gt.eps)then
        print*, 'Error: reguire eps >= eps2 >= 0.'
        go to 46
      endif
c  Parameter a | z=a*r.
c  In this version, fit is to f-eps2, rather than to f.
c   and z is figured from station h, but z > hmax to assure
c   z not above the lowest station.
c  Define datum, neighbor, and z(k) array.
      fsum=0
      do 51 k=1,kmax
      if(mod(k,50).eq.0) write(6,50) k
      fsum=fsum+f(k)
      rsqmin=1.e20
      do 55 j=1,kmax
      if(j.eq.k)go to 55
      xsq=(x(k)-x(j))**2
      ysq=(y(k)-y(j))**2
      rsq=xsq+ysq
      if(rsq.lt.rsqmin)rsqmin=rsq
55    continue
      z(k)=a*sqrt(rsqmin)+h(k)
      if(z(k).lt.hmax)z(k)=hmax
51    continue
50    format(1h+,'working on data point ',i5)
      datum=fsum/kmax
c  Remove datum.
      do 56 k=1,kmax
56    f(k)=f(k)-datum
c  Read grid specs and mask.
      print*,'Enter output grid file name:'
      print 1
      read(5,800)ofile
      print*,'Enter id:'
      print 1
      read(5,4)(id(i),i=1,14)
4     format(14a4)
      pgm(1)='e_so'
      pgm(2)='urce'
      print*,'Enter x0,y0,dx,dy,ncol,nrow'
      print 1
      read*,x0,y0,dx,dy,ncol,nrow
c  Mask output grid.
      do 400 i=1,ncol
400   g(i)=dval
      lrec=4*ncol
      open(12,status='scratch',access='direct',
     & recl=lrec,form='unformatted')
      do 401 j=1,nrow
401   call rrowio(g,ncol,12,j,2)
      del=0.5*(dx+dy)
      denom=del*a
      do 410 k=1,kmax
      mask=(z(k)-h(k))/denom
c  Note masking range = 2*distance to nearest neighbor.
      i=((x(k)-x0)/del)+1.5
      i1=i-mask
      if(i1.lt.1)i1=1
      i2=i+mask
      if(i2.gt.ncol)i2=ncol
      j=((y(k)-y0)/del)+1.5
      j1=j-mask
      if(j1.lt.1)j1=1
      j2=j+mask
      if(j2.gt.nrow)j2=nrow
      do 409 j=j1,j2
      call rrowio(g,ncol,12,j,1)
      do 408 i=i1,i2
408   g(i)=0.0e0
409   call rrowio(g,ncol,12,j,2)
410   continue
c  Iterations (kk denotes the source).
      iplotr= ask_int('Enter plot device code (0=no plot,5=HP,8=CGA,9=EG
     &A,10=VGA): ')
c      print*,'Hit ENTER to continue, and '
c      print*,' hit ENTER again to erase the graph and continue.'
c      read (*,*)
      n=0
c      call graphicsmode()
c      call frame()
220   n=n+1
      af0=0
      do 210 k=1,kmax
      af=abs(f(k))
      if(af.lt.af0)go to 210
      kk=k
      af0=af
210   continue
      if(n.eq.1)afmax=af0
      if(n.eq.1.and.iplotr.ne.0) then
        call rplot(0.,float(nmax), 0., afmax, iplotr)
        rad(1)=0.
        rad(2)=nmax
        ct(1)=eps
        ct(2)=eps
        call curv(rad, ct, 2, 300)
      endif
c      call plotpoint(n,af0,afmax)
      if(iplotr.ne.0) call vchar(float(n),af0,6,1,601,.08,0.,0.,0.)
      if(af0.le.eps)go to 300
c  Fit to f +- eps.
      if(f(kk).gt.0)then
        cnst=(f(kk)-eps2)*abs(z(kk)-h(kk))
      else
        cnst=(f(kk)+eps2)*abs(z(kk)-h(kk))
      endif
      c(kk)=c(kk)+cnst
c  Remove this source.
      do 211 k=1,kmax
      rsq=(x(k)-x(kk))**2+(y(k)-y(kk))**2+(z(kk)-h(k))**2
c      print*, k,f(k),cnst,rsq
      if(rsq.eq.0.0) then
        print*,'Duplicate data points at x=',x(k),' y=',y(k)
        print*,'Run G_SCREEN before ES!'
        stop
      endif
211   f(k)=f(k)-(cnst/sqrt(rsq))
      if(n.lt.nmax)go to 220
c  Iterations completed.
300   if(n.lt.nmax)nmax=n
c      call endgraphics()
      call endit()
      write(6,301)nmax,af0,datum
301   format(/,' After ',i4,' iterations, the maximum error was ',
     1       e14.7,'.'/,
     2       ' A datum of ',e14.7,' was subtracted from the data.')
      print*,'Make a note of datum value if you intend to use ES_CK,'
      print*,'ES_COF, or ES_FWD.'
      print*
c  Collapse source array and output sources file.
      print*,'Enter output sources file name:'
      print 1
      read(5,800)ofile2
      open(13,file=ofile2,form='unformatted',status='unknown')
      k=0
      j=0
310   k=k+1
311   j=j+1
      if(j.gt.max(kmax,nmax))go to 320
      if(c(j).eq.0.0)go to 311
      x(k)=x(j)
      y(k)=y(j)
      h(k)=h(j)
      c(k)=c(j)
      z(k)=z(j)
c      write(13)id(1),id(2),x(k),y(k),z(k),c(k),v,v,v,v
      v(1)=z(k)
      v(2)=c(k)
      write(13)id(1),id(2),x(k),y(k),v
      go to 310
320   kmax=k-1
      close(13)
      print*,kmax,' sources.'
c  Calculate grid.
      nz=1
      open(11,file=ofile,status='unknown',form='unformatted')
      write(11)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
c  Note redefinition of z.
      do 258 k=1,kmax
258   z(k)=(z(k)-hcalc)**2
      do 75 j=1,nrow
      yg=y0+(j-1)*dy
      call rrowio(g,ncol,12,j,1)
      do 74 i=1,ncol
      if(g(i).ne.0.0e0)go to 74
      xg=x0+(i-1)*dx
      g(i)=datum
      do 73 k=1,kmax
      rsq=(xg-x(k))**2+(yg-y(k))**2+z(k)
73    g(i)=g(i)+(c(k)/sqrt(rsq))
74    continue
      call rowio(g,yg,1,ncol,11,2)
75    continue
999   continue
      close(10)
      close(11)
      close(12)
      stop
      end
      subroutine rowio(a,y,nz,ncol,ld,key)
      dimension a(nz,ncol)
      go to (1,2),key
1     read(ld) y,a
      go to 90
2     write(ld)y,a
90    return
      end
      subroutine rrowio(a,ncol,ld,jrec,key)
c  To read/write row from random-access file.
      dimension a(ncol)
      go to (1,2),key
1     read(ld,rec=jrec) a
      go to 90
2     write(ld,rec=jrec) a
90    return
      end
c      subroutine graphicsmode()
c      include 'fgraph.fd'
c      integer*2 dummy,maxx,maxy
c      record /videoconfig/ screen
c      common maxx,maxy
cc Find graphics configuration.
c      call getvideoconfig(screen)
c      select case(screen.adapter)
c      case($cga)
c       dummy=setvideomode($hresbw)
cc       print*,'CGA'
c      case($ocga)
c       dummy=setvideomode($orescolor)
c      case($ega,$oega)
c       if(screen.monitor.eq.$mono)then
c        dummy=setvideomode($eresnocolor)
cc        print*,'EGA-MON0'
c       else
c        dummy=setvideomode($erescolor)
cc        print*,'EGA-COLOR'
c       endif
c      case($vga,$ovga,$mcga)
cc       dummy=setvideomode($vres256color)
c       dummy=setvideomode($vres16color)
cc       print*,'VGA'
c      case($hgc)
c       dummy=setvideomode($hercmono)
cc       print*,'HERCULES'
c      case default
c       dummy=0
c      end select
c      if(dummy.eq.0)stop 'Error: cannot set graphics mode.'
cc Find min/max screen dimensions.
c      call getvideoconfig(screen)
c      maxx=screen.numxpixels-1
c      maxy=screen.numypixels-1
cc      print*,'SCREEN',maxx,maxy
c      end
c      integer*2 function newx(x)
c Function converts x from [0,1000] to [0,maxx] range.
c      integer*2 x,maxx,maxy
c      real*4 xtmp
c      common maxx,maxy
c      xtmp=maxx*.001
c      xtmp=xtmp*x+.5
c      newx=xtmp
c      end
c      integer*2 function newy(y)
cc Function converts y from [0,1000] to [0,maxy] range.
c      integer*2 y,maxx,maxy
c      real*4 ytmp
c      common maxx,maxy
c      ytmp=maxy*.001
c      ytmp=ytmp*y+.5
c      newy=ytmp
c      end
c      subroutine frame()
c      include 'fgraph.fd'
c      external newx,newy
c      integer*2 dummy,newx,newy,maxx,maxy
c      record /xycoord/ xy
c      common maxx,maxy
cc Border.
c      dummy=setcolor(int4(4))
c      dummy=rectangle($gborder,0,0,maxx,maxy)
c      dummy=setcolor(int4(10))
c      call setvieworg(0,newy(int2(950)),xy)
cc 0 line.
cc      call setlinestyle(#AA3C)
c      call setlinestyle(0)
c      call moveto(0,newy(int2(0)),xy)
c      dummy=lineto(newx(int2(1000)),0)
c      end
c      subroutine plotpoint(i,f,fmax)
c      include 'fgraph.fd'
c      integer*2 dummy,newx,newy,ix,iy,ii
c      imax=1000
c      ii=int2(i)
c      ff=-(f/fmax)*1000
c      ix=newx(ii)
c      iy=newy(int2(ff))
c      dummy=setpixel(ix,iy)
c      return
c      end
c      subroutine endgraphics()
c      include 'fgraph.fd'
c      integer*2 dummy
c      read(*,*)            !Wait for ENTER
c      dummy=setvideomode($defaultmode)
c      end
      subroutine rplot(rmin, rmax, hmin, hmax, iplotr)
      dimension xp(4), yp(4), dxp(2), dyp(2)
      dimension x(n), y(n)
c      integer str(n)
      call pltset(iplotr, xp(4), yp(4), 1)
      xp(4) = min(xp(4),10.)
      yp(4) = min(yp(4),8.)
      xp(4) = xp(4) - .25
      yp(4) = yp(4) - .25
      dxp(1) = rmin
      dxp(2) = (rmax + .01) - mod(rmax,.01)
      dyp(1) = hmin
      dyp(2) = hmax
      xp(1) = xp(4) - 1.5
      xp(2) = 0.
      xp(3) = 1.25
      yp(1) = yp(4) - 1.25
      yp(2) = 0.
      yp(3) = 1.
      call scale(dxp, dyp, xp, yp, 4, ier)
      if (ier .ne. 0) goto 90
c      goto (51, 52, 53), ityp
c   51 call yaxis(dyp, dxp, yp, 0.1, 2, .12, '(f4.1)', 4)
c      goto 54
c   52 call yaxis(dyp, dxp, yp, 1.0, 2, .12, '(f4.0)', 4)
c      goto 54
   53 call yaxis(dyp, dxp, yp, 0.5, 2, .12, '(f4.1)', 4)
c   54 goto (55, 55, 56), ityp
c   55 call xaxis(dxp, dyp, xp, 0.1, 2, .12, '(f4.1)', 4)
c      goto 57
   56 call xaxis(dxp, dyp, xp, 5.0, 2, .12, '(f4.0)', 4)
   57 call neatl
c      goto (58, 58, 59), ityp
   58 call vchar((dxp(1) + dxp(2)) * .5, dyp(1), 'iterations', 10, 2,
     &.12, 0., -.54, -.64)
c      goto 60
c   59 call vchar((dxp(1) + dxp(2)) * .5, dyp(1), 'azimuth', 7, 2, .12,
c     &0., -.54, -.54)
c   60 goto (61, 62, 62), ityp
   61 call vchar(dxp, (dyp(1) + dyp(2)) * .5, 'max error', 9, 2, .12,
     &1.5706, -.48, .79)
c      goto 64
c   62 call vchar(dxp, (dyp(1) + dyp(2)) * .5, 'log power', 9, 2, .12,
c     &1.5706, -.48, .54)
   64 return 
   90 stop 
      entry curv(x, y, n, m)
      call line(x, y, n, 0, m)
c      call line(x, y, 1, 0, m)
c      do i=2,n-1
c      call line(x(i),y(i),2,1,m)
c      end do
      return 
c      entry text(str,n,off)
c      call vchar((dxp(1) + dxp(2)) * .1, dyp(2) - off, str, n, 2, .10,
c     & 0., 0., 0.)
c      call home
c      return      
      entry endit()
      call endpt(ie)
      return 
      end
      function ask_int(request)
      character request*(*)
    2 write(unit=6, fmt=100) request
  100 format(1x,a,$)
      read(unit=5, fmt=*, err=1) ask_int
      return 
    1 write(unit=6, fmt=101) 
  101 format(/,20h **ERROR - try again)
      call bell
      goto 2
      end
      subroutine bell()
      character ding*1
      ding = char(7)
      write(unit=*, fmt=*) ding
      return 
      end
