invar.f
83 lines
| 2.7 KiB
| text/x-fortran
|
FortranFixedLexer
| r0 | 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 | ||||
