!***************************************************************************************************
! Copyright 1987, D. Bilitza, 2007 S. Bourdarie
!
! This file is part of IRBEM-LIB.
!
! IRBEM-LIB is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! IRBEM-LIB is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with IRBEM-LIB. If not, see .
!
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE IGRF(XXX,YYY,ZZZ,BXXX,BYYY,BZZZ)
C-------------------------------------------------------------------
C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL
C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61,
C 1970.
C--------------------------------------------------------------------
C CHANGES (D. BILITZA, NOV 87):
C - FIELD COEFFICIENTS IN BINARY DATA FILES INSTEAD OF BLOCK DATA
C - CALCULATES DIPOL MOMENT
C--------------------------------------------------------------------
C INPUT: ENTRY POINT FELDG
C GLAT GEODETIC LATITUDE IN DEGREES (NORTH)
C GLON GEODETIC LONGITUDE IN DEGREES (EAST)
C ALT ALTITUDE IN KM ABOVE SEA LEVEL
C
C ENTRY POINT FELDC
C V(3) CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM)
C X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE
C Y-AXIS POINTING TO EQUATOR AT 90 LONG.
C Z-AXIS POINTING TO NORTH POLE
C
C COMMON BLANK AND ENTRY POINT FELDI ARE NEEDED WHEN USED
C IN CONNECTION WITH L-CALCULATION PROGRAM SHELLG.
C
C COMMON /MODEL/ AND /GENER/
C UMR = ATAN(1.0)*4./180. *UMR=
C ERA EARTH RADIUS FOR NORMALIZATION OF CARTESIAN
C COORDINATES (6371.2 KM)
C AQUAD, BQUAD SQUARE OF MAJOR AND MINOR HALF AXIS FOR
C EARTH ELLIPSOID AS RECOMMENDED BY INTERNATIONAL
C ASTRONOMICAL UNION (6378.160, 6356.775 KM).
C NMAX MAXIMUM ORDER OF SPHERICAL HARMONICS
C TIME YEAR (DECIMAL: 1973.5) FOR WHICH MAGNETIC
C FIELD IS TO BE CALCULATED
C G(M) NORMALIZED FIELD COEFFICIENTS (SEE FELDCOF)
C M=NMAX*(NMAX+2)
C------------------------------------------------------------------------
C OUTPUT: BABS MAGNETIC FIELD STRENGTH IN GAUSS
C BNORTH, BEAST, BDOWN COMPONENTS OF THE FIELD WITH RESPECT
C TO THE LOCAL GEODETIC COORDINATE SYSTEM, WITH AXIS
C POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
C AND DOWNWARD.
C-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER*4 NMAX,IHMAX,LAST,I,K,IH,M,IMAX,IL
REAL*8 XXX,YYY,ZZZ
REAL*8 BXXX,BYYY,BZZZ
REAL*8 F,RQ,X,Y,Z,S,T
REAL*8 XI(3),H(144)
REAL*8 G(144)
COMMON/MODEL/ G
C
C-- IS RECORDS ENTRY POINT
C
C*****ENTRY POINT FELDG TO BE USED WITH GEODETIC CO-ORDINATES
RQ=1.D0/(XXX*XXX+YYY*YYY+ZZZ*ZZZ)
XI(1)=XXX*RQ
XI(2)=YYY*RQ
XI(3)=ZZZ*RQ
C
NMAX = 10
IHMAX=NMAX*NMAX+1
LAST=IHMAX+NMAX+NMAX
IMAX=NMAX+NMAX-1
DO 8 I=IHMAX,LAST
8 H(I)=G(I)
DO 6 K=1,3,2
I=IMAX
IH=IHMAX
1 IL=IH-I
F=2.D0/((I-K+2)*1.D0)
X=XI(1)*F
Y=XI(2)*F
Z=XI(3)*(F+F)
I=I-2
IF(I-1)5,4,2
2 DO 3 M=3,I,2
H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)-H(IH+M-1))
A -Y*(H(IH+M+2)+H(IH+M-2))
3 H(IL+M)=G(IL+M)+Z*H(IH+M)+X*(H(IH+M+2)-H(IH+M-2))
A +Y*(H(IH+M+3)+H(IH+M-1))
4 H(IL+2)=G(IL+2)+Z*H(IH+2)+X*H(IH+4)-Y*(H(IH+3)+H(IH))
H(IL+1)=G(IL+1)+Z*H(IH+1)+Y*H(IH+4)+X*(H(IH+3)-H(IH))
5 H(IL)=G(IL)+Z*H(IH)+2.*(X*H(IH+1)+Y*H(IH+2))
IH=IL
IF(I.GE.K)GOTO1
6 CONTINUE
S=.5D0*H(1)+2.D0*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2))
T=(RQ+RQ)*SQRT(RQ)
BXXX=T*(H(3)-S*XXX)
BYYY=T*(H(4)-S*YYY)
BZZZ=T*(H(2)-S*ZZZ)
RETURN
END
C