c----------( AAPTDP : point-dipole model fitting )----------
      parameter (LXF=1201, LYF=1201, MXPT=2601)
c
      parameter (RAD=180./3.14159265)
      dimension v(LXF,LYF), t(LXF,LYF), w(LXF,LYF)
      dimension c(3), sol(3), sol1(3), em(4,MXPT), f(3)
      character*80 wdr, flog, fnam1, fnam3, fnam4
      character*8 area, area1, rarea, mname, rname
      character buf*80, strt*25, ps*1, salt*8, out*72
c
      call premsg('AAPTDP : point-dipole model fitting')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnam1 = 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(1,file=fnam1,status='old')
    1 read(1,'(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(1,'(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.
      ixt = ixs + mszx*(mx-1)
      iyt = iys + mszy*(my-1)
      fmx = float(mszx)/1000.
      fmy = float(mszy)/1000.
      if(mx.gt.LXF.or.my.gt.LYF) call abendm('too large array size')
      read(1,*) ((v(i,j),i=1,mx),j=1,my)
c
      write(out,'(a,a8,a,i4,4x,4i8)')
     *     '--- Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
      call premsg(out)
      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)
c
      if(alt.eq.0.) goto 2
      if(alt.lt.0.) then
        call gparmf('    Altitude [in m] for this data ==> ', alt)
        if (alt.le.0.) call abendm('invalid Alt. value')
        write(salt,'(1x,f7.0)') alt
      endif
      write(strt,'(f7.0,a)') alt, ' m above sea level'
      tnul = 9999.9
      do 100 j=1,my
      do 100 i=1,mx
        t(i,j) = alt
  100 continue
      goto 3
    2 read(1,'(a)') buf
      if(buf(1:1).eq.'#') goto 2
      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('Header inconsistent')
      read(1,*) ixs1, iys1, mszx1, mszy1, mx1, my1, tnul
      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('Header inconsistent')
      read(1,*) ((t(i,j),i=1,mx),j=1,my)
      strt = ' given as 2nd set data   '
    3 close(1)
      call premsg('  obs.surface : ' // strt)
c
      do 200 j=1,my
      do 200 i=1,mx
        if(t(i,j).eq.tnul) then
          if(v(i,j).ne.vnul) call abendm('input file broken')
          w(i,j) = vnul
        else
          w(i,j) = 0.
        endif
  200 continue
c
      call premsg('')
      call premsg('Specify Direction Parameters (in degrees)')
   12 call gparmf2('  Inclin., Declin. of Ambient field : ', dip, dec)
      if(dip.lt.-90..or.dip.gt.90.) goto 12
      if(dec.lt.-180..or.dec.ge.360.) goto 12
      c(1) = cos(dip/RAD)*cos(dec/RAD)
      c(2) = cos(dip/RAD)*sin(dec/RAD)
      c(3) = sin(dip/RAD)
c
      call premsg('')
      call gparma('Enter Name-label for Model anom.  ==> ', 8, mname)
      call gparma('Enter Model anom. output filename ==> ',
     *            81-ldr, fnam3(ldr:80))
      open(11,file=fnam3,status='new')
      write(11,'(a)') '# AAPTDP model anomaly'
      write(11,'(a,a)') '# Source: ', fnam1(1:lrtrim(fnam1))
      write(11,'(a2,a8,i4,4x,4i8)')
     *     '# ', area, ncc, iorg, korg, ispa, ispb
      call premsg('')
      call gparma('Enter Name-label for Residual     ==> ', 8, rname)
      call gparma('Enter Residual output filename    ==> ',
     *            81-ldr, fnam4(ldr:80))
      open(12,file=fnam4,status='new')
      write(12,'(a)') '# AAPTDP residuals'
      write(12,'(a,a)') '# Source: ', fnam1(1:lrtrim(fnam1))
      write(12,'(a2,a8,i4,4x,4i8)')
     *     '# ', area, ncc, iorg, korg, ispa, ispb
      write(12,'(a,a)') '# Model : ', fnam3(1:lrtrim(fnam3))
      write(12,'(a2,a8,i4)') '# ', mname, ncc
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'AAPTDP : point-dipole model fitting'
      write(6,'(a,a)') ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
      write(6,'(a,a)') ' Model Anom. outfile : ', fnam3(1:lrtrim(fnam3))
      write(6,'(a,a)') ' Residual outfile    : ', fnam4(1:lrtrim(fnam4))
      write(6,1000) area, ncc, iorg, korg, ispa, ispb,
     *              xs, ys, mszx, mszy, mx, my, salt
      write(6,3000) strt, dip, dec
c
      open(9,status='scratch')
      call premsg(' ')
      call premsg('Which way to define window,')
      call premsg('    UTM coordinates in km ("u") or')
      call premsg('    mesh-count within the area ("m"),')
   13 call gparma('  Select  "u" or "m" : ', 1, ps)
      if(ps.ne.'u'.and.ps.ne.'m') goto 13
      goto 4
   43 call premsg('  too wide window')
      goto 4
   42 call premsg('  too narrow window')
      goto 4
   41 call premsg('  out of area')
    4 if(ps.eq.'u') then
        call premsg ('Select fitting window (in km)')
        call gparmf2('    SW-corner coord. (Xsw, Ysw)'
     *               // ' [both 0 to End] ==> ', xsw, ysw)
        if(xsw.eq.0..and.ysw.eq.0.) goto 8
        is = nint((xsw-xs)/fmx)
        js = nint((ysw-ys)/fmy)
        call gparmf2('    window size (N-S, E-W) ==> ', xsz, ysz)
        ni = nint(xsz/fmx)
        nj = nint(ysz/fmy)
      else
        call premsg ('Select fitting window'
     *               // ' (in mesh-unit starting from 0)')
        call gparmi2('    window size (N-S, E-W)'
     *               // ' [both 0 to End] ==> ', ni, nj)
        if(ni.eq.0.and.nj.eq.0) goto 8
        call gparmi2('    SW-corner location   (IXs, JYw)   ==> ',
     *               is, js)
      endif
      xsz = fmx*float(ni)
      ysz = fmy*float(nj)
      ni = ni + 1
      nj = nj + 1
      if(is.lt.0.or.js.lt.0) goto 41
      if((is+ni).gt.mx.or.(js+nj).gt.my) goto 41
      if(ni.le.5.or.nj.le.5) goto 42
      if((ni*nj).gt.MXPT) goto 43
      xis = fmx*float(is)
      yjs = fmy*float(js)
      xsw = xs + xis
      ysw = ys + yjs
c
      xsol = xsz/2.
      ysol = ysz/2.
      zsol = xsz/4.
      n = 0
      av = 0.
      do 300 j=1,nj
        y = ysol - fmy*float(j-1)
        do 300 i=1,ni
          if(v(is+i,js+j).eq.vnul) goto 300
          n = n + 1
          x = xsol - fmx*float(i-1)
          z = zsol + t(is+i,js+j)/1000.
          call adipol(x, y, z, c, em(1,n))
          em(4,n) = v(is+i,js+j)
          av = av + em(4,n)
  300 continue
      if(n.le.12) call abendm('too few available data')
      av = av / float(n)
      do 310 k=1,n
        em(4,k) = em(4,k) - av
  310 continue
      call lsqsol(n, em, sol, drms)
      step = 0.1
      ky = 0
      kyr = 0
    5 xsol1 = xsol
      ysol1 = ysol
      zsol1 = zsol
      if(ky.eq.0) then
        zsol1 = zsol + step
      else if(ky.eq.1) then
        zsol1 = zsol - step
      else if(ky.eq.2) then
        xsol1 = xsol + step
      else if(ky.eq.3) then
        xsol1 = xsol - step
      else if(ky.eq.4) then
        ysol1 = ysol + step
      else if(ky.eq.5) then
        ysol1 = ysol - step
      endif
      n = 0
      av = 0.
      do 400 j=1,nj
        y = ysol1 - fmy*float(j-1)
        do 400 i=1,ni
          if(v(is+i,js+j).eq.vnul) goto 400
          n = n + 1
          x = xsol1 - fmx*float(i-1)
          z = zsol1 + t(is+i,js+j)/1000.
          call adipol(x, y, z, c, em(1,n))
          em(4,n) = v(is+i,js+j)
          av = av + em(4,n)
  400 continue
      if(n.le.12) call abendm('too few available data')
      av = av / float(n)
      do 410 k=1,n
        em(4,k) = em(4,k) - av
  410 continue
      call lsqsol(n, em, sol1, drms1)
      if(drms1.ge.drms) goto 6
      kyr = ky
      drms = drms1
      xsol = xsol1
      ysol = ysol1
      zsol = zsol1
      sol(1) = sol1(1)
      sol(2) = sol1(2)
      sol(3) = sol1(3)
      goto 5
    6 ky = ky + 1
      if(ky.eq.6) ky = 0
      if(ky.ne.kyr) goto 5
      step = step / 2.
      ky = 0
      kyr = 0
      if(step.gt.0.001) goto 5
c
      ah = sqrt(sol(1)*sol(1) + sol(2)*sol(2))
      at = sqrt(sol(1)*sol(1) + sol(2)*sol(2) + sol(3)*sol(3))
      the = atan2(sol(3),ah)*RAD
      phi = atan2(sol(2),sol(1))*RAD
      write(6,4000) xsw, ysw, xsz, ysz, xsol, ysol, zsol,
     *              (sol(k),k=1,3), at, the, phi
      write(11,4000) xsw, ysw, xsz, ysz, xsol, ysol, zsol,
     *               (sol(k),k=1,3), at, the, phi
      write(9,*) xsw, ysw, xsz, ysz, xsol, ysol, zsol,
     *           (sol(k),k=1,3), at, the, phi
      do 500 j=1,my
        y = ysol + yjs - fmy*float(j-1)
        do 500 i=1,mx
          if(t(i,j).eq.tnul) goto 500
          x = xsol + xis - fmx*float(i-1)
          z = zsol + t(i,j)/1000.
          call adipol(x, y, z, c, f)
          ww = sol(1)*f(1) + sol(2)*f(2) + sol(3)*f(3)
          w(i,j) = w(i,j) + ww
  500 continue
      goto 4
c
    8 call clspin()
      do 600 j=1,my
      do 600 i=1,mx
        if(v(i,j).eq.vnul) goto 600
        v(i,j) = v(i,j) - w(i,j)
        if(w(i,j).eq.vnul) v(i,j) = vnul
  600 continue
c
      write(11,'(a8,i4,4x,4i8)') mname, ncc, iorg, korg, ispa, ispb
      write(11,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixs,iys, mszx,mszy, mx,my, vnul, alt
      do 700 j=1,my
        write(11,'((f7.1,9(1x,f7.1)))') (w(i,j),i=1,mx)
  700 continue
      if(alt.eq.0.) then
        write(11,'(a)') '# Var.Obs.Alt. data are copied from Source'
        write(11,'(a8,i4,4x,4i8)') mname, ncc, iorg, korg, ispa, ispb
        write(11,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *       ixs,iys, mszx,mszy, mx,my, tnul, -1.
        do 710 j=1,my
          write(11,'((f7.1,9(1x,f7.1)))') (t(i,j),i=1,mx)
  710   continue
      endif
      close(11)
      endfile(9)
      rewind(9)
c
      write(12,'(a8,i4,4x,4i8)') rname, ncc, iorg, korg, ispa, ispb
      write(12,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixs,iys, mszx,mszy, mx,my, vnul, alt
      do 800 j=1,my
        write(12,'((f7.1,9(1x,f7.1)))') (v(i,j),i=1,mx)
  800 continue
      if(alt.eq.0.) then
        write(12,'(a)') '# Var.Obs.Alt. data are copied from Source'
        write(12,'(a8,i4,4x,4i8)') rname, ncc, iorg, korg, ispa, ispb
        write(12,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *       ixs,iys, mszx,mszy, mx,my, tnul, -1.
        do 810 j=1,my
          write(12,'((f7.1,9(1x,f7.1)))') (t(i,j),i=1,mx)
  810   continue
      endif
      close(12)
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)') 'AAPTDP : point-dipole model fitting'
        write(99,'(a,a)')
     *     ' Mag. Anom. infile   : ', fnam1(1:lrtrim(fnam1))
        write(99,'(a,a)')
     *     ' Model Anom. outfile : ', fnam3(1:lrtrim(fnam3))
        write(99,'(a,a)')
     *     ' Residual outfile    : ', fnam4(1:lrtrim(fnam4))
        write(99,1000) area, ncc, iorg, korg, ispa, ispb,
     *                 xs, ys, mszx, mszy, mx, my, salt
        write(99,3000) strt, dip, dec
   92   read(9,*,end=93) xsw, ysw, xsz, ysz, xsol, ysol, zsol,
     *                   (sol(k),k=1,3), at, the, phi
        write(99,4000) xsw, ysw, xsz, ysz, xsol, ysol, zsol,
     *                 (sol(k),k=1,3), at, the, phi
        goto 92
   93   close(99)
      endif
      close(9)
      stop
c
 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/
     *       3x, '   (mesh)   ', '      xs  ', '      ys  ',
     *       '   msz-x', '   msz-y', '    mx', '    my', '    alt '/
     *       3x, 'obs. mag.   ', 2f10.2, 2i8, 2i6, a8)
 3000 format('  obs. surface  : ', a25/
     *       '  ambient field :  dip=', f6.1, 3x, 'dec=', f7.1)
 4000 format('#W# ', 2f9.3, 2f8.3, '  : ', 3f8.3/
     *       '##', 4x, 3f10.3, 2x, f10.3, 2f8.2)
      end
c---
      subroutine lsqsol(n, em, sol, drms)
      dimension em(4,*), sol(3)
      dimension q(4,3)
      do 1 i=1,3
      do 1 ii=1,4
        q(ii,i) = 0.
    1 continue
      do 2 k=1,n
        do 2 i=1,3
          do 2 ii=i,4
            q(ii,i) = q(ii,i) + em(i,k)*em(ii,k)
    2 continue
      q(1,2) = q(2,1)
      q(1,3) = q(3,1)
      q(2,3) = q(3,2)
      do 3 k=1,3
      do 3 ii=k+1,4
        q(ii,k) = q(ii,k) / q(k,k)
        do 3 i=1,3
          if(i.ne.k) q(ii,i) = q(ii,i) - q(ii,k)*q(k,i)
    3 continue
      sol(1) = q(4,1)
      sol(2) = q(4,2)
      sol(3) = q(4,3)
      sum = 0.
      do 4 k=1,n
        g = em(1,k)*q(4,1) + em(2,k)*q(4,2) + em(3,k)*q(4,3)
        sum = sum + (em(4,k)-g)*(em(4,k)-g)
    4 continue
      drms = sqrt(sum/float(n))
      return
      end
c---
      subroutine adipol(x, y, z, c, f)
      dimension c(3), f(3)
      xx = x * x
      yy = y * y
      zz = z * z
      xy = x * y
      yz = y * z
      zx = z * x
      rr = xx + yy + zz
      r = sqrt(rr)
      r5 = rr * rr * r / 100.
      f(1) = ((xx*3.-rr)*c(1) + xy*c(2)*3. + zx*c(3)*3.) / r5
      f(2) = ((yy*3.-rr)*c(2) + yz*c(3)*3. + xy*c(1)*3.) / r5
      f(3) = ((zz*3.-rr)*c(3) + zx*c(1)*3. + yz*c(2)*3.) / r5
      return
      end
