c ESMAG.FOR
c
c Finds equivalent dipole sources that fit magnetic data in grid file format.
c The equivalent dipoles are oriented in the direction of the regional field,
c and are located on a specified surface (ie. the topography).
c
c Jeff Phillips - Oct 1993
c restricted to 10000 sources and 3000 columns
      parameter (MAXP = 10000)
      parameter (MAXF = 3000)
      character*50 ifile,ofile
      character*1 ans,drive
      dimension id(2)
c      real*8 x(MAXP),y(MAXP),h(MAXP)
c      real*8 xp(MAXP),yp(MAXP),f(MAXP),c(MAXP),aj(MAXP)
c      real*8 a,xd,yd,fmean,fx2,fsum,cnst,dipole,datum
      dimension x(MAXF),y(MAXF),h(MAXF),f(MAXF),ah(MAXF)
      dimension xp(MAXP),yp(MAXP),c(MAXP),aj(MAXP),a(MAXP)
      dimension ic(MAXP),jc(MAXP)
      dimension title(14),pgm(2)
c      real*8 f(MAXP)
      real*8 fd(MAXP),ddval
      real mx,my,mz
      common /normal/mx,my,mz
      data pi/3.14159265/
      data dval/1.e38/,ddval/1.d38/
c  Read data.
799   print*,'Enter magnetic field grid file name:'
      print 1
1     format(' *',$)
      read(5,800)ifile
800   format(a)
      open(10,file=ifile,form='unformatted',status='old')
      read(10) title,pgm,nc,nr,nz,xo,dx,yo,dy
      print*,'Enter observation surface grid file name:'
      print 1
      read(5,800)ifile
      open(11,file=ifile,form='unformatted',status='old')
      read(11) title,pgm,nc1,nr1,nz,xo1,dx1,yo1,dy1
      print*,'Enter source surface grid file name:'
      print 1
      read(5,800)ifile
      open(12,file=ifile,form='unformatted',status='old')
      read(12) title,pgm,nc2,nr2,nz,xo2,dx2,yo2,dy2
      if(nc1.ne.nc.or.nc1.ne.nc2) stop 'grids not registered'
      if(nr1.ne.nr.or.nr1.ne.nr2) stop 'grids not registered'
      if(xo1.ne.xo.or.xo1.ne.xo2) stop 'grids not registered'
      if(dx1.ne.dx.or.dx1.ne.dx2) stop 'grids not registered'
      if(yo1.ne.yo.or.yo1.ne.yo2) stop 'grids not registered'
      if(dy1.ne.dy.or.dy1.ne.dy2) stop 'grids not registered'
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
      print'(a\)',' Enter the (RAM) drive letter for temporary files: '
      read(5,800) drive
      open(20,file=drive//':\f.tmp',form='unformatted',
     1 status='unknown',access='direct',recl=nc*8)
      open(21,file=drive//':\h.tmp',form='unformatted',
     1 status='unknown',access='direct',recl=nc*4)
      open(22,file=drive//':\hh.tmp',form='unformatted',
     1 status='unknown',access='direct',recl=nc*4)
c     hmin=1.e20
c     hmax=-1.e20
      fmean=0.
      nf=0
      do 15 i=1,nc
      x(i)=xo+(i-1)*dx
   15 continue
      icount=0
      do 20 j=1,nr
      y(j)=yo+(j-1)*dy
      read(10) dlt,(f(i),i=1,nc)
      read(11) dlt,(h(i),i=1,nc)
      read(12) dlt,(ah(i),i=1,nc)
      do 10 i=1,nc
      fd(i)=ddval
      if(h(i).ge.dval.or.ah(i).ge.dval) then
c       f(i)=dval
        go to 10
c     endif
      else
      if(f(i).ge.dval) then
        h(i)=dval
        ah(i)=dval
        go to 10
      endif
      endif
      fmean=fmean+f(i)
      nf=nf+1
c  Elevation converted to +down system, in km.
      fd(i)=f(i)
      h(i)=-h(i)*efact
      ah(i)=-ah(i)*efact
c      if(ah(i).le.h(i)) stop 'source above observation'
      if(ah(i).le.h(i)) then
        h(i)=ah(i)-.001
        icount=icount+1
      endif
c     if(h(i).gt.hmax) hmax=h(i)
c     if(h(i).lt.hmin) hmin=h(i)
   10 continue
      write(20,rec=j) (fd(i),i=1,nc)
      write(21,rec=j) (h(i),i=1,nc)
      write(22,rec=j) (ah(i),i=1,nc)
   20 continue
      close(10)
      close(11)
      close(12)
      if(icount.gt.0) then
        print*
        print*,'WARNING: the observation surface fell below the source'
        print*,'surface at ',icount,' points.  The observation surface'
        print*,'was raised slightly above the source surface at these '
        print*,'locations.'
      endif
      fmean=fmean/nf
