lines.f
168 lines
| 4.6 KiB
| text/x-fortran
|
FortranFixedLexer
| r0 | C $Id: lines.f 3304 2011-01-17 15:25:59Z brideout $ | |||
| C | ||||
| SUBROUTINE LINES(R1,R2,R3,B,ARC,ERR,J,VP,VN,TM) | ||||
| C | ||||
| C Private/Internal subroutine. Part of Apex coordinate computation | ||||
| C package. See COORD for public API. Makes repeated calls to the | ||||
| C IGRF, tracing magnetic field line to minimum B. | ||||
| C | ||||
| C Input: | ||||
| C ERR - tolerance factor (see INVAR) | ||||
| C TM - Time in floating point years (e.g. 1995.7) | ||||
| C | ||||
| C Input, Output: | ||||
| C R1,R2,R3 - field strength in geocentric directions. | ||||
| C B - Magnitude of field (array) | ||||
| C ARC - Altitudes in earth radii (array) | ||||
| C VP - Geocentric latitude (array) | ||||
| C VN - Geocentric longitude (array) | ||||
| C | ||||
| C Output: | ||||
| C J - Number of points in trace. | ||||
| C | ||||
| C .. Scalar Arguments .. | ||||
| DOUBLE PRECISION ERR,TM | ||||
| INTEGER J | ||||
| C .. | ||||
| C .. Array Arguments .. | ||||
| DOUBLE PRECISION ARC(200),B(200),R1(3),R2(3),R3(3),VN(3),VP(3) | ||||
| C .. | ||||
| C .. Local Scalars .. | ||||
| DOUBLE PRECISION A1,A2,A3,AA,AAB,AD,AM,AO6,ARCJ,ASUM,BB,BD,BP,BR, | ||||
| * BT,CC,CD,COP,COT,CRE,DD,DN,PRE1,PRE2,PRE3,QRT, | ||||
| * RBAR,RKM,RT,SIP,SIT,SNA,X | ||||
| INTEGER I,ILP,IS | ||||
| C .. | ||||
| C .. Local Arrays .. | ||||
| DOUBLE PRECISION RA(3) | ||||
| C .. | ||||
| C .. External Subroutines .. | ||||
| EXTERNAL MILMAG | ||||
| C .. | ||||
| C .. Intrinsic Functions .. | ||||
| INTRINSIC ABS,DCOS,DSIN,DSQRT | ||||
| C .. | ||||
| C .. Data statements .. | ||||
| DATA PRE1,PRE2,PRE3,RA/0.0D0,0.0D0,0.0D0,0.0D0,0.0D0,0.0D0/ | ||||
| C .. initialize | ||||
| AA=0.0D0 | ||||
| AD=0.0D0 | ||||
| ARCJ=0.0D0 | ||||
| ASUM=0.0D0 | ||||
| BB=0.0D0 | ||||
| BD=0.0D0 | ||||
| CC=0.0D0 | ||||
| CD=0.0D0 | ||||
| CRE = 0.25D0 | ||||
| IF (ERR.LT.0.15625D0) CRE = (ERR**0.333333333D0) | ||||
| A3 = ARC(3) | ||||
| AAB = ABS(A3) | ||||
| SNA = A3/AAB | ||||
| A1 = ARC(1) | ||||
| A2 = ARC(2) | ||||
| AO6 = A3*A3/6.0D0 | ||||
| J = 3 | ||||
| ILP = 1 | ||||
| IS = 1 | ||||
| 10 CONTINUE | ||||
| IF (VN(2).LT.0) VN(2) = -VN(2) | ||||
| 20 CONTINUE | ||||
| IF (VN(2).GT.3.141592653D0) THEN | ||||
| VN(2) = 6.283185307D0 - VN(2) | ||||
| GO TO 20 | ||||
| END IF | ||||
| 30 CONTINUE | ||||
| IF (VN(3).LT.0) THEN | ||||
| VN(3) = VN(3) + 6.283185307D0 | ||||
| GO TO 30 | ||||
| END IF | ||||
| 40 CONTINUE | ||||
| IF (VN(3).GT.6.283185307D0) THEN | ||||
| VN(3) = VN(3) - 6.283185307D0 | ||||
| GO TO 40 | ||||
| END IF | ||||
| IF (IS.EQ.2) THEN | ||||
| SIT = ABS(DSIN(VN(2))) | ||||
| B(J) = B(J)*((PRE1/VN(1))**3) | ||||
| QRT = 0.5D0*ABS(R3(1))/(0.1D0+ABS(R3(2)*VN(1))) | ||||
| X = (ABS(VN(1)-PRE1)+QRT*ABS(VN(1)*VN(2)-PRE2)+ | ||||
| * ABS(VN(1)*SIT*VN(3)-PRE3))/(AAB*ERR*DSQRT(1.D0+QRT*QRT)) | ||||
| IF (.NOT.(ILP.EQ.1.OR.ILP.EQ.3)) THEN | ||||
| IF (X.GE.3.3D0) THEN | ||||
| A3 = A3*0.2D0*(8.0D0+X)/(0.8D0+X) | ||||
| J = J - 1 | ||||
| ILP = 3 | ||||
| ASUM = A2 + A1 | ||||
| AA = ASUM*A1 | ||||
| BB = A2*A1 | ||||
| CC = ASUM*A2 | ||||
| DO 50 I = 1,3 | ||||
| VN(I) = VP(I) | ||||
| R3(I) = R2(I) | ||||
| R2(I) = R1(I) | ||||
| R1(I) = RA(I) | ||||
| 50 CONTINUE | ||||
| GO TO 60 | ||||
| END IF | ||||
| END IF | ||||
| IF (J.GE.200) GO TO 80 | ||||
| A1 = A2 | ||||
| IF (B(J).GT.B(2)) GO TO 80 | ||||
| ILP = 2 | ||||
| A2 = A3 | ||||
| A3 = A3*.2D0*(8.D0+X)/(.8D0+X) | ||||
| AM = (2.D0-R3(2)*VN(1))*VN(1)*CRE | ||||
| IF (ABS(A3).GT.AM) A3 = SNA*AM | ||||
| IF (SNA*R3(1)+.5D0.LE.0) THEN | ||||
| AM = -.5D0*SNA*VN(1)/R3(1) | ||||
| IF (ABS(A3).GT.AM) A3 = SNA*AM | ||||
| END IF | ||||
| 60 CONTINUE | ||||
| ARC(J+1) = A3 | ||||
| AAB = ABS(A3) | ||||
| IS = 1 | ||||
| J = J + 1 | ||||
| AO6 = A3*A3/6.0D0 | ||||
| ARCJ = A1 + A2 + A3 | ||||
| AD = (ASUM+A1)/AA | ||||
| BD = ASUM/BB | ||||
| CD = A1/CC | ||||
| ELSE | ||||
| SIT = ABS(DSIN(VN(2))) | ||||
| PRE1 = VN(1) | ||||
| PRE2 = PRE1*VN(2) | ||||
| PRE3 = PRE1*SIT*VN(3) | ||||
| RKM = VN(1)*6371.2D0 | ||||
| COT = DCOS(VN(2)) | ||||
| SIP = DSIN(VN(3)) | ||||
| COP = DCOS(VN(3)) | ||||
| CALL MILMAG(TM,RKM,SIT,COT,SIP,COP,BR,BT,BP,B(J)) | ||||
| R3(1) = BR/B(J) | ||||
| DN = B(J)*VN(1) | ||||
| R3(2) = BT/DN | ||||
| R3(3) = BP/(DN*SIT) | ||||
| ASUM = A3 + A2 | ||||
| AA = ASUM*A2 | ||||
| BB = A3*A2 | ||||
| CC = ASUM*A3 | ||||
| IS = 2 | ||||
| END IF | ||||
| DO 70 I = 1,3 | ||||
| DD = R1(I)/AA - R2(I)/BB + R3(I)/CC | ||||
| IF (IS.NE.2) THEN | ||||
| RT = R1(I) - (AD*R1(I)-BD*R2(I)+CD*R3(I)-DD*ARCJ)*ARCJ | ||||
| RA(I) = R1(I) | ||||
| R1(I) = R2(I) | ||||
| R2(I) = R3(I) | ||||
| R3(I) = RT | ||||
| VP(I) = VN(I) | ||||
| END IF | ||||
| RBAR = (R2(I)+R3(I))/2.D0 - DD*AO6 | ||||
| VN(I) = VP(I) + A3*RBAR | ||||
| 70 CONTINUE | ||||
| GO TO 10 | ||||
| 80 CONTINUE | ||||
| RETURN | ||||
| C | ||||
| END | ||||
