C $Id: coordf.f 6710 2019-01-24 14:12:14Z brideout $ C SUBROUTINE COORDF(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 a faster version of coord.f used by getMag and faraday C because coord is slower calculating parms not needed for C those methods. RCOR(36-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 IF ((QDIAG .NE. 0) ) THEN WRITE (6,FMT='(''b = '', 9F9.5)') BR,BT,BP,XB,YB,ZB,BX,BY,BZ END IF C C .....Calculate inclination..... HB = DSQRT(XB*XB+YB*YB) AINC = DATAN2(ZB,HB)/DTR C C .....Calculate declination..... DEC = DATAN2(YB,XB)/DTR C C .....Convert southward, eastward, outward components of magnetic C field at observation point to earth centered cartesian C coordinates..... CALL VCTCNV(BX,BY,BZ,PX,PY,PZ,BR,BT,BP,RKM,90.D0-GCLAT,GLON,2) C C .....Compute unit vectors perpendicular to magnetic field and C southward (bperpsouth) and eastward (bperpeast) DECGC = DATAN2(BP,-BT)/DTR IF (QSPHERICAL .NE. 0) THEN AINC = 90.0D0 DECGC = 0.0D0 END IF IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''ainc,dec,decgc = '', 3f10.5)') AINC,DEC,DECGC END IF SOUTHR = 0.0D0 SOUTHT = COS(DTR*DECGC) SOUTHP = -SIN(DTR*DECGC) EASTR = 0.0D0 EASTT = SIN(DTR*DECGC) EASTP = COS(DTR*DECGC) UPR = 1.0D0 UPT = 0.0D0 UPP = 0.0D0 CALL VCTCNV(SOUTHX,SOUTHY,SOUTHZ,PX,PY,PZ,SOUTHR,SOUTHT,SOUTHP, * RKM,90.D0-GCLAT,GLON,2) CALL VCTCNV(EASTX,EASTY,EASTZ,PX,PY,PZ,EASTR,EASTT,EASTP,RKM, * 90.D0-GCLAT,GLON,2) CALL VCTCNV(UPX,UPY,UPZ,PX,PY,PZ,UPR,UPT,UPP,RKM,90.D0-GCLAT,GLON, * 2) SOUTHM = VMAG(SOUTH) SOUTHX = SOUTHX/SOUTHM SOUTHY = SOUTHY/SOUTHM SOUTHZ = SOUTHZ/SOUTHM EASTM = VMAG(EAST) EASTX = EASTX/EASTM EASTY = EASTY/EASTM EASTZ = EASTZ/EASTM CALL VPROD(BF,EAST,BPERPSOUTH) BPERPSOUTHM = VMAG(BPERPSOUTH) PERPSOUTH(1) = BPERPSOUTH(1)/BPERPSOUTHM PERPSOUTH(2) = BPERPSOUTH(2)/BPERPSOUTHM PERPSOUTH(3) = BPERPSOUTH(3)/BPERPSOUTHM CALL VPROD(SOUTH,BF,BPERPEAST) BPERPEASTM = VMAG(BPERPEAST) PERPEAST(1) = BPERPEAST(1)/BPERPEASTM PERPEAST(2) = BPERPEAST(2)/BPERPEASTM PERPEAST(3) = BPERPEAST(3)/BPERPEASTM UPB(1) = -BF(1)/VMAG(BF) UPB(2) = -BF(2)/VMAG(BF) UPB(3) = -BF(3)/VMAG(BF) C IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''southr,southt,southp,eastr,eastt,eastp = '')') WRITE (6,FMT='(3e13.5,3x,3e13.5)') SOUTHR,SOUTHT,SOUTHP,EASTR, * EASTT,EASTP WRITE (6,FMT='(''south = '', 4e13.5)') SOUTH,SPROD(SOUTH,SOUTH) WRITE (6,FMT='('' east = '', 4e13.5)') EAST,SPROD(EAST,EAST) WRITE (6,FMT='('' up = '', 4e13.5)') UP,SPROD(UP,UP) SP = SPROD(SOUTH,EAST) WRITE (6,FMT='(''sp = '', e13.5)') SP CALL VPROD(SOUTH,EAST,ZEROM) WRITE (6,FMT='(''zerom = '', 4e13.5)') ZEROM,SPROD(ZEROM,ZEROM) CALL VPROD(UP,SOUTH,EAST2) CALL VPROD(EAST,UP,SOUTH2) CALL VPROD(SOUTH,EAST,UP2) WRITE (6,FMT='(''south = '', 3e13.5)') SOUTH WRITE (6,FMT='(''south2 = '', 3e13.5)') SOUTH2 WRITE (6,FMT='(''east = '', 3e13.5)') EAST WRITE (6,FMT='(''east2 = '', 3e13.5)') EAST2 WRITE (6,FMT='(''up = '', 3e13.5)') UP WRITE (6,FMT='(''up2 = '', 3e13.5)') UP2 WRITE (6,FMT='(''bperpsouth,bperpsouthm,perpsouth = '')') WRITE (6,FMT='(7e13.5)') BPERPSOUTH,BPERPSOUTHM,PERPSOUTH WRITE (6,FMT='(''bperpeast,bperpeastm,perpeast = '')') WRITE (6,FMT='(7e13.5)') BPERPEAST,BPERPEASTM,PERPEAST WRITE (6,FMT='(''bf,bfm,upb = '', 7e13.5)') BF,VMAG(BF),UPB WRITE (6,FMT='(''perpsouth,perpeast,upb = '')') WRITE (6,FMT='(2f8.2,3f8.4,2x,3f8.4,2x,3f8.4)') GCLAT,GLON, * PERPSOUTH,PERPEAST,UPB CALL VCTCNV(PERPSOUTH(1),PERPSOUTH(2),PERPSOUTH(3),PX,PY,PZ, * PERPSOUTHR,PERPSOUTHT,PERPSOUTHP,RKM,90.0D0-GCLAT, * GLON,1) CALL VCTCNV(PERPEAST(1),PERPEAST(2),PERPEAST(3),PX,PY,PZ, * PERPEASTR,PERPEASTT,PERPEASTP,RKM,90.0D0-GCLAT, * GLON,1) CALL VCTCNV(UPB(1),UPB(2),UPB(3),PX,PY,PZ,UPBR,UPBT,UPBP,RKM, * 90.0D0-GCLAT,GLON,1) WRITE (6,FMT='(3f8.4,2x,3f8.4,2x,3f8.4)') PERPSOUTHR, * PERPSOUTHT,PERPSOUTHP WRITE (6,FMT='(3f8.4,2x,3f8.4,2x,3f8.4)') PERPEASTR,PERPEASTT, * PERPEASTP WRITE (6,FMT='(3f8.4,2x,3f8.4,2x,3f8.4)') UPBR,UPBT,UPBP CALL VPROD(UPB,PERPSOUTH,PERPEAST2) WRITE (6,FMT='(''perpeast = '', 3e13.5)') PERPEAST WRITE (6,FMT='(''perpeast2 = '', 3e13.5)') PERPEAST2 CALL VPROD(PERPEAST,UPB,PERPSOUTH2) WRITE (6,FMT='(''perpsouth = '', 3e13.5)') PERPSOUTH WRITE (6,FMT='(''perpsouth2 = '', 3e13.5)') PERPSOUTH2 CALL VPROD(PERPSOUTH,PERPEAST,UPB2) WRITE (6,FMT='(''upb = '', 3e13.5)') UPB WRITE (6,FMT='(''upb2 = '', 3e13.5)') UPB2 END IF C C .....Calculate l-shell parameter..... GCALT = RKM - 6378.16D0/DSQRT(1.D0+.0067397D0*ST*ST) CALL INVAR(TM1,GCLAT,GLON,GCALT,0.01D0,BB,RL) RL = MAX(RL,1.0D0) C C .....Calculate dip latitude..... RLATM = ATAN(.5D0*TAN(DTR*AINC))/DTR C C .....Calculate invariant latitude..... RLATI = ACOS(DSQRT(1.0D0/RL))/DTR C C .....Convert radar propagation vector and observation point C position to earth centered cartesian coordinates..... CALL RPCART(SR,SLATGC,SLON,AZ,EL,RANGE,RFX,RFY,RFZ,PFX,PFY,PFZ) C C .....Calculate earth-centered cartestian components of unit vector C perpendicular to radar propagation vector and geomagnetic C geomagnetic field at observation point..... CALL VPROD(RF,BF,EF) EFMAG = VMAG(EF,EF) EFX = EFX/EFMAG EFY = EFY/EFMAG EFZ = EFZ/EFMAG C C .....Calculate aspect angle..... CASPCT = SPROD(RF,BF)/(VMAG(RF)*VMAG(BF)) CASPCT = MIN(MAX(CASPCT,-1.0D0),1.0D0) ASPCT = ACOS(CASPCT)/DTR C C .....Calculate direction cosines of radar beam with respect to C geocentric south, east, up..... CALL VCTCNV(RFX,RFY,RFZ,PFX,PFY,PFZ,RFR,RFT,RFP,RR,RT,RP,1) CST = RFT/RANGE CSP = RFP/RANGE CSR = RFR/RANGE C C .....Compute geodetic direction cosines. the direction cosines C are with respect to south, east, up rather than x (north), C y (east), z (down) as used in gdv..... CALL GDV(GDLAT,GCLAT,CSR,CST,CSP,CX,CY,CZ) CTAB(1) = -CX CTAB(2) = CY CTAB(3) = -CZ C IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''rl,rlatm,rlati,aspct = '', 4f9.2)') RL,RLATM, * RLATI,ASPCT WRITE (6,FMT='(''rfx,rfy,rfz = '', 3f10.2)') RFX,RFY,RFZ WRITE (6,FMT='(''cst,csp,csr = '', 3f9.5)') CST,CSP,CSR WRITE (6,FMT='(''ctab = '', 3f9.5)') CTAB END IF C C .....Calculate observation point apex coordinates..... ISTOP = 0 NPR = 0 INIT = 0 D = 1.0D-6 CALL LINTRA(TM1,GCLAT,GLON,RKM,GDALT,0.0D0,PLAT,PLON,PRKM,ARC, * ARAD,ALAT,ALON,ISTOP,NPR,INIT,IER) Q0(1) = DTR*(90.0D0-ALAT) Q0(2) = DTR*ALON C C .....Calculate arc length along field line from surface to C observation point..... ISTOP = -1 CALL LINTRA(TM1,GCLAT,GLON,RKM,GDALT,0.0D0,PLAT,PLON,PRKM,ARC, * ARAD,ALAT,ALON,ISTOP,NPR,INIT,IER) Q0(3) = ARC/6370.0D0 C C .....Get conjugate point using lintra..... HI = DBLE(GDALT) C C .....Tell lintra to find mag conj, not apex..... ISTOP = 1 CALL LINTRA(TM,GCLAT,GLON,RKM,HI,GDALT,CPLAT,CPLON,CPRKM,CARC, * CARAD,CALAT,CALON,ISTOP,NPR,INIT,IER) C .....Get geodetic latitude of conjugate point..... CALL CONVRT(2,CGDLAT,CGDALT,CPLAT,CPRKM) IF (CPLON.LT.0.0D0) THEN CPLON = CPLON + 360.0D0 END IF 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) = RLATM RCOR(12) = RLATI RCOR(13) = RL RCOR(14) = ALAT RCOR(15) = ALON c RCOR(16) = GF(1) c RCOR(17) = GF(2) c RCOR(18) = GF(3) c RCOR(19) = GF(4) RCOR(16) = 0.0 RCOR(17) = 0.0 RCOR(18) = 0.0 RCOR(19) = 0.0 RCOR(20) = CTAB(1) RCOR(21) = CTAB(2) RCOR(22) = CTAB(3) 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) = AINC RCOR(30) = DEC RCOR(31) = GCLAT RCOR(32) = ASPCT RCOR(33) = CPLAT RCOR(34) = CGDLAT RCOR(35) = CPLON 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