c     write(6,42)hmin,hmax
c2    format(' hmin,max = ',2e12.3,' km')
c  Begin computations.
      print*
c     print*,' a = equivalent source depth'
      print*,' nmax = maximum iterations (eg. 200).'
      print*,' eps = desired maximum error (eg. 1.0)'
      print*
c     print*,' Enter a,nmax,eps:'
      print*,' Enter nmax,eps:'
      print 1
c     read*, a,nmax,eps
      read*, nmax,eps
      print*,' Enter inclination, declination:'
      print 1
      read*, ainc,decl
      mz=ainc*pi/180.
      my=decl*pi/180.
      mx=cos(mz)*cos(my)
      my=cos(mz)*sin(my)
      mz=sin(mz)
      print*,'Enter output sources file name:'
      print 1
      read(5,800)ofile
c     open(13,file=ofile,form='unformatted',status='unknown')
c      print*,'Enter smoothing factor (0<alpha<1):'
c      read(5,*) alpha
c  Iterations (kk denotes the source).
      print*,' iteration   sources   max_err   std_err'
      n=0
      ns=0
      afmin=1.e20
      afstd=1.e20
      datum=0.
c     smoth=0.
220   n=n+1
      af0=0
      datum=datum+fmean
      do 40 j=1,nr
      read(20,rec=j) (fd(i),i=1,nc)
      read(21,rec=j) (h(i),i=1,nc)
      read(22,rec=j) (ah(i),i=1,nc)
      do 30 i=1,nc
      if(fd(i).ge.ddval) go to 30
      fd(i)=fd(i)-fmean
      af=abs(fd(i))
      if(af.lt.af0) go to 30
      ii=i
      jj=j
      hh=h(i)
      ahh=ah(i)
      ff=fd(i)
      af0=af
   30 continue
      write(20,rec=j) (fd(i),i=1,nc)
   40 continue
      if(n.eq.1)afmax=af0
      if(af0.le.eps)go to 300
      if(ns.eq.0) go to 60
      do 50 k=1,ns
      if(ic(k).eq.ii.and.jc(k).eq.jj) go to 70
   50 continue
c  Find dipole location
   60 continue
        ns=ns+1
        if(ns.gt.MAXP) then
          print*,'exceeded source array size'
          go to 300
        endif
        call dipoloc(y(jj),x(ii),hh,yp(ns),xp(ns),ahh,aj(ns))
        a(ns)=ahh
        k=ns
        ic(k)=ii
        jc(k)=jj
   70 xd=xp(k)
      yd=yp(k)
      cnst=ff/aj(k)
c      cnst=aj(k)*ff/(alpha+aj(k)**2)
c  Remove this source.
      af0=0
      fmean=0
      fx2=0
c      do 211 k=1,kmax
      do 90 j=1,nr
      read(20,rec=j) (fd(i),i=1,nc)
      read(21,rec=j) (h(i),i=1,nc)
      do 80 i=1,nc
      if(fd(i).ge.ddval) go to 80
      fd(i)=fd(i)-dipole(yd,xd,ahh,cnst,y(j),x(i),h(i))
      af=abs(fd(i))
      if(af.lt.af0) go to 312
      af0=af
312   fx2=fx2+fd(i)*fd(i)
      fmean=fmean+fd(i)
   80 continue
      write(20,rec=j) (fd(i),i=1,nc)
   90 continue
      c(k)=c(k)+cnst
c      smoth=0.0
c      do 212 k=1,ns
c      smoth=smoth+abs(c(k))
c212   continue
c      smoth=smoth/float(ns)
      fmean=fmean/float(nf)
      afstd=sqrt((fx2-float(nf)*fmean*fmean)/(float(nf)-1.))
        if(af0.lt.afmin) then
c          write(6,412) n,ns,af0,afstd,fmean,smoth
          write(6,412) n,ns,af0,afstd
        endif
c       if(af0.lt.afmin) print *
c  Collapse source array and output sources file.
      if(af0.lt.afmin) then
        open(13,file=ofile,form='unformatted',status='unknown')
        call writeit(ns,c,xp,yp,id,a,n,datum,ainc,decl)
        close(13)
        afmin=af0
        itmin=n
      endif
412   format(1h ,2i10,4f10.3)
      if(n.lt.nmax)go to 220
c  Iterations completed.
300   if(n.lt.nmax)nmax=n
c      call endit()
      print*
      print*,'The smallest maximum error ',afmin,' occurred at'
      print*,'iteration ',itmin
c      print*
c      print*,'Choose different parameters? (y or n): '
c      read(5,800)ans
c      if(ans.eq.'y'.or.ans.eq.'Y') go to 799
      print*
      print*,ns,' sources written to file ',ofile
      print*
      stop
      end
      subroutine writeit(ksor,c,xp,yp,id,a,n,datum,ainc,decln)
      dimension c(ksor+1),xp(ksor+1),yp(ksor+1),a(ksor+1)
      dimension v(6),id(2)
