C $Id: coord.f 5103 2015-06-11 12:38:33Z brideout $ C SUBROUTINE COORD(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 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 The grad* parameters are gradients of the geomagnetic field B. C They measure the separation between field lines and hence the C variation of the electric field as a function of position. C this permits us to represent the perpendicular E field and C hence the E X B drift along a field line by a single C parameter. Note that only the relative value of these C parameters is significant. C along the field line. C RCOR(36) - gradapex1(1) - apex latitude component of gradient C of B C RCOR(37) - gradapex2(2) - apex lonfitude component gradient C of B C RCOR(38) - gradapex3(3) - up b component of gradient of B C gradmag = gradient of B converted to standard magnetic south, C east, up C RCOR(39) - gradmag1(1) - gradapex1 magnetic south component C RCOR(40) - gradmag1(2) - gradapex1 magnetic east component C RCOR(41) - gradmag1(3) - gradapex1 up B component C RCOR(42) - gradmag2(1) - gradapex2 magnetic south component C RCOR(43) - gradmag2(2) - gradapex2 magnetic east component C RCOR(44) - gradmag3(3) - gradapex2 up B component C RCOR(45) - gradmag3(1) - gradapex3 magnetic south component C RCOR(46) - gradmag3(2) - gradapex3 magnetic east component C RCOR(47) - gradmag3(3) - gradapex3 up B component C The following three parameters are the direction cosines of C the radar beam wrt standard magnetic south, magnetic east C and up B directions. These are useful for expressing final C results but not for calculating the variation of E along a C field line because they do not define a coordinate system C (in which the gradient of a potential can be calculated) but C only a set of orthogonal directions. C RCOR(48) - ctab(7) - direction cosine of beam wrt vector C perpendicular to B in magnetic meridian C RCOR(49) - ctab(8) - direction cosine of beam wrt vector C perpendicular to ctab(7) and the magnetic C meridian C RCOR(50) - ctab(9) - direction cosine of beam wrt up B C direction cosines of measured component of electric field wrt C apex contravariant unit vectors C RCOR(51) - efcoscn1 - wrt apex latitude C RCOR(52) - efcoscn1 - wrt apex longitude C RCOR(53) - efcoscn1 - wrt up B C RCOR(54) = ARC (?) C RCOR(55-64) - undefined 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,COSCN1,COSCN2,COSCN3, * COSCO1,COSCO2,COSCO3,COST1,COST2,COST3,CP,CPLAT, * CPLON,CPRKM,CSP,CSR,CST,CT,CTABM,CX,CY,CZ,D,DEC, * DECGC,DET,DGCLAT,DGDALT,DGDLAT,DGLON,DRKM,DTR, * EAARC,EALAT,EALON,EASTM,EASTP,EASTR,EASTT,EASTX, * EASTY,EASTZ,ECN1M,ECN2M,ECN3M,ECO1M,ECO2M,ECO3M, * EFCOSCN1,EFCOSCN2,EFCOSCN3,EFMAG,EFX,EFY,EFZ, * GCALT,GCLAT,GDLATS,GRADAPEXM,GRADMAGM,GRADXYZM, * HB,HI,P,PERPEASTP,PERPEASTR,PERPEASTT,PERPSOUTHP, * PERPSOUTHR,PERPSOUTHT,PFX,PFY,PFZ,PHI,PLAT,PLON, * PRKM,PX,PY,PZ,Q1,Q2,Q3,R,RFP,RFR,RFT,RFX,RFY,RFZ, * RKM,RL,RLATI,RLATM,RMAG,RP,RR,RT,SOUTHM,SOUTHP, * SOUTHR,SOUTHT,SOUTHX,SOUTHY,SOUTHZ,SP,ST,T,THETA, * TM1,TP,TR,TT,UPBP,UPBR,UPBT,UPP,UPR,UPT,UPX,UPY, * UPZ,VCN,VCO,X1,X2,X3,XB,XDUM,YB,ZB INTEGER I,IER,INIT,ISTOP,J,K,NPR C .. C .. Local Arrays .. DOUBLE PRECISION BF(3),BPERPEAST(3),BPERPSOUTH(3),CTAB(9), * DQDX(3,3),DQDX1(3,3),DXDQ(3,3),DXDQSP(3,3), * EAST(3),EAST2(3),ECN1(3),ECN1X(3),ECN2(3), * ECN2X(3),ECN3(3),ECN3X(3),ECNU1(3),ECNU2(3), * ECNU3(3),ECO1(3),ECO1X(3),ECO2(3),ECO2X(3), * ECO3(3),ECO3X(3),EF(3),GCN(3,3),GCNSPHERE(3,3), * GCO(3,3),GCO1(3,3),GF(4),GRADAPEX(3), * GRADAPEXSP1(3),GRADAPEXSP2(3),GRADAPEXSP3(3), * GRADMAG(3),GRADXYZ(3),PERPEAST(3),PERPEAST2(3), * PERPSOUTH(3),PERPSOUTH2(3),PF(3),Q(3),Q0(3), * RF(3),SOUTH(3),SOUTH2(3),SV(3),TV(3),UP(3), * UP2(3),UPB(3),UPB2(3),X(3),X0(3),ZEROM(3) INTEGER LW(3),MW(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 IF (QSPHERICAL .NE. 0) THEN Q0(1) = DTR*(90.0D0-GCLAT) Q0(2) = DTR*GLON Q0(3) = RKM/6370.0D0 THETA = Q0(1) PHI = Q0(2) R = Q0(3) END IF C C .....Calculate apex-cartesian direction cosines..... X0(1) = PX/6370.0D0 X0(2) = PY/6370.0D0 X0(3) = PZ/6370.0D0 X(1) = X0(1) X(2) = X0(2) X(3) = X0(3) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(3f7.1,2x,2f10.5,f9.3,2x,3f10.6,2x,3f10.6)') AZ, * EL,RANGE,ALAT,ALON,ARC,Q0(1),Q0(2),Q0(3),X(1),X(2),X(3) END IF DO 20 J = 1,3 X(J) = X(J) + D X1 = X(1) X2 = X(2) X3 = X(3) CALL CSCONV(6370.0D0*X1,6370.0D0*X2,6370.0D0*X3,DRKM,DGCLAT, * DGLON,1) DGCLAT = 90.0D0 - DGCLAT CALL CONVRT(2,DGDLAT,DGDALT,DGCLAT,DRKM) ISTOP = 0 CALL LINTRA(TM1,DGCLAT,DGLON,DRKM,DGDALT,0.0D0,PLAT,PLON,PRKM, * ARC,ARAD,ALAT,ALON,ISTOP,NPR,INIT,IER) ISTOP = -1 CALL LINTRA(TM1,DGCLAT,DGLON,DRKM,GDALT,0.0D0,PLAT,PLON,PRKM, * ARC,ARAD,ALAT,ALON,ISTOP,NPR,INIT,IER) Q(1) = DTR*(90.0D0-ALAT) Q(2) = DTR*ALON Q(3) = ARC/6370.0D0 C IF (QSPHERICAL .NE. 0) THEN Q(1) = DTR*(90.0D0-DGCLAT) Q(2) = DTR*DGLON Q(3) = DRKM/6370.0D0 END IF IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(3f7.1,2x,2f10.5,f9.3,2x,3f10.6,2x,3f10.6)') * AZ,EL,RANGE,ALAT,ALON,ARC,Q(1),Q(2),Q(3),X(1),X(2),X(3) END IF X(J) = X0(J) DO 10 I = 1,3 DQDX(I,J) = (Q(I)-Q0(I))/D 10 CONTINUE 20 CONTINUE IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''dqdx (finite differences) = '')') DO 30 I = 1,3 WRITE (6,FMT='(3E13.5)') (DQDX(I,J),J=1,3) 30 CONTINUE END IF C DO 50 I = 1,3 DO 40 J = 1,3 DXDQ(I,J) = DQDX(I,J) 40 CONTINUE 50 CONTINUE CALL MINV(DXDQ,3,DET,LW,MW) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''dxdq (minv) = '')') DO 60 I = 1,3 WRITE (6,FMT='(3E13.5)') (DXDQ(I,J),J=1,3) 60 CONTINUE WRITE (6,FMT='(''det = '', E13.5)') DET END IF C IF (QSPHERICAL .NE. 0) THEN C x = r*sin(theta)*cos(phi) C y = r*sin(theta)*sin(phi) C z = r*cos(theta) DXDQSP(1,1) = R*COS(THETA)*COS(PHI) DXDQSP(2,1) = R*COS(THETA)*SIN(PHI) DXDQSP(3,1) = -R*SIN(THETA) DXDQSP(1,2) = -R*SIN(THETA)*SIN(PHI) DXDQSP(2,2) = R*SIN(THETA)*COS(PHI) DXDQSP(3,2) = 0.0D0 DXDQSP(1,3) = SIN(THETA)*COS(PHI) DXDQSP(2,3) = SIN(THETA)*SIN(PHI) DXDQSP(3,3) = COS(THETA) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''dxdq (spherical coordinates) = '')') DO 70 I = 1,3 WRITE (6,FMT='(3E13.5)') (DXDQSP(I,J),J=1,3) 70 CONTINUE WRITE (6,FMT='(''dxdqsp should be almost equal to dxdq'')') K = 0 DO 90 I = 1,3 DO 80 J = 1,3 IF (ABS(DXDQ(I,J)-DXDQSP(I,J)).GT.1.0D-05) K = K + 1 80 CONTINUE 90 CONTINUE IF (K.EQ.0) THEN WRITE (6,FMT='(''Test result: OK'')') ELSE WRITE (6,FMT='(''Test result: FAILED'')') END IF WRITE (6,FMT='('' '')') END IF END IF C DO 110 I = 1,3 DO 100 J = 1,3 DQDX1(I,J) = DXDQ(I,J) 100 CONTINUE 110 CONTINUE CALL MINV(DQDX1,3,DET,LW,MW) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''dqdx1 (minv) = '')') DO 120 I = 1,3 WRITE (6,FMT='(3E13.5)') (DQDX1(I,J),J=1,3) 120 CONTINUE WRITE (6,FMT='(''det = '', E13.5)') DET END IF C C .....Construct base vectors (= rows dqdx, dxdq)..... DO 130 I = 1,3 ECN1(I) = DQDX(1,I) ECN2(I) = DQDX(2,I) ECN3(I) = DQDX(3,I) 130 CONTINUE IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''ecn1 = '', 3E13.5)') (ECN1(J),J=1,3) WRITE (6,FMT='(''ecn2 = '', 3E13.5)') (ECN2(J),J=1,3) WRITE (6,FMT='(''ecn3 = '', 3E13.5)') (ECN3(J),J=1,3) END IF C DO 140 I = 1,3 ECO1(I) = DXDQ(I,1) ECO2(I) = DXDQ(I,2) ECO3(I) = DXDQ(I,3) 140 CONTINUE IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''eco1 = '', 3E13.5)') (ECO1(J),J=1,3) WRITE (6,FMT='(''eco2 = '', 3E13.5)') (ECO2(J),J=1,3) WRITE (6,FMT='(''eco3 = '', 3E13.5)') (ECO3(J),J=1,3) END IF C C .....Compute covariant base vectors from contravariant C base vectors..... VCN = TPROD(ECN1,ECN2,ECN3) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''vcn = '', E13.5)') VCN WRITE (6,FMT='('' '')') CALL VRECIP(ECN1,ECN2,ECN3,ECO1X,ECO2X,ECO3X) WRITE (6,FMT='(''eco1x = '', 3E13.5)') (ECO1X(J),J=1,3) WRITE (6,FMT='(''eco2x = '', 3E13.5)') (ECO2X(J),J=1,3) WRITE (6,FMT='(''eco3x = '', 3E13.5)') (ECO3X(J),J=1,3) END IF C VCO = TPROD(ECO1,ECO2,ECO3) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''vco = '', E13.5)') VCO WRITE (6,FMT='('' '')') END IF CALL VRECIP(ECO1,ECO2,ECO3,ECN1X,ECN2X,ECN3X) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''ecn1x = '', 3E13.5)') (ECN1X(J),J=1,3) WRITE (6,FMT='(''ecn2x = '', 3E13.5)') (ECN2X(J),J=1,3) WRITE (6,FMT='(''ecn3x = '', 3E13.5)') (ECN3X(J),J=1,3) C WRITE (6,FMT='('' '')') SP = SPROD(ECN1,ECO1) WRITE (6,FMT='(''sp = '', E13.5)') SP SP = SPROD(ECN1,ECO2) WRITE (6,FMT='(''sp = '', E13.5)') SP SP = SPROD(ECN1,ECO3) WRITE (6,FMT='(''sp = '', E13.5)') SP SP = SPROD(ECN2,ECO1) WRITE (6,FMT='(''sp = '', E13.5)') SP SP = SPROD(ECN2,ECO2) WRITE (6,FMT='(''sp = '', E13.5)') SP SP = SPROD(ECN2,ECO3) WRITE (6,FMT='(''sp = '', E13.5)') SP SP = SPROD(ECN3,ECO1) WRITE (6,FMT='(''sp = '', E13.5)') SP SP = SPROD(ECN3,ECO2) WRITE (6,FMT='(''sp = '', E13.5)') SP SP = SPROD(ECN3,ECO3) WRITE (6,FMT='(''sp = '', E13.5)') SP END IF C C .....Calculate covariant apex metric tensor from C derivative matrix..... CALL GMET(DXDQ,GCO) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''gco [gmet(dxdq)] = '')') DO 150 I = 1,3 WRITE (6,FMT='(3E13.5)') (GCO(I,J),J=1,3) 150 CONTINUE END IF C C .....Calculate contravariant apex metric tensor from C derivative matrix..... CALL MTRAN3(DQDX) CALL GMET(DQDX,GCN) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''gcn [gmet(dqdx)] = '')') DO 160 I = 1,3 WRITE (6,FMT='(3E13.5)') (GCN(I,J),J=1,3) 160 CONTINUE END IF C C .....Calculate covariant apex metric tensor from C base vectors..... CALL GMETBV(ECO1,ECO2,ECO3,GCO) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''gco [gmetbv(eco)] = '')') DO 170 I = 1,3 WRITE (6,FMT='(3E13.5)') (GCO(I,J),J=1,3) 170 CONTINUE END IF C C .....Calculate contravariant apex metric tensor from C base vectors..... CALL GMETBV(ECN1,ECN2,ECN3,GCN) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''gcn [gmetbv(ecn)] = '')') DO 180 I = 1,3 WRITE (6,FMT='(3E13.5)') (GCN(I,J),J=1,3) 180 CONTINUE END IF C C .....Calculate apex contravariant metric tensor..... DO 200 I = 1,3 DO 190 J = 1,3 GCN(I,J) = GCO(I,J) 190 CONTINUE 200 CONTINUE CALL MINV(GCN,3,DET,LW,MW) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''gcn [minv(gco)] = '')') DO 210 I = 1,3 WRITE (6,FMT='(3E13.5)') (GCN(I,J),J=1,3) 210 CONTINUE WRITE (6,FMT='(''det = '', e13.5)') DET END IF C DO 230 I = 1,3 DO 220 J = 1,3 GCO1(I,J) = GCN(I,J) 220 CONTINUE 230 CONTINUE CALL MINV(GCO1,3,DET,LW,MW) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''gco1 [minv(gcn)] = '')') DO 240 I = 1,3 WRITE (6,FMT='(3E13.5)') (GCO1(I,J),J=1,3) 240 CONTINUE WRITE (52,FMT='(f8.1,4e13.5)') GDALT,GCN(1,1),GCN(1,2), * GCN(2,1),GCN(2,2) WRITE (6,FMT='(''det = '', e13.5)') DET WRITE (6,FMT='('' '')') END IF C C .....Compute metric tensor for spherical coordinates..... IF (QSPHERICAL .NE. 0) THEN Q1 = Q0(3) Q2 = Q0(3)*SIN(Q0(1)) Q3 = 1.0d0 DO 260 I = 1,3 DO 250 J = 1,3 GCNSPHERE(I,J) = 0.0D0 250 CONTINUE 260 CONTINUE GCNSPHERE(1,1) = 1.0d0/Q1**2 GCNSPHERE(2,2) = 1.0d0/Q2**2 GCNSPHERE(3,3) = 1.0d0/Q3**2 IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''Q**2 = '', 3e13.5)') Q1**2,Q2**2,Q3**2 WRITE (6,FMT='(''1/Q**2 = '', 3e13.5)') 1.0d0/Q1**2, * 1.0d0/Q2**2,1.0d0/Q3**2 WRITE (6,FMT='(''gcn [spherical coordinates] = '')') DO 270 I = 1,3 WRITE (6,FMT='(3E13.5)') (GCNSPHERE(I,J),J=1,3) 270 CONTINUE WRITE (6,FMT='(''Gcnsphere should be almost equal to gcn'')' * ) K = 0 DO 290 I = 1,3 DO 280 J = 1,3 IF (ABS(GCN(I,J)-GCNSPHERE(I,J)).GT.1.0D-05) K = K + 1 280 CONTINUE 290 CONTINUE IF (K.EQ.0) THEN WRITE (6,FMT='(''Test result: OK'')') ELSE WRITE (6,FMT='(''Test result: FAILED'')') END IF WRITE (6,FMT='('' '')') END IF END IF C C .....Compute direction cosines..... ECO1M = VMAG(ECO1) ECO2M = VMAG(ECO2) ECO3M = VMAG(ECO3) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''eco1m,eco2m,eco3m = '', 3e13.5)') ECO1M,ECO2M, * ECO3M ECO1M = VMAG(DXDQ(1,1)) ECO2M = VMAG(DXDQ(1,2)) ECO3M = VMAG(DXDQ(1,3)) WRITE (6,FMT='(''eco1m,eco2m,eco3m = '', 3e13.5)') ECO1M,ECO2M, * ECO3M END IF COSCO1 = SPROD(ECO1,RF) COSCO2 = SPROD(ECO2,RF) COSCO3 = SPROD(ECO3,RF) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''cosco1,cosco2,cosco3 = '', 3e13.5)') COSCO1, * COSCO2,COSCO3 COSCO1 = SPROD(DXDQ(1,1),RF) COSCO2 = SPROD(DXDQ(1,2),RF) COSCO3 = SPROD(DXDQ(1,3),RF) WRITE (6,FMT='(''cosco1,cosco2,cosco3 = '', 3e13.5)') COSCO1, * COSCO2,COSCO3 END IF C ECN1M = VMAG(ECN1) ECN2M = VMAG(ECN2) ECN3M = VMAG(ECN3) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''ecn1m,ecn2m,ecn3m = '', 3e13.5)') ECN1M,ECN2M, * ECN3M ECN1M = VMAG(DQDX(1,1)) ECN2M = VMAG(DQDX(1,2)) ECN3M = VMAG(DQDX(1,3)) WRITE (6,FMT='(''ecn1m,ecn2m,ecn3m = '', 3e13.5)') ECN1M,ECN2M, * ECN3M END IF CALL UVECT(ECN1,ECNU1) CALL UVECT(ECN2,ECNU2) CALL UVECT(ECN3,ECNU3) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''sprod(ecnu1,ecnu1) = '', e13.5)') SPROD(ECNU1, * ECNU1) WRITE (6,FMT='(''sprod(ecnu1,ecnu2) = '', e13.5)') SPROD(ECNU1, * ECNU2) WRITE (6,FMT='(''sprod(ecnu1,ecnu3) = '', e13.5)') SPROD(ECNU1, * ECNU3) WRITE (6,FMT='(''sprod(ecnu2,ecnu1) = '', e13.5)') SPROD(ECNU2, * ECNU1) WRITE (6,FMT='(''sprod(ecnu2,ecnu2) = '', e13.5)') SPROD(ECNU2, * ECNU2) WRITE (6,FMT='(''sprod(ecnu2,ecnu3) = '', e13.5)') SPROD(ECNU2, * ECNU3) WRITE (6,FMT='(''sprod(ecnu3,ecnu1) = '', e13.5)') SPROD(ECNU3, * ECNU1) WRITE (6,FMT='(''sprod(ecnu3,ecnu2) = '', e13.5)') SPROD(ECNU3, * ECNU2) WRITE (6,FMT='(''sprod(ecnu3,ecnu3) = '', e13.5)') SPROD(ECNU3, * ECNU3) WRITE (6,FMT='(''sprod1 = '', e13.5)') SPROD(ECN1,BF) WRITE (6,FMT='(''sprod2 = '', e13.5)') SPROD(ECN2,BF) WRITE (6,FMT='(''sprod3 = '', e13.5)') SPROD(ECN3,BF) WRITE (6,FMT='(''sprod1 = '', e13.5)') SPROD(ECO1,BF) WRITE (6,FMT='(''sprod2 = '', e13.5)') SPROD(ECO2,BF) WRITE (6,FMT='(''sprod3 = '', e13.5)') SPROD(ECO3,BF) END IF COSCN1 = SPROD(ECN1,RF) COSCN2 = SPROD(ECN2,RF) COSCN3 = SPROD(ECN3,RF) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''coscn1,coscn2,coscn3 = '', 3e13.5)') COSCN1, * COSCN2,COSCN3 COSCN1 = SPROD(DQDX(1,1),RF) COSCN2 = SPROD(DQDX(1,2),RF) COSCN3 = SPROD(DQDX(1,3),RF) WRITE (6,FMT='(''coscn1,coscn2,coscn3 = '', 3e13.5)') COSCN1, * COSCN2,COSCN3 END IF C RMAG = VMAG(RF) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''rmag = '', e13.5)') RMAG END IF C CTAB(4) = COSCN1/(ECN1M*RMAG) CTAB(5) = COSCN2/(ECN2M*RMAG) CTAB(6) = COSCN3/(ECN3M*RMAG) CTAB(4) = SPROD(ECN1,RF)/(VMAG(ECN1)*VMAG(RF)) CTAB(5) = SPROD(ECN2,RF)/(VMAG(ECN2)*VMAG(RF)) CTAB(6) = SPROD(BF,RF)/(VMAG(BF)*VMAG(RF)) CTABM = SQRT(CTAB(4)**2+CTAB(5)**2+CTAB(6)**2) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='('' '')') WRITE (6,FMT='(''ctab(4),ctab(5),ctab(6),ctabm = '', 4e13.5)') * CTAB(4),CTAB(5),CTAB(6),CTABM WRITE (6,FMT='('' '')') END IF C C Compute direction cosines with respect to contravariant unit C vectors of measured component of electric field. EFCOSCN1 = SPROD(ECN1,EF)/(VMAG(ECN1)*VMAG(EF)) EFCOSCN2 = SPROD(ECN2,EF)/(VMAG(ECN2)*VMAG(EF)) EFCOSCN3 = SPROD(BF,EF)/(VMAG(BF)*VMAG(EF)) C GF(1) = GCN(1,1)*ECO1M GF(2) = GCN(2,1)*ECO1M GF(3) = GCN(1,2)*ECO2M GF(4) = GCN(2,2)*ECO2M C C .....Here we assert that the field line is an equipotential..... EALAT = 1.0D0 EALON = 1.0D0 EAARC = 1.0D0 C IF (QSPHERICAL .NE. 0) THEN GRADAPEXSP1(1) = 1/Q1 GRADAPEXSP1(2) = 0.0D0 GRADAPEXSP1(3) = 0.0D0 GRADAPEXSP2(1) = 0.0D0 GRADAPEXSP2(2) = 1/Q2 GRADAPEXSP2(3) = 0.0D0 GRADAPEXSP3(1) = 0.0D0 GRADAPEXSP3(2) = 0.0D0 GRADAPEXSP3(3) = 0.0D0 GRADAPEXSP3(3) = 1/Q3 IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''gradapexsp = '')') WRITE (6,FMT='(4f9.5,2x,4f9.5)') GRADAPEXSP1, * VMAG(GRADAPEXSP1) WRITE (6,FMT='(4f9.5,2x,4f9.5)') GRADAPEXSP2, * VMAG(GRADAPEXSP2) WRITE (6,FMT='(4f9.5,2x,4f9.5/)') GRADAPEXSP3, * VMAG(GRADAPEXSP3) END IF END IF C GRADAPEX(1) = 1.0D0 GRADAPEX(2) = 1.0D0 GRADAPEX(3) = 0.0D0 GRADXYZ(1) = ECN1(1)*GRADAPEX(1) + ECN2(1)*GRADAPEX(2) + * ECN3(1)*GRADAPEX(3) GRADXYZ(2) = ECN1(2)*GRADAPEX(1) + ECN2(2)*GRADAPEX(2) + * ECN3(2)*GRADAPEX(3) GRADXYZ(3) = ECN1(3)*GRADAPEX(1) + ECN2(3)*GRADAPEX(2) + * ECN3(3)*GRADAPEX(3) GRADMAG(1) = SPROD(PERPSOUTH,GRADXYZ) GRADMAG(2) = SPROD(PERPEAST,GRADXYZ) GRADMAG(3) = SPROD(UPB,GRADXYZ) GRADAPEXM = VMAG(GRADAPEX) GRADXYZM = VMAG(GRADXYZ) GRADMAGM = VMAG(GRADMAG) IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''gclat,glon,sr,gdalt = '', 4f10.2/)') GCLAT, * GLON,SR,GDALT WRITE (6,FMT='(''eco1,ecn1 = '', 3f9.5,2x,3f9.5,2x,f9.5)') * ECO1,ECN1,SPROD(ECO1,ECN1) WRITE (6,FMT='(''eco2,ecn2 = '', 3f9.5,2x,3f9.5,2x,f9.5)') * ECO2,ECN2,SPROD(ECO2,ECN2) WRITE (6,FMT='(''eco3,ecn3 = '', 3f9.5,2x,3f9.5,2x,f9.5/)') * ECO3,ECN3,SPROD(ECO3,ECN3) WRITE (6,FMT='(''gradapex,gradxyz = '')') WRITE (6,FMT='(4f9.5,2x,4f9.5)') GRADAPEX,VMAG(GRADAPEX), * GRADXYZ,VMAG(GRADXYZ) END IF IF ((QSPHERICAL .NE. 0) .AND. (QDIAG .NE. 0)) THEN WRITE (6,FMT='(''theta,phi = '', 2f9.5)') THETA/DTR,PHI/DTR WRITE (6,FMT='(''gradxyzsp = '')') WRITE (6,FMT='(3f9.5)') COS(THETA)*COS(PHI), * COS(THETA)*SIN(PHI),-SIN(THETA) WRITE (6,FMT='(3f9.5)') - SIN(PHI),COS(PHI),0.0D0 WRITE (6,FMT='(3f9.5/)') SIN(THETA)*COS(PHI), * SIN(THETA)*SIN(PHI),COS(THETA) END IF C IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''sprod11 = '', f10.5)') SPROD(GRADXYZ,GRADXYZ) WRITE (6,FMT='(3f10.5,2x,4f10.5,2x,4f10.5)') GRADAPEX,GRADXYZ, * VMAG(GRADXYZ),GRADMAG,VMAG(GRADMAG) WRITE (6,FMT='(''bf = '', 4f10.5)') BF,VMAG(BF) WRITE (6,FMT='(''sprod(gradxyz,bf) = '', f10.5)') * SPROD(GRADXYZ,BF) END IF IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''decgc,ainc = '', 2f10.5)') DECGC,AINC WRITE (6,FMT='(''perpsouth = '', 4f10.5)') PERPSOUTH, * VMAG(PERPSOUTH) WRITE (6,FMT='(''perpeast = '', 4f10.5)') PERPEAST, * VMAG(PERPEAST) WRITE (6,FMT='(''upb = '', 4f10.5)') UPB,VMAG(UPB) WRITE (6,FMT='(''gradmag1 = '', 3f10.5)') GRADMAG WRITE (6,FMT= *'(''Vectors 1, 2 and 3 must have the same length in '', * ''all three coordinate systems'')') WRITE (6,FMT='(3f9.5,2x,3f9.5,2x,3f9.5,2x,3f9.5)') * SPROD(GRADAPEX,GRADAPEX),SPROD(GRADXYZ,GRADXYZ) END IF CTAB(7) = SPROD(PERPSOUTH,RF) CTAB(8) = SPROD(PERPEAST,RF) CTAB(9) = SPROD(UPB,RF) CTABM = VMAG(CTAB(7)) CTAB(7) = CTAB(7)/CTABM CTAB(8) = CTAB(8)/CTABM CTAB(9) = CTAB(9)/CTABM C IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''CTAB(7),CTAB(8),CTAB(9) = '')') WRITE (6,FMT='(4E13.5)') CTAB(7),CTAB(8),CTAB(9),CTABM END IF C C .....Spherical coordinate test case..... IF (QSPHERICAL .NE. 0) THEN GF(1) = GCN(1,1)*ECO2M GF(2) = GCN(2,1)*ECO2M GF(3) = GCN(1,2)*ECO3M GF(4) = GCN(2,2)*ECO3M IF (QDIAG .NE. 0) THEN WRITE (6,FMT='(''gf = '')') WRITE (6,FMT='(2e13.5)') GF(1),GF(2),GF(3),GF(4) WRITE (6,FMT='(''1/Q = '', 3e13.5)') 1.0d0/Q1,1.0d0/Q2, * 1.0d0/Q3 END IF END IF C C .....Calculate direction cosines of a vector perpendicular to C the line-of-sight with respect to the apex coordinates..... TT = -DSIN(DTR*(180.D0-AZ)) TP = DCOS(DTR*(180.D0-AZ)) TR = 0.D0 CALL VCTCNV(TV(1),TV(2),TV(3),SV(1),SV(2),SV(3),TR,TT,TP,SR, * 90.D0-SLATGC,SLON,2) COST1 = SPROD(DQDX(1,1),TV)/ECN1M COST2 = SPROD(DQDX(1,2),TV)/ECN2M COST3 = SPROD(DQDX(1,3),TV)/ECN3M 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 RCOR(16) = GF(1) RCOR(17) = GF(2) RCOR(18) = GF(3) RCOR(19) = GF(4) RCOR(20) = CTAB(1) RCOR(21) = CTAB(2) RCOR(22) = CTAB(3) RCOR(23) = CTAB(4) RCOR(24) = CTAB(5) RCOR(25) = CTAB(6) RCOR(26) = COST1 RCOR(27) = COST2 RCOR(28) = COST3 RCOR(29) = AINC RCOR(30) = DEC RCOR(31) = GCLAT RCOR(32) = ASPCT RCOR(33) = CPLAT RCOR(34) = CGDLAT RCOR(35) = CPLON RCOR(36) = GRADAPEX(1) RCOR(37) = GRADAPEX(2) RCOR(38) = GRADAPEX(3) RCOR(39) = GRADXYZ(1) RCOR(40) = GRADXYZ(2) RCOR(41) = GRADXYZ(3) RCOR(42) = GRADMAG(1) RCOR(43) = GRADMAG(2) RCOR(44) = GRADMAG(3) RCOR(45) = GRADMAG(1)*B RCOR(46) = GRADMAG(2)*B RCOR(47) = GRADMAG(3)*B RCOR(45) = GRADMAG(1) RCOR(46) = GRADMAG(2) RCOR(47) = GRADMAG(3) RCOR(48) = CTAB(7) RCOR(49) = CTAB(8) RCOR(50) = CTAB(9) RCOR(51) = EFCOSCN1 RCOR(52) = EFCOSCN2 RCOR(53) = EFCOSCN3 RCOR(54) = ARC C RETURN END