#define _FILE_OFFSET_BITS 64
#define _LARGEFILE_SOURCE
c----------( EMAGF : preparation for <amag> inversion )----------
      parameter (LHW=100, LXO=1201, LYO=1201)
      parameter (LXT=LXO+LHW+LHW, LYT=LYO+LHW+LHW)
c
      dimension g(-LHW:LHW, -LHW:LHW)
      dimension t(LXT,LYT), f(LXO,LYO), h(LXO,LYO), e(6)
      character*80 wdr, flog, fnam1, fnam2, fnam3, fnam4
      character*8 area, tarea, area1, salt
      character buf*80, strt*25, strs*20, out*72
c
      call premsg('EMAGF : preparation for <amag> inversion')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnam1 = wdr
      fnam2 = wdr
      fnam3 = wdr
      fnam4 = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call gparma('Enter Mag.Anomaly data filename   ==> ',
     *            81-ldr, fnam1(ldr:80))
      open(12,file=fnam1,status='old')
    1 read(12,'(a)') buf
      if(buf(1:1).eq.'#') goto 1
      read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      nc = ncc
      if(ncc.ge.200) nc = ncc - 200
      if(nc.ge.1.and.nc.le.62) then
        iorg = 0
        ispa = 0
        ispb = 0
        korg = (nc-30)*360 - 180
        if(nc.gt.60) korg = 0
        if(nc.eq.61) iorg = +5400
        if(nc.eq.62) iorg = -5400
      endif
      read(12,'(a)') buf
      read(buf,*) ixs, iys, mszx, mszy, mx, my, vnul, alt
      write(salt,'(1x,f7.0)') alt
      if (alt.eq.0.) salt = '    var.'
      if (alt.lt.0.) salt = '   undef'
      xs = float(ixs)/1000.
      ys = float(iys)/1000.
      if(mx.gt.LXO.or.my.gt.LYO) call abendm('too large array size')
      read(12,*) ((f(i,j),i=1,mx),j=1,my)
c
      call gparma('Enter Src.Altitude data filename  ==> ',
     *            81-ldr, fnam2(ldr:80))
      open(10,file=fnam2,status='old')
    2 read(10,'(a)') buf
      if(buf(1:1).eq.'#') goto 2
      read(buf,'(a8,i4,4x,4i8)') tarea, ncct, iorgt,korgt, ispat,ispbt
      nct = ncct
      if(ncct.ge.200) nct = ncct - 200
      if(nct.ge.1.and.nct.le.62) then
        iorgt = 0
        ispat = 0
        ispbt = 0
        korgt = (nct-30)*360 - 180
        if(nct.gt.60) korgt = 0
        if(nct.eq.61) iorgt = +5400
        if(nct.eq.62) iorgt = -5400
      endif
      read(10,*) ixst, iyst, mszxt, mszyt, mxt, myt, tnul
      xst = float(ixst)/1000.
      yst = float(iyst)/1000.
      if(mxt.gt.LXT.or.myt.gt.LYT) call abendm('too large array size')
      read(10,*) ((t(i,j),i=1,mxt),j=1,myt)
      close(10)
