coordff.f
289 lines
| 10.6 KiB
| text/x-fortran
|
FortranFixedLexer
r0 | C $Id: coordff.f 6720 2019-01-29 21:39:36Z brideout $ | |||
C | ||||
SUBROUTINE COORDFF(SLATGD,SLON,SR,SLATGC,TM,AZ,EL,RANGE, | ||||
* GDLAT,GLON,GDALT,QDIAG,QSPHERICAL,RCOR) | ||||
C Calculates the listed coordinates of a specified point. the point | ||||
C may be specified either by slatgd, slon, sr, slatgc, tm, az, el, | ||||
C range (range .ge. 0.) or by gdlat, glon, gdalt, tm (range .lt. | ||||
C 0.). Most of the output parameters involve the geomagnetic field | ||||
C as specified by the International Geomagnetic Referenc eField | ||||
C (IGRF) fo the specified epoch (tm). Apex latitude and longitude, | ||||
C along with either gdalt or arc form a geomagnetic coordinate | ||||
C system within either hemisphere. This system is necessarily | ||||
C non-orthogonal which accounts for the complexity of many of the | ||||
C calculations in coord. Apex latitude is always positive which is | ||||
C appropriate when using apex latitude and longitude to label | ||||
C field lines. | ||||
C | ||||
C This is an even faster version of coord.f used by getMag and faraday | ||||
C because coord is slower calculating parms not needed for | ||||
C those methods. RCOR(11-64) always return 0. | ||||
C | ||||
C Input: | ||||
C SLATGD - station geodetic latitude | ||||
C SLON - station longitude | ||||
C SR - radial distance of station from center of earth | ||||
C SLATGC - station geocentric latitude | ||||
C TM - time in years (e.g. 1975.2) | ||||
C AZ - radar azimuth | ||||
C EL - radar elevation | ||||
C RANGE - range to observation point | ||||
C QDIAG - Print diagnostic output on terminal if non-zero | ||||
C QSPHERICAL - Assume spherical coordinates instead of apex coordinates | ||||
C if non-zero. This helps verify the numerical | ||||
C nonorthogonal coordinate calculations because the | ||||
C spherical case also has an analytic solution. | ||||
C QDIAG and QSPHERICAL should be set to 0 except for development | ||||
C Input, output: | ||||
C GDLAT - geodetic latitude of observation point | ||||
C GLON - longitude of observation point | ||||
C GDALT - altitude above spheroid of observation point | ||||
C Output: | ||||
C RCOR(1) - az - radar azimuth | ||||
C RCOR(2) - el - radar elevation | ||||
C RCOR(3) - range - range to observation point | ||||
C RCOR(4) - gdlat - geodetic latitude of observation point | ||||
C RCOR(5) - glon - longitude of observation point | ||||
C RCOR(6) - gdalt - altitude above spheroid of observation point | ||||
C RCOR(7) - b - magnitude of geomagnetic field | ||||
C RCOR(8) - br - radial component of geomagnetic field | ||||
C RCOR(9) - bt - southward component of geomagnetic field | ||||
C RCOR(10) - bp - eastward component of geomagnetic field | ||||
C RCOR(11) - rlatm - dip latitude | ||||
C RCOR(12) - rlati - invariant latitude | ||||
C RCOR(13) - rl - magnetic l parameter | ||||
C RCOR(14) - alat - apex latitude | ||||
C RCOR(15) - alon - apex longitude | ||||
C | ||||
C RCOR(16) - g(1,1) magnetic coordinate system metric tensor, | ||||
C upper half stored row-wise | ||||
C RCOR(17) - g(2,1) " " | ||||
C RCOR(18) - g(2,1) " " | ||||
C RCOR(19) - g(2,1) " " | ||||
C | ||||
C RCOR(20) - ctab(1) - direction cosine of beam wrt geodetic | ||||
C south | ||||
C RCOR(21) - ctab(2) - direction cosine of beam wrt geodetic | ||||
C east | ||||
C RCOR(22) - ctab(3) - direction cosine of beam wrt geodetic | ||||
C altitude (up) | ||||
C | ||||
C RCOR(23) - ctab(4) - direction cosine of beam wrt apex | ||||
C latitude (northward) | ||||
C RCOR(24) - ctab(5) - direction cosine of beam wrt apex | ||||
C longitude (eastward) | ||||
C RCOR(25) - ctab(6) - direction cosine of beam wrt fieldline | ||||
C geomagnetic field (up) | ||||
C | ||||
C RCOR(26) - cost1 - x-direction cosine of a vector | ||||
C perpendicularto l.o.s. w/respect to apex | ||||
C coords. | ||||
C RCOR(27) - cost2 - y-direction cosine " " | ||||
C RCOR(28) - cost3 - z-direction cosine " " | ||||
C | ||||
C RCOR(29) - ainc - inclination of geomagnetic field | ||||
C RCOR(30) - dec - declination of geomagnetic field | ||||
C RCOR(31) - gclat - geocentric latitude | ||||
C RCOR(32) - aspct - aspect angle | ||||
C RCOR(33) - cplat - conjugate geocentric latitude | ||||
C RCOR(34) - gcdlat - conjugate geodetic latitude | ||||
C RCOR(35) - cplon conjugate longitude | ||||
C RCOR(36 - 64) are set to zero. Use coord.f for those | ||||
C | ||||
C .. Scalar Arguments .. | ||||
DOUBLE PRECISION AZ,EL,GDALT,GDLAT,GLON,RANGE,SLATGC,SLATGD,SLON, | ||||
* SR,TM | ||||
INTEGER QDIAG,QSPHERICAL | ||||
C .. | ||||
C .. Array Arguments .. | ||||
DOUBLE PRECISION RCOR(64) | ||||
C .. | ||||
C .. Local Scalars .. | ||||
DOUBLE PRECISION AINC,ALAT,ALON,ARAD,ARC,ASPCT,B,BB,BP,BPERPEASTM, | ||||
* BPERPSOUTHM,BR,BT,BX,BY,BZ,CALAT,CALON,CARAD, | ||||
* CARC,CASPCT,CGDALT,CGDLAT, | ||||
* COST1,COST2,COST3,CP,CPLAT, | ||||
* CPLON,CPRKM,CSP,CSR,CST,CT,CX,CY,CZ,D,DEC, | ||||
* DECGC,DTR, | ||||
* EASTM,EASTP,EASTR,EASTT,EASTX, | ||||
* EASTY,EASTZ, | ||||
* EFMAG,EFX,EFY,EFZ, | ||||
* GCALT,GCLAT,GDLATS, | ||||
* HB,HI,P,PERPEASTP,PERPEASTR,PERPEASTT,PERPSOUTHP, | ||||
* PERPSOUTHR,PERPSOUTHT,PFX,PFY,PFZ,PLAT,PLON, | ||||
* PRKM,PX,PY,PZ,RFP,RFR,RFT,RFX,RFY,RFZ, | ||||
* RKM,RL,RLATI,RLATM,RP,RR,RT,SOUTHM,SOUTHP, | ||||
* SOUTHR,SOUTHT,SOUTHX,SOUTHY,SOUTHZ,SP,ST,T, | ||||
* TM1,UPBP,UPBR,UPBT,UPP,UPR,UPT,UPX,UPY, | ||||
* UPZ,XB,XDUM,YB,ZB | ||||
INTEGER IER,INIT,ISTOP,NPR | ||||
C .. | ||||
C .. Local Arrays .. | ||||
DOUBLE PRECISION BF(3),BPERPEAST(3),BPERPSOUTH(3),CTAB(9), | ||||
* EAST(3),EAST2(3), | ||||
* EF(3), | ||||
* PERPEAST(3),PERPEAST2(3), | ||||
* PERPSOUTH(3),PERPSOUTH2(3),PF(3),Q0(3), | ||||
* RF(3),SOUTH(3),SOUTH2(3),UP(3), | ||||
* UP2(3),UPB(3),UPB2(3),ZEROM(3) | ||||
C .. | ||||
C .. External Functions .. | ||||
DOUBLE PRECISION SPROD,TPROD,VMAG | ||||
EXTERNAL SPROD,TPROD,VMAG | ||||
C .. | ||||
C .. External Subroutines .. | ||||
EXTERNAL CONVRT,CSCONV,GDV,GMET,GMETBV,INVAR,LINTRA,LOOK,MILMAG, | ||||
* MINV,MTRAN3,POINT,RPCART,UVECT,VCTCNV,VPROD,VRECIP | ||||
C .. | ||||
C .. Intrinsic Functions .. | ||||
INTRINSIC ABS,ACOS,ATAN,COS,DATAN2,DBLE,DSQRT,MAX,MIN, | ||||
* SIN,SQRT | ||||
C .. | ||||
C .. Statement Functions .. | ||||
DOUBLE PRECISION TAN | ||||
C .. | ||||
C .. Equivalences .. | ||||
EQUIVALENCE (RFX,RF(1)),(RFY,RF(2)),(RFZ,RF(3)) | ||||
EQUIVALENCE (PFX,PF(1)),(PFY,PF(2)),(PFZ,PF(3)) | ||||
EQUIVALENCE (BX,BF(1)),(BY,BF(2)),(BZ,BF(3)) | ||||
EQUIVALENCE (EFX,EF(1)),(EFY,EF(2)),(EFZ,EF(3)) | ||||
EQUIVALENCE (SOUTHX,SOUTH(1)),(SOUTHY,SOUTH(2)),(SOUTHZ,SOUTH(3)) | ||||
EQUIVALENCE (EASTX,EAST(1)),(EASTY,EAST(2)),(EASTZ,EAST(3)) | ||||
EQUIVALENCE (UPX,UP(1)),(UPY,UP(2)),(UPZ,UP(3)) | ||||
C .. | ||||
C .. Data statements .. | ||||
DATA DTR/0.0174532925199D0/ | ||||
C .. | ||||
C .. Constants | ||||
C err1 and err2 are error flags from geo-cgm.f | ||||
C .. Statement Function definitions .. | ||||
TAN(XDUM) = DSIN(XDUM)/(DCOS(XDUM)+1.D-38) | ||||
C .. | ||||
C .....Initialize..... | ||||
COST1 = 0.0D0 | ||||
COST2 = 0.0D0 | ||||
COST3 = 0.0D0 | ||||
TM1 = ABS(TM) | ||||
IF (RANGE.LT.0.0D0) THEN | ||||
C | ||||
C .....Entry when gdlat,glon,gdalt are given..... | ||||
C .....Calculation fails when gdlat is too close to 0 or 90..... | ||||
GDLATS = GDLAT | ||||
GDLAT = MIN(GDLAT,89.99D0) | ||||
GDLAT = MAX(GDLAT,-89.99D0) | ||||
IF (ABS(GDLAT).LT.0.01D0) THEN | ||||
GDLAT = 0.01D0 | ||||
END IF | ||||
CALL CONVRT(1,GDLAT,GDALT,GCLAT,RKM) | ||||
CALL LOOK(SR,SLATGC,SLON,RKM,GCLAT,GLON,AZ,EL,RANGE) | ||||
ELSE | ||||
C | ||||
C .....Entry when az, el, range are given..... | ||||
CALL POINT(SR,SLATGC,SLON,AZ,EL,RANGE,RKM,GCLAT,GLON) | ||||
CALL CONVRT(2,GDLAT,GDALT,GCLAT,RKM) | ||||
GDLATS = GDLAT | ||||
GDLAT = MIN(GDLAT,89.99D0) | ||||
GDLAT = MAX(GDLAT,-89.99D0) | ||||
IF (ABS(GDLAT).LT.0.01D0) THEN | ||||
GDLAT = 0.01D0 | ||||
END IF | ||||
END IF | ||||
C | ||||
C .....We now have all needed information on station, radar beam | ||||
C and location of measurement point..... | ||||
C | ||||
IF (QSPHERICAL .NE. 0) THEN | ||||
SLATGD = SLATGC | ||||
GDLAT = GCLAT | ||||
GDALT = RKM - SR | ||||
END IF | ||||
C | ||||
IF (QDIAG .NE. 0) THEN | ||||
WRITE (6,'(''Locations of radar and measurement:'')') | ||||
WRITE (6,FMT='(''slatgd,slon,sr,slatgc,tm,az,el,range = '')') | ||||
WRITE (6,FMT='(8f9.2)') SLATGD,SLON,SR,SLATGC,TM,AZ,EL,RANGE | ||||
WRITE (6,FMT='(''gclat,glon,rkm = '')') | ||||
WRITE (6,FMT='(3f9.2)') GCLAT,GLON,RKM | ||||
WRITE (6,FMT='(''gdlat,glon,gdalt = '')') | ||||
WRITE (6,FMT='(3f9.2)') GDLAT,GLON,GDALT | ||||
WRITE (6,FMT='('' '')') | ||||
END IF | ||||
C | ||||
C .....Calculate magnetic field at observation point..... | ||||
T = DTR*(90.D0-GCLAT) | ||||
CT = DCOS(T) | ||||
ST = DSIN(T) | ||||
P = DTR*GLON | ||||
CP = DCOS(P) | ||||
SP = DSIN(P) | ||||
CALL MILMAG(TM1,RKM,ST,CT,SP,CP,BR,BT,BP,B) | ||||
IF (QSPHERICAL .NE. 0) THEN | ||||
BR = -1.0D0 | ||||
BT = 0.0D0 | ||||
BP = 0.0D0 | ||||
END IF | ||||
CALL GDV(GDLAT,GCLAT,BR,BT,BP,XB,YB,ZB) | ||||
C | ||||
C | ||||
C .....Fill rcor with calculated coordinates..... | ||||
GDLAT = GDLATS | ||||
RCOR(1) = AZ | ||||
RCOR(2) = EL | ||||
RCOR(3) = RANGE | ||||
RCOR(4) = GDLAT | ||||
RCOR(5) = GLON | ||||
RCOR(6) = GDALT | ||||
RCOR(7) = B | ||||
RCOR(8) = BR | ||||
RCOR(9) = BT | ||||
RCOR(10) = BP | ||||
RCOR(11) = 0.0 | ||||
RCOR(12) = 0.0 | ||||
RCOR(13) = 0.0 | ||||
RCOR(14) = 0.0 | ||||
RCOR(15) = 0.0 | ||||
RCOR(16) = 0.0 | ||||
RCOR(17) = 0.0 | ||||
RCOR(18) = 0.0 | ||||
RCOR(19) = 0.0 | ||||
RCOR(20) = 0.0 | ||||
RCOR(21) = 0.0 | ||||
RCOR(22) = 0.0 | ||||
RCOR(23) = 0.0 | ||||
RCOR(24) = 0.0 | ||||
RCOR(25) = 0.0 | ||||
RCOR(26) = 0.0 | ||||
RCOR(27) = 0.0 | ||||
RCOR(28) = 0.0 | ||||
RCOR(29) = 0.0 | ||||
RCOR(30) = 0.0 | ||||
RCOR(31) = 0.0 | ||||
RCOR(32) = 0.0 | ||||
RCOR(33) = 0.0 | ||||
RCOR(34) = 0.0 | ||||
RCOR(35) = 0.0 | ||||
RCOR(36) = 0.0 | ||||
RCOR(37) = 0.0 | ||||
RCOR(38) = 0.0 | ||||
RCOR(39) = 0.0 | ||||
RCOR(40) = 0.0 | ||||
RCOR(41) = 0.0 | ||||
RCOR(42) = 0.0 | ||||
RCOR(43) = 0.0 | ||||
RCOR(44) = 0.0 | ||||
RCOR(45) = 0.0 | ||||
RCOR(46) = 0.0 | ||||
RCOR(47) = 0.0 | ||||
RCOR(45) = 0.0 | ||||
RCOR(46) = 0.0 | ||||
RCOR(47) = 0.0 | ||||
RCOR(48) = 0.0 | ||||
RCOR(49) = 0.0 | ||||
RCOR(50) = 0.0 | ||||
RCOR(51) = 0.0 | ||||
RCOR(52) = 0.0 | ||||
RCOR(53) = 0.0 | ||||
RCOR(54) = 0.0 | ||||
C | ||||
RETURN | ||||
END | ||||