|
|
C $Id: startr.f 3304 2011-01-17 15:25:59Z brideout $
|
|
|
C
|
|
|
SUBROUTINE STARTR(R1,R2,R3,B,ARC,V,TM)
|
|
|
C
|
|
|
C Private/Internal subroutine. Part of Apex coordinate computation
|
|
|
C package. See COORD for public API.
|
|
|
C
|
|
|
C Input:
|
|
|
C TM - time in years for desired field (e.g. 1971.25)
|
|
|
C
|
|
|
C Input, Output:
|
|
|
C ARC - Altitudes in earth radii (array).
|
|
|
C V - floating point array.
|
|
|
C
|
|
|
C Output:
|
|
|
C R1,R2,R3 - field strength in geocentric directions.
|
|
|
C B - Magnitude of field (array)
|
|
|
C
|
|
|
C .. Scalar Arguments ..
|
|
|
DOUBLE PRECISION TM
|
|
|
C ..
|
|
|
C .. Array Arguments ..
|
|
|
DOUBLE PRECISION ARC(200),B(200),R1(3),R2(3),R3(3),V(3,3)
|
|
|
C ..
|
|
|
C .. Local Scalars ..
|
|
|
DOUBLE PRECISION AER,BP,BR,BT,COP,COT,DN,OER,RKM,SIP,SIT,SSQ
|
|
|
INTEGER I,IS
|
|
|
C ..
|
|
|
C .. External Subroutines ..
|
|
|
EXTERNAL MILMAG
|
|
|
C ..
|
|
|
C .. Intrinsic Functions ..
|
|
|
INTRINSIC ABS,DCOS,DSIN
|
|
|
C ..
|
|
|
SIT = ABS(DSIN(V(2,2)))
|
|
|
AER = V(1,2)
|
|
|
SSQ = SIT*SIT
|
|
|
OER = (6356.912D0+SSQ*(21.3677D0+.108D0*SSQ))/6371.2D0
|
|
|
V(1,2) = AER + OER
|
|
|
10 CONTINUE
|
|
|
IF (V(3,2).LT.0) THEN
|
|
|
V(3,2) = V(3,2) + 6.283185307D0
|
|
|
GO TO 10
|
|
|
END IF
|
|
|
RKM = V(1,2)*6371.2D0
|
|
|
C if(model.eq.6) rkm=rkm+14.288-ssq*(21.3677+.108*ssq)
|
|
|
COT = DCOS(V(2,2))
|
|
|
SIP = DSIN(V(3,2))
|
|
|
COP = DCOS(V(3,2))
|
|
|
CALL MILMAG(TM,RKM,SIT,COT,SIP,COP,BR,BT,BP,B(2))
|
|
|
R2(1) = BR/B(2)
|
|
|
DN = B(2)*V(1,2)
|
|
|
R2(2) = BT/DN
|
|
|
R2(3) = BP/(DN*SIT)
|
|
|
IS = 0
|
|
|
20 CONTINUE
|
|
|
DO 30 I = 1,3
|
|
|
V(I,1) = V(I,2) - ARC(2)*R2(I)
|
|
|
30 CONTINUE
|
|
|
SIT = ABS(DSIN(V(2,1)))
|
|
|
40 CONTINUE
|
|
|
RKM = V(1,1)*6371.2D0
|
|
|
SSQ = SIT*SIT
|
|
|
C if(model.eq.6) rkm=rkm+14.288-ssq*(21.3677+.108*ssq)
|
|
|
COT = DCOS(V(2,1))
|
|
|
SIP = DSIN(V(3,1))
|
|
|
COP = DCOS(V(3,1))
|
|
|
CALL MILMAG(TM,RKM,SIT,COT,SIP,COP,BR,BT,BP,B(1))
|
|
|
IF (B(1).GE.B(2)) THEN
|
|
|
R1(1) = BR/B(1)
|
|
|
ARC(3) = ARC(2)
|
|
|
DN = B(1)*V(1,1)
|
|
|
R1(2) = BT/DN
|
|
|
R1(3) = BP/(DN*SIT)
|
|
|
DO 50 I = 1,3
|
|
|
V(I,1) = V(I,2) - ARC(2)*(R1(I)+R2(I))/2.D0
|
|
|
50 CONTINUE
|
|
|
SIT = ABS(DSIN(V(2,1)))
|
|
|
IS = IS + 1
|
|
|
IF (IS.EQ.1) GO TO 40
|
|
|
GO TO 60
|
|
|
END IF
|
|
|
ARC(2) = -ARC(2)
|
|
|
GO TO 20
|
|
|
60 CONTINUE
|
|
|
DO 70 I = 1,3
|
|
|
V(I,3) = V(I,2) + ARC(3)*((1.5D0)*R2(I)-.5D0*R1(I))
|
|
|
70 CONTINUE
|
|
|
R3(1) = 0.0D0
|
|
|
R3(2) = 0.0D0
|
|
|
R3(3) = 0.0D0
|
|
|
RETURN
|
|
|
C
|
|
|
END
|
|
|
|