c
      write(out,'(a,a8,a,i4,4x,4i8)')
     *     '--- Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
      call premsg(out)
      if(ncct.ne.ncc.or.iorgt.ne.iorg.or.korgt.ne.korg
     *            .or.ispat.ne.ispa.or.ispbt.ne.ispb) then
        write(out,'(11x,a8,6x,i4,4x,4i8)')
     *       tarea, ncct, iorgt,korgt, ispat,ispbt
        call premsg(out)
        call abendm('Header inconsistent')
      endif
      if(tarea.ne.area) call premsg('           ' // tarea)
      write(out,'(a,2a10,2x,2a8,2a6,a7)')
     *     '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my',
     *     'alt'
      call premsg(out)
      write(out,'(a,2f10.2,2i8,2i6,a8)')
     *     '--- mag. obs.   ', xs, ys, mszx, mszy, mx, my, salt
      call premsg(out)
      write(out,'(a,2f10.2,2i8,2i6)')
     *     '--- source alt. ', xst, yst, mszxt, mszyt, mxt, myt
      call premsg(out)
      xt = xs + float(mszx*(mx-1))/1000.
      yt = ys + float(mszy*(my-1))/1000.
      fmx = float(mszx)
      fmy = float(mszy)
      xtt = xst + float(mszxt*(mxt-1))/1000.
      ytt = yst + float(mszyt*(myt-1))/1000.
      fmxt = float(mszxt)
      fmyt = float(mszyt)
      if(xst.gt.xs.or.yst.gt.ys .or.
     *   xtt.lt.xt.or.ytt.lt.yt) call abendm('ranges conflict')
c
      if(alt.eq.0.) goto 3
      if(alt.lt.0.) then
        call gparmf('    Mag. Obs. Altitude [in m]  ==> ', alt)
        if (alt.le.0.) call abendm('invalid value')
        write(salt,'(1x,f7.0)') alt
      endif
      write(strt,'(f7.0,a)') alt, ' m above sea level'
      goto 4
    3 read(12,'(a)') buf
      if(buf(1:1).eq.'#') goto 3
      read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      nc1 = ncc1
      if(ncc1.ge.200) nc1 = ncc1 - 200
      if(nc1.ge.1.and.nc1.le.62) then
        iorg1 = 0
        ispa1 = 0
        ispb1 = 0
        korg1 = (nc1-30)*360 - 180
        if(nc1.gt.60) korg1 = 0
        if(nc1.eq.61) iorg1 = +5400
        if(nc1.eq.62) iorg1 = -5400
      endif
      if(area1.ne.area.or.ncc1.ne.ncc .or.
     *   iorg1.ne.iorg.or.korg1.ne.korg .or.
     *   ispa1.ne.ispa.or.ispb1.ne.ispb)
     *     call abendm('Header1 inconsistent')
      read(12,*) ixs1, iys1, mszx1, mszy1, mx1, my1, hnul
      if(ixs1.ne.ixs.or.iys1.ne.iys .or.
     *   mszx1.ne.mszx.or.mszy1.ne.mszy .or.
     *   mx1.ne.mx.or.my1.ne.my)
     *     call abendm('Header2 inconsistent')
      read(12,*) ((h(i,j),i=1,mx),j=1,my)
      strt = ' given as 2nd set data   '
    4 close(12)
      call premsg('  obs.surface : ' // strt)
c
      call premsg('')
      call premsg('Specify Source distribution Parameters')
      call premsg('  Southwest corner Coord. in km ?')
      call gparmf('    Northing : ', xl)
      ixl = ixst + mszxt*nint((xl-xst)*1000./fmxt)
      xl = float(ixl)/1000.
      if(xl.lt.xst.or.xl.gt.xs) call abendm('ranges conflict')
      call gparmf('    Easting  : ', yl)
      iyl = iyst + mszyt*nint((yl-yst)*1000./fmyt)
      yl = float(iyl)/1000.
      if(yl.lt.yst.or.yl.gt.ys) call abendm('ranges conflict')
      call premsg('  Mesh size and numbers to N ?')
      call gparmi('    Unit Mesh size in m :  ', mszxm)
      if(mszxm.le.0) call abendm('invalid value')
      nmsx = mszxm / mszxt
      if((nmsx*mszxt).ne.mszxm) call abendm('mesh size illegal')
      hnmsx = float(nmsx-1)/2.
      fmxm = float(mszxm)
      call gparmi('    Number of Mesh to N :  ', mxm)
      xh = xl + float(mszxm*(mxm-1))/1000.
      if(xh.lt.xt.or.xh.gt.xtt) call abendm('ranges conflict')
      call premsg('  Mesh size and numbers to E ?')
      call gparmi('    Unit Mesh size in m :  ', mszym)
      if(mszym.le.0) call abendm('invalid value')
      nmsy = mszym / mszyt
      if((nmsy*mszyt).ne.mszym) call abendm('mesh size illegal')
      hnmsy = float(nmsy-1)/2.
      fmym = float(mszym)
      call gparmi('    Number of Mesh to E :  ', mym)
      yh = yl + float(mszym*(mym-1))/1000.
      if(yh.lt.yt.or.yh.gt.ytt) call abendm('ranges conflict')
      call premsg('  Truncation of Source effect ?')
      call gparmf('    Half width (in km) :  ', hw)
      nhx = nint(hw*1000./fmx)
      nhy = nint(hw*1000./fmy)
      if(nhx.gt.LHW.or.nhy.gt.LHW)
     *     call abendm('shortage of work area within truncation')
      if(nhx.lt.10.or.nhy.lt.10) call abendm('truncation too short')
c
c      frnfxy = float((nhx+1+nhx)*(nhy+1+nhy) + 4) * 4. / 1024.
c      frmxym = float(mxm*mym + 1) / 1024. / 1024.
c      if((frnfxy*frmxym).gt.2.) then
c        call abendm('COEF file-size exceeds 2GB.')
c      endif
c
      call premsg('')
      call premsg('Specify Source-Bottom configuration')
      call prompt('    Constant Thickness (1) or Flat Bottom (0) ')
      call gparmi(': ', ks)
      if(ks.lt.0.or.ks.gt.1) call abendm('invalid value')
      strs = ' below sea level    '
      if(ks.eq.1) strs = ' below terr.surface '
      call gparmf('  Bottom-depth (in m)' // strs // ': ', botm)
c
      call premsg('Specify Direction Parameters (in degrees)')
      call gparmf2('  Inclin., Declin. of Ambient field : ', dip, dec)
      if(dip.lt.-90..or.dip.gt.90. .or.
     *   dec.lt.-180..or.dec.ge.360.) call abendm('invalid value')
      call gparmf2('  Inclin., Declin. of Magnetization : ', the, phi)
      if(the.lt.-90..or.the.gt.90. .or.
     *   phi.lt.-180..or.phi.ge.360.) call abendm('invalid value')
c
      call premsg('')
      call gparma('Enter COEF matrix output filename ==> ',
     *            81-ldr, fnam3(ldr:80))
      open(20,file=fnam3,form='unformatted',status='new')
      call gparma('Enter AMAG model output filename  ==> ',
     *            81-ldr, fnam4(ldr:80))
      call clspin()
      open(21,file=fnam4,status='new')
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'EMAGF : preparation for <amag> inversion'
      write(6,'(a,a)') ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)') ' Src. Alt. infile    : ', fnam2(1:lrtrim(fnam2))
      write(6,'(a,a)') ' COEF matrix outfile : ', fnam3(1:lrtrim(fnam3))
      write(6,'(a,a)') ' AMAG model outfile  : ', fnam4(1:lrtrim(fnam4))
      write(6,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xs,  ys,  mszx,  mszy,  mx,  my, salt,
     *              xst, yst, mszxt, mszyt, mxt, myt,
     *              xl,  yl,  mszxm, mszym, mxm, mym
      write(6,2000) hw, nhx, nhy
      write(6,3000) strt, botm, strs, dip, dec, the, phi
c
      call prompt('... gen. COEF matrix data : ')
      call dpcent(0, 0)
      xlm = (xl-xst) * 1000.
      ylm = (yl-yst) * 1000.
      xsm = (xs-xst) * 1000.
      ysm = (ys-yst) * 1000.
      rad = 3.141593 / 180.
      ex = cos(dip*rad) * cos(dec*rad)
      ey = cos(dip*rad) * sin(dec*rad)
      ez = sin(dip*rad)
      px = cos(the*rad) * cos(phi*rad)
      py = cos(the*rad) * sin(phi*rad)
      pz = sin(the*rad)
      cx = ex * px
      cy = ey * py
      cz = ez * pz
      cxy = (ex*py + ey*px) / 2.
      cyz = (ey*pz + ez*py) / 2.
      czx = (ez*px + ex*pz) / 2.
      c3 = cx - cz
      c4 = cy - cz
      c5 = cx - cy
      c1 = c3 + c5
      c2 = c4 - c5
      write(20) area, ncc, iorg,korg, ispa,ispb,
     *          ixs, iys, mszx,  mszy,  mx,  my,  alt,
     *          ixl, iyl, mszxm, mszym, mxm, mym, nhx, nhy
      am = fmxt*fmyt
c
      do 200 j=1,mym
        yc = ylm + fmym*float(j-1)
        ja = nint(yc/fmyt - hnmsy)
        kjc = nint((yc-ysm)/fmy)
        do 200 i=1,mxm
          xc = xlm + fmxm*float(i-1)
          ia = nint(xc/fmxt - hnmsx)
          kic = nint((xc-xsm)/fmx)
          do 300 kj=-LHW,LHW
          do 300 ki=-LHW,LHW
            g(ki,kj) = 0.
  300     continue
          do 400 js=1,nmsy
            jt = ja + js
            if(jt.le.0.or.jt.gt.myt) goto 400
            ym = float(jt-1)*fmyt - ysm
            do 500 is=1,nmsx
              it = ia + is
              if(it.le.0.or.it.gt.mxt) goto 500
              xm = float(it-1)*fmxt - xsm
              if(t(it,jt).eq.tnul) goto 500
              do 600 kj=-nhy,nhy
               jr = kjc + kj
               y = ym - float(jr)*fmy
               do 600 ki=-nhx,nhx
                ir = kic + ki
                x = xm - float(ir)*fmx
                if(ir.lt.0.or.ir.ge.mx.or.jr.lt.0.or.jr.ge.my)
     *              goto 600
                if(f(ir+1,jr+1).eq.vnul) goto 600
                hc = alt
                if(alt.eq.0.) hc = h(ir+1,jr+1)
                za = hc - t(it,jt)
                zb = hc + botm
                if(ks.eq.1) zb = za + botm
                if(abs(za).lt.1.) za = 1.
                if(abs(zb).lt.1.) zb = 1.
                if(iabs(ki).le.10.and.iabs(kj).le.10) then
                  call prism(x, y, fmxt, fmyt, za, zb, e)
                  v = cyz*e(1) + czx*e(2) + cxy*e(3)
     *                + cx*e(4) + cy*e(5) + cz*e(6)
                else
                  zza = za * za
                  zzb = zb * zb
                  xx = x * x
                  xy = x * y
                  yy = y * y
                  ss = xx + yy
                  ssss = ss * ss
                  rra = ss + zza
                  rrb = ss + zzb
                  ra = sqrt(rra)
                  rb = sqrt(rrb)
                  dza = za / ra
                  dzb = zb / rb
                  u1a = dza / ssss
                  u1b = dzb / ssss
                  u2a = u1a * zza / rra
                  u2b = u1b * zzb / rrb
                  u3a = 2. / ra / rra
                  u3b = 2. / rb / rrb
                  v1 = c1*xx + c2*yy + cxy*xy*6.
                  v2 = c3*xx + c4*yy + cxy*xy*2.
                  v3 = czx*x + cyz*y
                  v = ((u1b-u1a)*v1 - (u2b-u2a)*v2 - (u3b-u3a)*v3) * am
                endif
                g(ki,kj) = g(ki,kj) + v
  600         continue
  500       continue
  400     continue
          write(20) i,j, kic,kjc, ((g(ki,kj),ki=-nhx,nhx),kj=-nhy,nhy)
          call dpcent((j-1)*mxm+i, mxm*mym)
  200 continue
      close(20)
c
      ixl = nint(xl*1000.)
      iyl = nint(yl*1000.)
      write(21,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      write(21,'(2i12,4i6,1x,f7.1,1x,f7.0,i6)')
     *     ixl, iyl, mszxm, mszym, mxm, mym, 9999.9, -1., 0
      do 40 j=1,mym
        write(21,'((f7.1,9(1x,f7.1)))') (0.,i=1,mxm)
   40 continue
      close(21)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,status='old',iostat=ios)
        if(ios.eq.0) then
   90     read(99,*,end=91)
          goto 90
   91     backspace(99)
        else
          open(99,file=flog,status='new')
        endif
        write(99,'(//a)') 'EMAGF : preparation for <amag> inversion'
        write(99,'(a,a)')
     *       ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)')
     *       ' Src. Alt. infile    : ', fnam2(1:lrtrim(fnam2))
        write(99,'(a,a)')
     *       ' COEF matrix outfile : ', fnam3(1:lrtrim(fnam3))
        write(99,'(a,a)')
     *       ' AMAG model outfile  : ', fnam4(1:lrtrim(fnam4))
        write(99,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xs,  ys,  mszx,  mszy,  mx,  my, salt,
     *                 xst, yst, mszxt, mszyt, mxt, myt,
     *                 xl,  yl,  mszxm, mszym, mxm, mym
        write(99,2000) hw, nhx, nhy
        write(99,3000) strt, botm, strs, dip, dec, the, phi
        close(99)
      endif
      stop
c
 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/
     *       3x, '   (mesh)   ', '      xs  ', '      ys  ',
     *       '   msz-x', '   msz-y', '    mx', '    my', '    alt '/
     *       3x, 'mag. obs.   ', 2f10.2, 2i8, 2i6, a8/
     *       3x, 'source alt. ', 2f10.2, 2i8, 2i6/
     *       3x, 'source model', 2f10.2, 2i8, 2i6)
 2000 format(' Half width of calc.window :', f6.2, ' km (', i4,
     *       ' *', i4, '  mag.mesh )')
 3000 format('  obs. surface  : ', a25/
     *       '  bottom depth  : ', f8.1, ' m', a20/
     *       '  ambient field :  dip=', f6.1, 3x, 'dec=', f7.1/
     *       '  magnetization :  dip=', f6.1, 3x, 'dec=', f7.1)
      end
c---
      subroutine prism(xc, yc, fmx, fmy, za, zb, e)
      dimension e(6)
      do 1 i=1,6
        e(i) = 0.
    1 continue
      xa = xc - fmx/2.
      xb = xc + fmx/2.
      ya = yc - fmy/2.
      yb = yc + fmy/2.
      do 2 l=0,7
        x = xb
        y = yb
        z = zb
        if(mod(l,2).ne.0) x = xa
        if(mod(l/2,2).ne.0) y = ya
        if(l.ge.4) z = za
        rr = x*x + y*y + z*z
        r = sqrt(rr)
        e(1) = -e(1) + alog((r-x)/(r+x))
        e(2) = -e(2) + alog((r-y)/(r+y))
        e(3) = -e(3) - alog(r+z)
        e(4) = -e(4) - atan2((x*y),(rr+r*z-y*y))
        e(5) = -e(5) - atan2((x*y),(rr+r*z-x*x))
        e(6) = -e(6) + atan((x*y)/(r*z))
        if(l.ne.1.and.l.ne.5) goto 2
        do 3 i=1,6
          e(i) = -e(i)
    3   continue
    2 continue
      return
      end
