c es_ck.for
c
c  To calculate field of a point_sources file at the original
c  data, to check final fit to actual (not including phony) data.
c  Input from post (.pos) data and sources files.
c  Lin Cordell, 7 Feb 91.
$LARGE
      character*50 sfile,dfile
      character*1 ans
      character*4 id1,id2
      dimension xs(5000),ys(5000),zs(5000),c(5000),v(6)
      dimension x(10000),y(10000),h(10000),f(10000)
c      data dval/0.1701412e+39/
      print*,'Enter equivalent sources file name:'
      print 1
1     format(1x,'*'$)
      read(5,800)sfile
      open(10,file=sfile,status='old',form='unformatted')
      k=1
c400   read(10,end=401)id1,id2,xs(k),ys(k),zs(k),c(k),v,v,v,v
400   read(10,end=401)id1,id2,xs(k),ys(k),v
      zs(k)=v(1)
      c(k)=v(2)
      k=k+1
      if(k.gt.5000) stop 'too many sources'
      go to 400
401   continue
      kmax=k-1
      close(10)
      print*,kmax,' sources.'
      print*,'Enter data post file name:'
c  input data file must have elevations that are positive upward.
      print 1
      read(5,800)dfile
800   format(a)
      open(11,file=dfile,status='old',form='unformatted')
      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
      hmin=1.e20
      i=1
c402   read(11,end=403)id1,id2,x(i),y(i),v,f(i),h(i),v,v,v
402   read(11,end=403)id1,id2,x(i),y(i),v
      f(i)=v(ifld)
      h(i)=v(ielv)
      if(h(i).lt.hmin)hmin=h(i)
      i=i+1
      if(i.gt.10000) stop 'too many data points'
      go to 402
403   imax=i-1
      print*,'Enter datum (from ES program) to be added to '
      print*,'calculated field:'
      print 1
      read*,datum
79    print*,'Data file elevations in feet, meters, or km?'
      print 1
      read(5,800)ans
      ef=0
      if(ans.eq.'f'.or.ans.eq.'F')ef=.0003048
      if(ans.eq.'m'.or.ans.eq.'M')ef=.001
      if(ans.eq.'k'.or.ans.eq.'K')ef=1.
      if(ef.eq.0.)go to 79
      sum=0
      ssq=0
      rmax=0
      count=0
      do 410 i=1,imax
      calc=datum
      do 409 k=1,kmax
      denom=(x(i)-xs(k))**2+(y(i)-ys(k))**2+(zs(k)+ef*h(i))**2
      calc=calc+(c(k)/sqrt(denom))
409   continue
      dif=f(i)-calc
      if(abs(dif).gt.rmax)rmax=abs(dif)
      count=count+1
      sum=sum+dif
      ssq=ssq+dif*dif
410   continue
      ave=sum/count
      sig=(ssq/(count-1.))-ave
      write(6,411)rmax,ave,sig
411   format(' Max error = ',e12.4,', Mean = ',e12.4,', Sigma =',e12.4)
      close(10)
      close(11)
      stop
      end
   
