|
|
C $Id: invar.f 3304 2011-01-17 15:25:59Z brideout $
|
|
|
C
|
|
|
SUBROUTINE INVAR(TM,FLAT,FLONG,ALT,ERR,BB,FL)
|
|
|
C
|
|
|
C Private/Internal subroutine. Part of Apex coordinate computation
|
|
|
C package. See COORD for public API. INVAR converts coordinates
|
|
|
C TM, FLAT, FLON and ALT to L-shell coordinates FL and BB. The
|
|
|
C uncertainty in FL is typically less than 10.*ERR*FL (percent)
|
|
|
C
|
|
|
C Input:
|
|
|
C TM - time in years for desired field (e.g. 1971.25)
|
|
|
C FLAT - geocentric latitude (degrees)
|
|
|
C FLONG - east longitude
|
|
|
C ALT - altitude (km)
|
|
|
C ERR - tolerance factor
|
|
|
C
|
|
|
C Output:
|
|
|
C BB - Magnetic Field strength at point.
|
|
|
C FL - McIlwain's L-shell parameter i.e.
|
|
|
C Invariant Latitude = ACOS(DSQRT(1.0D0/FL))/DTR
|
|
|
C
|
|
|
C .. Scalar Arguments ..
|
|
|
DOUBLE PRECISION ALT,BB,ERR,FL,FLAT,FLONG,TM
|
|
|
C ..
|
|
|
C .. Local Scalars ..
|
|
|
DOUBLE PRECISION ASUM,BCO,CCO,DCLT,DCO,DN,DX,FLINT,SA,SC
|
|
|
INTEGER I,J,JEP,JUP
|
|
|
C ..
|
|
|
C .. Local Arrays ..
|
|
|
DOUBLE PRECISION ARC(200),B(200),BEG(200),BEND(200),BLOG(200),
|
|
|
* ECO(200),R1(3),R2(3),R3(3),V(3,3),VN(3),VP(3)
|
|
|
C ..
|
|
|
C .. External Subroutines ..
|
|
|
EXTERNAL CARMEL,INTEG,LINES,STARTR
|
|
|
C ..
|
|
|
C .. Intrinsic Functions ..
|
|
|
INTRINSIC ABS,DCOS,DEXP,DLOG,DSQRT
|
|
|
C ..
|
|
|
V(1,2) = ALT/6371.2D0
|
|
|
V(2,2) = (90.D0-FLAT)/57.2957795D0
|
|
|
V(3,2) = FLONG/57.2957795D0
|
|
|
ARC(1) = 0.D0
|
|
|
ARC(2) = (1.0D0+V(1,2))*DSQRT(ERR)*0.3D0
|
|
|
DCLT = 1.5708D0 - 0.2007D0*DCOS(V(3,2)+1.239D0)
|
|
|
IF (V(2,2).GT.DCLT) ARC(2) = -ARC(2)
|
|
|
CALL STARTR(R1,R2,R3,B,ARC,V,TM)
|
|
|
DO 10 I = 1,3
|
|
|
VP(I) = V(I,2)
|
|
|
VN(I) = V(I,3)
|
|
|
10 CONTINUE
|
|
|
CALL LINES(R1,R2,R3,B,ARC,ERR,J,VP,VN,TM)
|
|
|
IF (J.LT.200) THEN
|
|
|
JUP = J
|
|
|
DO 20 J = 1,JUP
|
|
|
ARC(J) = ABS(ARC(J))
|
|
|
BLOG(J) = DLOG(B(J))
|
|
|
20 CONTINUE
|
|
|
JEP = JUP - 1
|
|
|
DO 30 J = 2,JEP
|
|
|
ASUM = ARC(J) + ARC(J+1)
|
|
|
DX = BLOG(J-1) - BLOG(J)
|
|
|
DN = ASUM*ARC(J)*ARC(J+1)
|
|
|
BCO = ((BLOG(J-1)-BLOG(J+1))*ARC(J)**2-DX*ASUM**2)/DN
|
|
|
CCO = (DX*ARC(J+1)-(BLOG(J)-BLOG(J+1))*ARC(J))/DN
|
|
|
SA = .75D0*ARC(J)
|
|
|
SC = SA + .25D0*ASUM
|
|
|
DCO = BLOG(J-1) - CCO*SA*SC
|
|
|
ECO(J) = BCO + CCO*(SA+SC)
|
|
|
BEG(J) = DEXP(DCO+ECO(J)*.5D0*ARC(J))
|
|
|
BEND(J) = DEXP(DCO+ECO(J)*.5D0*(ASUM+ARC(J)))
|
|
|
30 CONTINUE
|
|
|
BEG(JUP) = BEND(JEP)
|
|
|
BEND(JUP) = B(JUP)
|
|
|
ECO(JUP) = (2.0D0/ARC(JUP))*DLOG(BEND(JUP)/BEG(JUP))
|
|
|
CALL INTEG(ARC,BEG,BEND,B,JEP,ECO,FLINT)
|
|
|
CALL CARMEL(B(2),FLINT,FL)
|
|
|
ELSE
|
|
|
FL = -1.0D0
|
|
|
END IF
|
|
|
BB = B(2)
|
|
|
RETURN
|
|
|
C
|
|
|
END
|
|
|
|
