|
|
C $Id: iterat.f 3304 2011-01-17 15:25:59Z brideout $
|
|
|
C
|
|
|
SUBROUTINE ITERAT
|
|
|
C
|
|
|
C Private/Internal subroutine. Part of Apex coordinate computation
|
|
|
C package. See COORD for public API. ITERAT integrates magnetic
|
|
|
C field line using a 4-point adams formula after initialization.
|
|
|
C First 7 iterations advance point by 3*DS.
|
|
|
C
|
|
|
C Input (via common block ITER):
|
|
|
C L - step count. set l=1 first time thru,
|
|
|
C set l=l+1 thereafter.
|
|
|
C B,BR,BT,BP - field + components at point y
|
|
|
C ST - sine of geocentric colatitude
|
|
|
C SGN - sgn=+1 traces in direction of field
|
|
|
C sgn=-1 traces in negative field direction
|
|
|
C DS - integration stepsize (arc increment) in km
|
|
|
C
|
|
|
C Input, Output (via common block ITER):
|
|
|
C Y - R,DLAT,DLON: geocentric tracing point coordinates
|
|
|
C (km,deg)
|
|
|
C
|
|
|
C Output (via common block ITER):
|
|
|
C YOLD - Y at iteration L-1
|
|
|
C
|
|
|
C .. Scalars in Common ..
|
|
|
C ..
|
|
|
C .. Arrays in Common ..
|
|
|
C ..
|
|
|
C .. Local Scalars ..
|
|
|
C
|
|
|
DOUBLE PRECISION D12,D2,D24,D6,FAC,RAD
|
|
|
INTEGER I,J
|
|
|
C ..
|
|
|
C .. Local Arrays ..
|
|
|
DOUBLE PRECISION YP(3,4)
|
|
|
C ..
|
|
|
C .. Common blocks ..
|
|
|
COMMON /ITER/Y,YOLD,HP,BR,BT,BP,B,ST,SGN,DS,L
|
|
|
DOUBLE PRECISION B,BP,BR,BT,DS,HP,SGN,ST
|
|
|
INTEGER L
|
|
|
DOUBLE PRECISION Y(3),YOLD(3)
|
|
|
C ..
|
|
|
C .. Save statement ..
|
|
|
SAVE D12,D2,D24,D6,FAC,RAD,I,J,YP
|
|
|
C ..
|
|
|
C .. Data statements ..
|
|
|
DATA RAD/57.2957795D0/,D12/0.0D0/,D2/0.0D0/,D24/0.0D0/,D6/0.0D0/
|
|
|
C ..
|
|
|
C
|
|
|
YP(1,4) = SGN*BR/B
|
|
|
FAC = SGN*RAD/(B*Y(1))
|
|
|
YP(2,4) = -BT*FAC
|
|
|
YP(3,4) = BP*FAC/ST
|
|
|
IF (L.GT.7) THEN
|
|
|
DO 20 I = 1,3
|
|
|
YOLD(I) = Y(I)
|
|
|
Y(I) = YOLD(I) + D24*(55.D0*YP(I,4)-59.D0*YP(I,3)+
|
|
|
* 37.D0*YP(I,2)-9.D0*YP(I,1))
|
|
|
DO 10 J = 1,3
|
|
|
YP(I,J) = YP(I,J+1)
|
|
|
10 CONTINUE
|
|
|
20 CONTINUE
|
|
|
ELSE
|
|
|
DO 30 I = 1,3
|
|
|
IF (L.EQ.2) THEN
|
|
|
YP(I,2) = YP(I,4)
|
|
|
Y(I) = YOLD(I) + D2*(YP(I,2)+YP(I,1))
|
|
|
ELSE IF (L.EQ.3) THEN
|
|
|
Y(I) = YOLD(I) + D6*(2.D0*YP(I,4)+YP(I,2)+3.D0*YP(I,1))
|
|
|
ELSE IF (L.EQ.4) THEN
|
|
|
YP(I,2) = YP(I,4)
|
|
|
YOLD(I) = Y(I)
|
|
|
Y(I) = YOLD(I) + D2*(3.D0*YP(I,2)-YP(I,1))
|
|
|
ELSE IF (L.EQ.5) THEN
|
|
|
Y(I) = YOLD(I) + D12*(5.D0*YP(I,4)+8.D0*YP(I,2)-YP(I,1))
|
|
|
ELSE IF (L.EQ.6) THEN
|
|
|
YP(I,3) = YP(I,4)
|
|
|
YOLD(I) = Y(I)
|
|
|
Y(I) = YOLD(I) + D12*(23.D0*YP(I,3)-16.D0*YP(I,2)+
|
|
|
* 5.D0*YP(I,1))
|
|
|
ELSE IF (L.EQ.7) THEN
|
|
|
Y(I) = YOLD(I) + D24*(9.D0*YP(I,4)+19.D0*YP(I,3)-
|
|
|
* 5.D0*YP(I,2)+YP(I,1))
|
|
|
ELSE
|
|
|
D2 = DS/2.D0
|
|
|
D6 = DS/6.D0
|
|
|
D12 = DS/12.D0
|
|
|
D24 = DS/24.D0
|
|
|
YP(I,1) = YP(I,4)
|
|
|
YOLD(I) = Y(I)
|
|
|
Y(I) = YOLD(I) + DS*YP(I,1)
|
|
|
END IF
|
|
|
30 CONTINUE
|
|
|
END IF
|
|
|
RETURN
|
|
|
C
|
|
|
END
|
|
|
|
