      program field
c
c---------------------------------------------------------------------
c
c     example program for evaluating time-dependent field model
c
c     Jeremy Bloxham & Andrew Jackson
c
c
c---------------------------------------------------------------------
c
c     uses code by: David Gubbins, Kathy Whaler, David Barraclough,
c                   Rick O'Connell, and Carl de Boor
c
c---------------------------------------------------------------------
c     for details of B-splines see:
c            Carl de Boor "A Practical Guide to Splines"
c            Springer-Verlag, 1978.
c---------------------------------------------------------------------
c
c     lmax   is maximum degree of spherical harmonics
c     nspl   is number of B-splines
c
c     gt     is the full array of coefficients in the expansion of
c               the geomagnetic potential in spherical harmonics
c               and cubic B-splines
c
c     g      is an array of geomagnetic main field coefficients at a
c               particular point in time
c
c     gd     is an array of geomagnetic secular variation coefficients
c               at a particular point in time
c
c     tknts  is an array of knot points
c
c     p      is an array of real associated Legendre polynomials
c
c     dp     is an array of derivatives of associated Legendre
c               polynomials with respect to colatitude
c---------------------------------------------------------------------
c
c     inorm = 1 (selects Gauss-Schmidt normalization)
c     jord  = 4 (order of B-splines)
c
c---------------------------------------------------------------------
c
c     CALLS:    interv   - calculates which knot lies immediately left
c                          of the current time point
c
c               bspline  - calculates B-splines at current time point
c
c               bspline1 - calculates first time derivative of
c                          B-splines at current time point
c
c               b0, b1   - functions required by bspline1
c
c               plmbar   - calculates associated Legendre polynomials
c                          and derivatives
c
c               magfdz   - evaluates field model
c
c               icoord   - performs transformation between geodetic and
c                          geocentric coords
c----------------------------------------------------------------------

      parameter (lmax=14)
      parameter (nspl=163)

      parameter (n=lmax*(lmax+2))
      parameter (np=n*nspl)
      parameter (nl=(lmax+1)*(lmax+2)/2)

      dimension gt(n,nspl)
      dimension spl(nspl),tknts(nspl+4)

      dimension g(n),gd(n)

      dimension p(nl),dp(nl)
      dimension dx(lmax*(lmax+2)),dy(lmax*(lmax+2)),dz(lmax*(lmax+2))

      data inorm/1/
      data jord/4/
      data fac/1.74532925e-2/

c*********************************************************************
c
c     the next line will need to be changed by the user
c     to suit the particular user's environment
c

      open(1,file='gufm1')
      open(7,file='dec.dat')
      open(8,file='inc.dat')

c*********************************************************************
c
c     in this example we compute the field at a point on the Earth's
c     surface close to Harvard Yard in 1985.0
c     alat, alon are the latitude and longitude in degrees
c     alt is the altitude in km above the earth's surface,
c     taken as radius 6371.2km
c     icoord choses geodetic (icoord=1) or geocentric (icoord=0) coordinates
c     coordinates are geodetic in this example

c      alat =   59+55./60.
c      alon =    10+45./60.
c      alt  =    0.0

      print*,' Program writes declination and inclination in degrees '
      print*,' for a given site 1590-1990 to files dec.dat and inc.dat '

      print*,'  '
      print*,' Use decimal degrees [1] or degrees and minutes [2]?'
      read(5,*)iresp
      if(iresp.eq.1)then

      print*,' input latitude of site in range [-90,90] (e.g. -57.35) :'
      read(5,*)alat
      if(alat.lt.-90.or.alat.gt.90)then
        print*,' Incorrect latitude'
        stop
      endif
      print*,' input longitude of site  (e.g. 97.35) :'
      read(5,*)alon

      elseif(iresp.eq.2)then
      print*,' input latitude degrees of site in range [-90,90] :'
      read(5,*)alat
      if(alat.lt.-90.or.alat.gt.90)then
        print*,' Incorrect latitude'
        stop
      endif
      print*,' input latitude minutes of site in range [0,60] :'
      read(5,*)alatm
      if(alatm.lt.0.or.alatm.gt.60)then
        print*,' Incorrect latitude minutes'
        stop
      endif
      alat=alat+sign(1,alat)*alatm/60.
      print*,sign(1,alat),alat
      print*,' input longitude degrees of site  :'
      read(5,*)alon
      print*,' input longitude minutes of site in range [0,60] :'
      read(5,*)alonm
      if(alonm.lt.0.or.alonm.gt.60)then
        print*,' Incorrect longitude minutes'
        stop
      endif
      alon=alon+sign(1,alon)*alonm/60.

      else
      print*,' must chose 1 or 2'
      stop
      endif

      print*,' input altitude of site in km  :'
      read(5,*)alt
      icoord = 1
