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