iterat.f
98 lines
| 3.1 KiB
| text/x-fortran
|
FortranFixedLexer
r0 | 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 | ||||