c**********************************************************************

      theta = (90.0-alat)*fac
      phi   = alon*fac

c-----
c     input current model

      read(1,*)
      read(1,*) lm,ns,(tknts(k),k=1,ns+4)
      read(1,*) gt

      ts=tknts(4)
      tf=tknts(ns+1)
c-----
c     transform coordinates to geocentric colatitude,longitude,radius
c     if neccessary

      theta = (90.0-alat)*fac
      phi   = alon*fac
      if(icoord.eq.1)call coords(alt,theta,rad,sd,cd)

      sinth=sin(theta)
      costh=cos(theta)
      call plmbar(p,dp,costh,lmax,inorm)

c-----
c     calculate main field coefficients at time time
c
      do itime=0,200
      time=ts+itime/200.*(tf-ts)
c      print*,time

      call interv(tknts,time,nspl,nleft)
      call bspline(tknts,time,nspl,jord,nleft,spl(nleft-3))

      do  k=1,n
       g(k)=0.0
       do j=1,4
        g(k) = g(k) + spl(j+nleft-4)*gt(k,j+nleft-4)
       enddo
      enddo

c-----
c     calculate secular variation coefficients at time time
c
      call bspline1(tknts,time,nspl,nleft,spl)

      do  k=1,n
       gd(k)=0.0
       do j=1,4
        gd(k) = gd(k) + spl(j+nleft-4)*gt(k,j+nleft-4)
       enddo
      enddo

c     calculate main field elements

      call magfdz(p,dp,theta,phi,rad,lmax,g,dx,dy,dz,x,y,z,h,f,
     >ainc,d,icoord,sd,cd)

c     convert inclination and declination to degrees

      ainc=ainc/fac
      d=d/fac

c     calculate secular variation field elements

      xd = 0.
      yd = 0.
      zd = 0.
      do k=1,n
        xd = xd + dx(k)*gd(k)
        yd = yd + dy(k)*gd(k)
        zd = zd + dz(k)*gd(k)
      enddo

      hd = (x*xd + y*yd)/h
      fd = (x*xd + y*yd + z*zd)/f
      dd = (x*yd - y*xd)/h**2
      aincd = (h*zd - z*hd)/f**2

c     convert inclination and declination sv to degrees/yr

      aincd=aincd/fac
      dd=dd/fac

      write(7,*)time,d
      write(8,*)time,ainc
      enddo
      stop
      end

      subroutine interv(tknts,time,nspl,nleft)

      dimension tknts(nspl+1)

c  interv3 altered from subroutine interv to treat cubic splines where
c  knots equally spaced (no stripping of repeated knots needed)
c  interval starts at tknts(4) and ends at tknts(nspl+1)

c-----
c    calculate nleft:
c                  tknts(nleft) < tknts(nleft+1)
c                  tknts(nleft) <= time <= tknts(nleft+1)
c

c  check we are in-range

      if(time.lt.tknts(4).or.time.gt.tknts(nspl+1)) return

      do 200 n=5,nspl+1
       if(time.le.tknts(n)) then
        nleft=n-1
        goto 210
       endif
200   continue
210   continue

      return
      end
       subroutine bspline(tknts,t,nspl,jorder,nleft,spl)

c calculate splines of order jorder where 1 <= jorder <= 4

       dimension tknts(nspl+4)
       dimension spl(4)

       dimension deltal(4),deltar(4)

       spl(1)=1.0

       do 200 j=1,jorder-1

       deltar(j) = tknts(nleft+j) - t
       deltal(j) = t - tknts(nleft+1-j)
       saved=0.0

       do 100 i=1,j
        term = spl(i)/(deltar(i)+deltal(j+1-i))
        spl(i) = saved + deltar(i)*term
        saved = deltal(j+1-i)*term
100    continue

       spl(j+1) = saved

