startr.f
94 lines
| 2.5 KiB
| text/x-fortran
|
FortranFixedLexer
| r0 | 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 | ||||