c     rewind(13)
c     v(1)=a
      v(3)=n
      v(4)=datum
      v(5)=ainc
      decl=decln
      v(6)=decl
      do 10 j=1,ksor
      xs=xp(j)
      ys=yp(j)
      v(1)=a(j)
      v(2)=c(j)
      write(13)id(1),id(2),xs,ys,v
   10 continue
      return
      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 rplot(rmin, rmax, hmin, hmax, iplotr)
c      dimension xp(4), yp(4), dxp(2), dyp(2)
c      dimension x(n), y(n)
c      call pltset(iplotr, xp(4), yp(4), 1)
c      xp(4) = min(xp(4),10.)
c      yp(4) = min(yp(4),8.)
c      xp(4) = xp(4) - .25
c      yp(4) = yp(4) - .25
c      dxp(1) = rmin
c      dxp(2) = (rmax + .01) - mod(rmax,.01)
c      dyp(1) = hmin
c      dyp(2) = hmax
c      xp(1) = xp(4) - 1.5
c      xp(2) = 0.
c      xp(3) = 1.25
c      yp(1) = yp(4) - 1.25
c      yp(2) = 0.
c      yp(3) = 1.
c      call scale(dxp, dyp, xp, yp, 4, ier)
c      if (ier .ne. 0) goto 90
c   53 call yaxis(dyp, dxp, yp, (hmax-hmin)/20., 2, .12, '(f5.0)', 4)
c   56 call xaxis(dxp, dyp, xp, (rmax-rmin)/20., 2, .12, '(f5.0)', 4)
c   57 call neatl
c   58 call vchar((dxp(1) + dxp(2)) * .5, dyp(1), 'iterations', 10, 2,
c     &.12, 0., -.54, -.64)
c   61 call vchar(dxp, (dyp(1) + dyp(2)) * .5, 'max error', 9, 2, .12,
c     &1.5706, -.48, .79)
c   64 return
c   90 stop
c      entry curv(x, y, n, m)
c      call line(x, y, n, 0, m)
c      return
c      entry endit()
c      call endpt(ie)
c      return
c      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
c      real*8 function dipole(xd,yd,zd,cnst,x,y,z)
      function dipole(xd,yd,zd,cnst,x,y,z)
c     note: x is north, y is east, z is down
c      real*8 xd,yd,zd,x,y,z,rx,ry,rz,r2,r,cnst
      real mx,my,mz
      common /normal/mx,my,mz
      rx=(x-xd)
      ry=(y-yd)
      rz=(z-zd)
      r2=rx**2+ry**2+rz**2
      r=sqrt(r2)
      dipole=cnst*(3.*(mx*rx+my*ry+mz*rz)**2/r2-1.)/r**3
      end
      subroutine dipoloc(x,y,z,xd,yd,zd,cnst)
c     input x,y,z,zd
c     returns xd,yd,cnst
c     note: x is north, y is east, z is down
c      real*8 del,delx,dely,x,y,z,xs1,ys1,xs2,ys2,xs3,ys3,xd,yd,zd,q
c      real*8 cnst,cs1,cs2,cs3,dipole,sign
      real mx,my,mz
      common /normal/mx,my,mz
      del=abs(z-zd)/10.
      delx=mx*del
      dely=my*del
      sign=1.d0
c      print*,'source: ',x,y,z,' level: ',zd
      if(abs(mz).eq.1.0.or.mz.eq.0.0) then
        xd=x
        yd=y
        cnst=dipole(x,y,z,sign,xd,yd,zd)
        return
      else if(mz.lt.0.0) then
        sign=-1.d0
c        del=-del
c        delx=-delx
c        dely=-dely
      endif
      xs1=x-delx
      ys1=y-dely
      cs1=dipole(x,y,z,sign,xs1,ys1,zd)
c      print*,xs1,ys1,cs1
      xs2=x
      ys2=y
      cs2=dipole(x,y,z,sign,xs2,ys2,zd)
c      print*,xs2,ys2,cs2
  100 xs3=xs2+delx
      ys3=ys2+dely
      cs3=dipole(x,y,z,sign,xs3,ys3,zd)
c      print*,xs3,ys3,cs3
      call interp(cs1,cs2,cs3,del,q,cnst,ierr)
      if(ierr.eq.0) then
        cnst=cnst*sign
        if(mz.gt.0.0) then
          xd=xs2+mx*q
          yd=ys2+my*q
        else
          xd=xs2-mx*q
          yd=ys2-my*q
        endif
c        print*, x,y,xd,yd,cnst
c        print*
        return
      endif
      xs1=xs2
      ys1=ys2
      cs1=cs2
      xs2=xs3
      ys2=ys3
      cs2=cs3
      go to 100
      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
c      real*8 d,xmax,a,b,c,y1,y2,y3,ymax
      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