200    continue

       return
       end

      subroutine coords(h,theta,r,sd,cd)
      pi=3.14159265
      b1=40680925.
      b2=40408585.
      theta=pi/2-theta
      clat=cos(theta)
      slat=sin(theta)
      one=b1*clat*clat
      two=b2*slat*slat
      three=one+two
      four=sqrt(three)
      r=sqrt(h*(h+2.*four)+(b1*one+b2*two)/three)
      cd=(h+four)/r
      sd=(b1-b2)/four*slat*clat/r
      sinth=slat*cd-clat*sd
      costh=clat*cd+slat*sd
      theta=pi/2.-atan2(sinth,costh)
      return
      end
      subroutine magfdz(p,dp,theta,phi,r,lmax,g,dx,dy,dz,
     >x,y,z,h,f,i,d,
     >igeo,sd,cd)

c
c***************************************************************
c
c     j bloxham  8 nov 1982 & 11 oct 1983
c
c     modified version of dg13.sv.fort:magfd & jb62.sv.progs:magfds
c
c     gives field components at radius r
c
c     if igeo=1 then the field elements are rotated to
c     the local geodetic frame at the point specified by
c     (r,theta,phi) in the spherical coordinate system.
c
c***************************************************************
c
c
c     this version 16 jan 87
c
c     saves dx dy dz in computation
c
cc======================================================================

      dimension g(lmax*(lmax+2))
      dimension dx(lmax*(lmax+2)),dy(lmax*(lmax+2)),dz(lmax*(lmax+2))
      dimension p((lmax+1)*(lmax+2)/2),dp((lmax+1)*(lmax+2)/2)
      real i

      b=6371.2/r
      x=0.
      y=0.
      z=0.
      sinth=sin(theta)
      if(abs(sinth).lt.1.e-10) sinth=1.e-10

      do 20 l=1,lmax

      l1=l+1
      bb=b**(l+2)
      k=l*l
      k1=(l*l1)/2+1

      dx(k)=dp(k1)*bb
      dy(k)=0.
      dz(k)=-p(k1)*l1*bb
      x=x+g(k)*dx(k)
      z=z+g(k)*dz(k)

      do 20 m=1,l

      t=float(m)*phi
      k=l*l+2*m-1
      k1=(l*l1)/2+m+1
      sint=sin(t)
      cost=cos(t)

      dxd = dp(k1)*bb
      dx(k) = dxd*cost
      dx(k+1) = dxd*sint
      x = x + (g(k)*dx(k)) + (g(k+1)*dx(k+1))

      dxy = m*p(k1)*bb/sinth
      dy(k) = dxy*sint
      dy(k+1) = -dxy*cost
      y = y + (g(k)*dy(k)) + (g(k+1)*dy(k+1))

      dzd = -l1*p(k1)*bb
      dz(k) = dzd*cost
      dz(k+1) = dzd*sint
      z = z + (g(k)*dz(k)) + (g(k+1)*dz(k+1))

20    continue

      if(igeo.eq.1) then

      xs = x
      x = x*cd + z*sd
      z = z*cd - xs*sd

      do 50 k=1,lmax*(lmax+2)
      dxk = dx(k)
      dzk = dz(k)
      dx(k) = dxk*cd + dzk*sd
      dz(k) = dzk*cd - dxk*sd
50    continue

      endif

      h=sqrt(x*x+y*y)
      f=sqrt(h*h+z*z)
      i=asin(z/f)
      d=atan2(y,x)

      return
      end

      subroutine plmbar(p,dp,z,lmax,inorm)
c
c  evaluates normalized associated legendre function p(l,m) as function of
c   z=cos(colatitude) using recurrence relation starting with p(l,l)
c   and then increasing l keeping m fixed.  normalization is:
c   integral(y(l,m)*y(l,m))=4.*pi, where y(l,m) = p(l,m)*exp(i*m*longitude),
c   which is incorporated into the recurrence relation. p(k) contains p(l,m)
c   with k=(l+1)*l/2+m+1; i.e. m increments through range 0 to l before
c   incrementing l. routine is stable in single and double precision to
c   l,m = 511 at least; timing proportional to lmax**2
c   r.j.o'connell 7 sept. 1989

c   a.jackson 19 october 1989  code added at end:
c   (1) choice of normalisation added ---
c       if inorm = 1 schmidt normalisation is chosen, where
c       p[schmidt](l,m) = 1/sqrt(2*l+1)*p[normalised](l,m)
c       inorm = 2 for fully normalised harmonics, otherwise error
c   (2) derivatives added and stored in dp(k)
c       using same arrangement as for p(k)
c
      dimension p(*),dp(*)
c     --dimension of p, dp must be (lmax+1)*(lmax+2)/2 in calling program
      if (lmax.lt.0.or.abs(z).gt.1.d0) pause 'bad arguments'
      if(inorm.ne.1.and.inorm.ne.2)then
        write(6,99)
