subroutine magfld( jd, p, b )
!
!     Given date, time, and a geocentric position (direction,
!     distance) return the Earth's magnetic (main) field vector b
!     (direction, magnitude).
!
!     B. Knapp, 1998-11-19, 2000-12-28, 2001-04-18 (wmm2000),
!               2004-03-25 (external cosd, sind, acosd, atan2d),
!               2006-03-27 (wmm2005)
!
!     Input
!
!        jd - Julian Day Number and fractional day (UT).  Note that
!             precision must be sufficient for an accurate UT time
!             of day, because the geomagnetic field depends on the
!             geographic longitude.
!
!        p  - Geocentric position, in Earth-centered inertial frame
!             rectangular coordinates, an array of length 4, with the
!             first three elements a unit vector giving the direction
!             (x,y,z) and the fourth element the distance from the
!             center of the Earth (in km).
!
!     Output
!
!        b  - The Earth's magnetic field vector, an array of length
!             4, with the first three elements a unit vector giving
!             the field direction (in the ECI frame) and the fourth
!             the field magnitude (in nTesla).
!
!
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      implicit none
!
!     Input
      real*8 jd, p(4)
!
!     Output
      real*8 b(4)
!
!     Local
      real*8 gst, lst, lat, lon, xyz(3), bmag, clat, slat, clon, slon
      integer*4 i
!
!     External functions, subroutines
      real*8 cosd, sind, acosd, atan2d, sidtim
      external cosd, sind, acosd, atan2d, wmm2005, sidtim
!
!
!     Local sidereal time (i.e., Right Ascension)
      lst = atan2d( p(2), p(1) )
      if ( lst .lt. 0.d0 ) lst = lst + 360.d0
!
!     Sidereal time at Greenwich
      gst = sidtim( jd )

!     Geocentric latitude (i.e., declination) and longitude
      lat = acosd( -p(3) ) - 90.d0
      lon = mod( gst - lst, 360.d0 )
      if ( lon .lt.   0.d0 ) lon = lon + 360.d0
      if ( lon .gt. 180.d0 ) lon = lon - 360.d0
!
!     Get the magnetic field in the local horizon coordinate frame
!     (x = north, y = east, z = nadir)
      call wmm2005( jd, lat, -lon, p(4), xyz )  !(pass east longitude)
!
!     Get the field magnitude
      bmag = 0.d0
      do i=1,3
         bmag = bmag + xyz(i)**2
      enddo
      bmag = sqrt( bmag )
!
!     Normalize xyz
      do i=1,3
         xyz(i) = xyz(i)/bmag
      enddo
!
!     Rotate to the ECI (equatorial) frame
      clon = cosd( lst )
      slon = sind( lst )
      clat = cosd( lat )
      slat = sind( lat )
      b(1) = -slat*clon*xyz(1) - slon*xyz(2) - clat*clon*xyz(3)
      b(2) = -slat*slon*xyz(1) + clon*xyz(2) - clat*slon*xyz(3)
      b(3) =       clat*xyz(1)               -      slat*xyz(3)
      b(4) = bmag
!
      return
      end