99      format('inorm incorrect: use inorm=1 for schmidt normalisation'
     >  ,/,    '                     inorm=2 for full normalisation')
        stop
      endif

c       --case for p(l,0)
	pm2=1.d0
	p(1)=1.d0
	dp(1)=0.d0
	if (lmax .eq. 0) return
	pm1=z
	p(2)=dsqrt(3.d0)*pm1
	k=2
	do 4 l=2,lmax
	  k=k+l
	  plm=(dfloat(2*l-1)*z*pm1-dfloat(l-1)*pm2)/dfloat(l)
	  p(k)=dsqrt(dfloat(2*l+1))*plm
	  pm2=pm1
4	  pm1=plm
c       --case for m > 0
	pmm = 1.d0
        sintsq = (1.d0-z)*(1.d0+z)
        fnum = -1.0d0
        fden = 0.0d0
	kstart = 1
	do 20 m =1 ,lmax
c         --case for p(m,m)
	  kstart = kstart+m+1
          fnum = fnum+2.0d0
          fden = fden+2.0d0
          pmm = pmm*sintsq*fnum/fden
	  pm2 = dsqrt(dfloat(4*m+2)*pmm)
	  p(kstart) = pm2
	  if (m .eq. lmax) goto 100
c         --case for p(m+1,m)
	  pm1=z*dsqrt(dfloat(2*m+3))*pm2
	  k = kstart+m+1
	  p(k) = pm1
c         --case for p(l,m) with l > m+1
	  if (m .lt. (lmax-1)) then
	   do 10 l = m+2,lmax
	    k = k+l
	    f1=dsqrt(dfloat((2*l+1)*(2*l-1))/dfloat((l+m)*(l-m)))
            f2=dsqrt(dfloat((2*l+1)*(l-m-1)*(l+m-1))
     &              /dfloat((2*l-3)*(l+m)*(l-m)))
            plm=z*f1*pm1-f2*pm2
	    p(k) = plm
	    pm2 = pm1
10	    pm1 = plm
	  endif
20	  continue

100     continue

c       choice of normalisation:
        if(inorm.eq.1)then
        k=1
        do 30 l=1,lmax
        fac=1.d0/dsqrt(dfloat(2*l+1))
        do 30 m=0,l
        k=k+1
        p(k)=p(k)*fac
30      continue
        endif

c       now find derivatives of p(z) wrt theta, where z=cos(theta)
c       recurrence is same regardless of normalisation, since l=constant
        dp(2)=-p(3)
        dp(3)=p(2)
        k=3
        do 200 l=2,lmax

          k=k+1
c         treat m=0 and m=l separately
          dp(k)=-dsqrt(dfloat(l*(l+1))/2.d0)*p(k+1)
          dp(k+l)=dsqrt(dfloat(l)/2.d0)*p(k+l-1)
          do 300 m=1,l-1
            k=k+1
            fac1=dsqrt( dfloat( (l-m)*(l+m+1) ) )
            fac2=dsqrt( dfloat( (l+m)*(l-m+1) ) )
            if(m.eq.1)fac2=fac2*dsqrt(2.d0)
            dp(k)=0.5d0*( fac2*p(k-1) - fac1*p(k+1) )
300       continue
          k=k+1

200     continue
	return
	end
      subroutine bspline1(tknts,t,nspl,nl,spl)

      parameter(nsplmx=200)
      dimension tknts(nspl+4)
      dimension spl(nspl),spl1(nsplmx)

      data jord/3/

      if(nspl.gt.nsplmx)then
      write(6,*)' increase dimensions of spl1 in bspline1'
      stop
      endif

      do 100 is=1,nspl
       spl(is)=0.0
100   continue

      call bspline(tknts,t,nspl,jord,nl,spl1(nl-2))

      spl(nl) = b0(tknts,nl)*spl1(nl)

      spl(nl-1) = b0(tknts,nl-1)*spl1(nl-1) + b1(tknts,nl-1)*spl1(nl)
      spl(nl-2) = b0(tknts,nl-2)*spl1(nl-2) + b1(tknts,nl-2)*spl1(nl-1)
      spl(nl-3) = b1(tknts,nl-3)*spl1(nl-2)

      return
      end
      function b0(tknts,i)
      dimension tknts(*)

      b0 = 3.0/( tknts(i+3) - tknts(i) )

      return
      end

      function b1(tknts,i)
      dimension tknts(*)

      b1 = -3.0/( tknts(i+4) - tknts(i+1) )

      return
      end