SUBROUTINE NFTHERY ( DATMC, DIONC, DLPGR, EARTH, EPBASE, SITEP, & & SITEV, SITEA, SUN, STAR, STARdt, XMOON, AT, T0_T1, & & R1, R1dt, R1mag, R1magdt, R2, R2dt, R2mag, R2magdt, & & STAR12, STAR12dt ) Implicit none ! ! Routine THERY computes the theoretical values of delay and delay rate ! using the information generated by the modules. Subroutine CONSEN is ! called to perform the computations using the Eubanks "Consensus" ! relativity model. The Shapiro and Hellings relativity model computations ! have been removed with Calc 9.0 and are no longer used here. ! ! References: None ! ! Program Interface ! Input variables: ! 1. DATMC(2,2) - The contributions to the delay and delay rate due to ! tropospheric refraction at each site. (sec, sec/sec) !** 2. DAXOC(2,2) - The contributions to the delay and rate due to the !** antenna axis offsets. First index runs over sites, !** the second runs over delay and rate (sec,sec/sec). ! 3. DIONC(2) - The contributions to the delay and rate due to ! ionospheric refraction at each site. (sec, sec/sec) ! 4. DLPGR - The CT time derivative of the long period terms in ! 'AT MINUS CT' offset. (sec/sec) (Not used) ! 5. EARTH(3,3) - The solar system barycentric Earth position ! velocity, and acceleration vectors. The first index ! runs over the vector components and the second runs ! over the time derivatives. (m, m/sec, m/sec**2) ! 6. EPBASE(3,2) - The J2000.0 baseline position and velocity vectors. ! (m, m/sec) ! 7. SITEA(3,2) - The J2000.0 geocentric acceleration vectors of each ! observing site. (m/sec**2) ! 8. SITEP(3,2) - The J2000.0 geocentric position vectors of each ! observing site. (m) ! 9. SITEV(3,2) - The J2000.0 geocentric velocity vectors of each ! observing site. (m/sec) ! 10. SUN(3,2) - The J2000.0 geocentric Sun position and velocity ! vectors. (m, m/sec) ! 11. XMOON(3,2) - The J2000.0 geocentric Moon position and velocity ! vectors. (m, m/sec) ! 12. STAR(3) - The J2000.0 source unit vector. (unitless) ! 13. AT - The Atomic Time fraction of the Atomic Time Day ! (TAI). (days) (Not used) ! ! Common blocks used - ! INCLUDE 'ccon.i' ! Variables 'from': ! 1. KTHEC - The Theory routine flow control flag. ! (No longer has any function.) ! 2. KTHED - The theory routine debug control flag. ! 3. KRELC - The relativity module flow control flag. ! 0 --> Gravitational bending used. ! NOT 0 --> Gravitational bending not used. ! INCLUDE 'put2s.i' ! Variables from: ! 2. DAXOC(2,2) - The contributions to the delay and rate due to the ! antenna axis offsets. First index runs over sites, ! the second runs over delay and rate (sec,sec/sec). ! Variables to: ! 1. CONDEL(2) - The theoretical delay from the Consensus ! model in two pieces in units of MICROSECONDS. ! The 1st is the integer microseconds and the ! 2nd is the submicroseconds portion. ! 2. CONRAT - The theoretical delay rate from the Consensus ! model/ (sec/sec). ! 3. SUN_CNTRB(2) - The solar gravitational bending delay and rate ! terms from the Consensus model. (sec, sec/sec) ! 4. CON_CNTRB(2) - The total gravitational bending delay and rate ! terms from the Consensus model. (sec, sec/sec) ! 5. CON_PART(2) - The total Consensus delay and rate partials ! with respect to Gamma. (sec, sec/sec) ! 6. Sunplus(2) - Higher order solar bending delay and rate ! contributions, as defined in IERS Conventions ! (1996), page 91, eqn. 14. (sec, sec/sec) It ! has also been added to the total solar bending ! delay and rate contibutions. ! 7. Bend_par(2) - The relativistic bending delay and rate partials ! with respect to Gamma. (sec, sec/sec) ! ! Program specifications - INTEGER*4 I REAL*8 DATMC(2,2), DIONC(2), SUN(3,2), EARTH(3,3), & & EPBASE(3,2), SITEA(3,2), STAR(3), SITEP(3,2), SITEV(3,2), & & XMOON(3,2), AT, SUNHAT(3), DLPGR, STARdt(3), STAR12(3,2), & & STAR12dt(3,2) Real*8 delta_t_grav, d_delta_t_grav, & & tg2_tg1, dtg2_tg1, xtg2_tg1, & & delta_t_grav_Sun, d_delta_t_grav_Sun, CONSENSUS(2), & & con_cont(2) Real*8 R1(3), R2(3), R1dt(3), R2dt(3), R1mag, R2mag, R1magdt, & & R2magdt, T0_T1 Real*8 tr2_tr1, dtr2_tr1 ! ! 4.2.4 DATA BASE ACCESS - ! !** 'PUT' VARIABLES: !** 1. CONDEL(2) - The theoretical delay from the Consensus !** model in two pieces in units of MICROSECONDS. !** The 1st is the integer microseconds and the !** 2nd is the submicroseconds portion. !** 2. CONRAT - The theoretical delay rate from the Consensus !** model/ (sec/sec). !** 3. SUN_CNTRB(2) - The solar gravitational bending delay and rate !** terms from the Consensus model. (sec, sec/sec) !** 4. CON_CNTRB(2) - The total gravitational bending delay and rate !** terms from the Consensus model. (sec, sec/sec) !** 5. CON_PART(2) - The total Consensus delay and rate partials !** with respect to Gamma. (sec, sec/sec) !** 6. Sunplus(2) - Higher order solar bending delay and rate !** contributions, as defined in IERS Conventions !** (1996), page 91, eqn. 14. (sec, sec/sec) It !** has also been added to the total solar bending !** delay and rate contibutions. !** 7. Bend_par(2) - The relativistic bending delay and rate partials !** with respect to Gamma. (sec, sec/sec) ! ! ACCESS CODES: ! 1. 'CONSNDEL' - The database access code for the theoretical ! delay using the Eubanks' Consensus model. ! 2. 'CONSNRAT' - The database access code for the theoretical ! delay rate using the Eubanks' Consensus model. ! 3. 'SUN CONT' - The database access code for the Consensus model ! Solar gravitational bending delay and rate terms. ! 4. 'CON CONT' - The database access code for the Consensus model ! TOTAL gravitational bending delay and rate terms. ! 5. 'CONSPART' - Database access code for the partial derivatives ! of the Consensus model delays and rates with ! respect to Gamma. ! 6. 'SUN2CONT' - Database access code for the additional solar ! bending due to higher order relativistic effects. ! This term is already in the theoretical, so to ! remove it, SUBTRACT the values in this access code. ! 7. 'BENDPART' - Database access code for the partial derivatives ! of the Consensus model gravitational bending ! delay and rate contributions with respect to Gamma. ! ! Subroutine interface - ! Caller subroutine: DRIVR ! Called subroutines: PUT4, CONSEN ! ! Program variables - ! 1. tg2_tg1 - The geometric delay corrected for relativistic ! effects using the Consensus model. ! 2. dtg2_tg1 - The time derivative of the geometric delay ! corrected for relativistic effects using the ! Consensus model. ! 3. delta_t_grav - The total differential gravitational time delay, ! or "bending delay" from the Consensus model. ! 4. d_delta_t_grav-The time derivative of the total differential ! gravitational time delay, or "bending delay," ! from the Consensus model. !** 5. SUN_CNTRB(2) - The solar gravitational bending delay and rate !** contributions from the Consensus model. (s, s/s). !** 6. CON_CNTRB(2) - The total gravitational bending delay and rate !** contributions from the Consensus model. (s, s/s). ! 7. CON_PART(2) - The partial derivatives of the Consensus model ! delay and rate with respect to Gamma (s and s/s). !** 8. Sunplus(2) - Higher order solar bending delay and rate !** contributions, as defined in IERS Conventions !** (1996), page 91, eqn. 14 (sec, sec/sec). It !** has also been added to the total solar bending !** delay and rate contibutions. !** 9. Bend_par(2) - The relativistic bending delay and rate partials !** with respect to Gamma. (sec, sec/sec) ! ! 4.2.9 PROGRAMMER - DALE MARKHAM 01/17/77 ! PETER DENATALE 07/20/77 ! CHOPO MA / DAVID GORDON 04/12/84 ! DAVID GORDON 06/19/84 REMOVED ATMOSPHERE. ! DAVID GORDON 07/31/84 RELATIVISTIC CORRECTIONS. ! DAVID GORDON 01/03/85 ADDED ATMOSPHERE AFTER MODS TO ! ATMOSPHERE FLAGS. ! SAVITA GOEL 06/03/87 CDS FOR A900. ! GREGG COOKE 05/01/89 NEW MODEL FROM I.SHAPIRO. ! JIM RYAN 06/06/89 Some bugs fixed and documentation ! modified. ! Jim Ryan 09/13/89 Robertson code removed and Hellings ! added. ! Jim Ryan 89.10.05 CPHYS common made an include file ! Jim Ryan 89.10.09 Relativity partials and contributions ! moved here. ! Jim Ryan 89.11.20 Shapiro algorithm modified to make ! it reflect Ryan's memo ! Jim Ryan 89.12.12 UNIX-like database interface ! implemented. ! Jim Ryan 89.12.14 Addional delay information stored ! and Helling rate implemented. ! C Ma 90.08.10 Corrections to documentation. ! Jim Ryan 90.11.20 Debug statement fixed and comments ! cleaned up. ! T. Marshall Eubanks & Brent Archinal 91.05.10 HELL EMS ! fixed, debug statements and comments modified. ! A900 and HP-UX versions consolidated except for ! dbh calls and ILUOUT. ! Jim Ryan 91.05.28 A few, mostly cosmetic changes made. ! David Gordon 93.04.27 GMEARTH put into cphys.i and ! defined in cinit.f; definition removed here. ! David Gordon 93.08.02 Thery modified to call subroutine ! Consen (previously called by DRIVR), and to ! do the puts for the Consensus delay and rate. ! David Gordon 93.12.22 Fixed up debug output. ! David Gordon 93.12.27 Added access code CONSCONT, correction ! to convert Hellings to Comsensus theoreticals. ! David Gordon 94.01.21 Added access code SHAPCONT, correction ! to convert Hellings to Shapiro theoreticals. ! David Gordon 94 Feb/March - Modified axis offset correction ! for use here. ! David Gordon 94.10.05 Corrected error in computing second ! half of Lcode for Shapiro delay contribution. ! Many unused variables removed. ! David Gordon 95.10.11 Minor correction to Hellings model ! when gravitational bending turned OFF (KRELC=1). ! David Gordon 96.01.30 Added Consensus model solar bending ! term (Sun_cntrb(2) and the L-code 'SUN CONT'). ! Removed 'CON PART' L-code (Solar bending ! partials w.r.t. Gamma) and replaced it with ! 'CONSPART' (Total Consensus delay and rate ! partials w.r.t. Gamma). ! David Gordon 98.08.18 Shapiro and Hellings models removed. ! David Gordon 98.11.16 Added Sunplus to CALL CONSEN argument ! list. Contains the delay and rate contributions ! due to higher order solar bending, as defined ! in the 1996 IERS Conventions, p. 91, eqn. 14. ! Added PUTR of 'SUN2CONT' to put it in the ! data bases. It is already included in the ! theoretical, therefore to remove its effects, ! you must SUBTRACT it from the theoretical. ! David Gordon 98.12.16 Added Bend_par(2), the partials of the ! gravitional bending contributions with respect ! to Gamma. ! Jim Ryan Sept 2002 Integer*2/4 mods. ! David Gordon Jan. 2013 Moved PUT's into subroutine PUT_TH. ! ! THERY program structure !-------------------------------------------------------------------------- ! ! write(6,*) 'NFTHERY/R1mag,R1magdt ', R1mag,R1magdt ! write(6,*) 'NFTHERY/R2mag,R2magdt ', R2mag,R2magdt ! ! Call subroutine Sekido to compute the delay and delay rate based on the ! Eubanks' Consensus relativity model. We put it in its own subroutine ! because it is a very long and detailed set of computations. 93JUL29, DG ! Call SEKIDO ( DATMC, EARTH, EPBASE, SITEP, SITEV, SITEA, SUN, & & XMOON, STAR, STARdt, T0_T1, R1, R1dt, R1mag, R1magdt, & & R2, R2dt, R2mag, R2magdt, STAR12, STAR12dt, & & tg2_tg1, dtg2_tg1 ) ! ! print *,' Delay,Rate: ', tg2_tg1*1.D12, dtg2_tg1*1.D12 ! Add the axis offset corrections ! TG2_TG1 = tg2_tg1 + DAXOC(1,1) + DAXOC(2,1) ! dtg2_tg1 = dtg2_tg1 + DAXOC(1,2) + DAXOC(2,2) ! print *,' Delay,Rate: ', tg2_tg1*1.D12, dtg2_tg1*1.D12 ! print *,' NFTHERY ' ! print *,'Delay(ps),Rate(ps/sec): ', tg2_tg1*1.D12,dtg2_tg1*1.D12 ! print *,'Delay,Rate(AXOF added): ', tg2_tg1,dtg2_tg1 ! tg2_tg1 = tg2_tg1 + DATMC(1,1) + DATMC(2,1) ! dtg2_tg1 = dtg2_tg1 + DATMC(1,2) + DATMC(2,2) ! print *,'Delay,Rate(AXOF,ATMC added): ', tg2_tg1,dtg2_tg1 ! print *,'DATMC (sec): ', DATMC ! print *, ' ' ! CONSENSUS(1) = tg2_tg1 ! Sekido delay CONSENSUS(2) = dtg2_tg1 ! Sekido rate ! ! print *,' THERY/Delay (psec): ', CONSENSUS(1)*1.D12 ! print *,' THERY/Rate (psec/sec): ', CONSENSUS(2)*1.D12 ! ! Convert the theoretical delay to microseconds and split the double ! precision results into an integer microseconds portion and a ! submicroseconds portion. ! xtg2_tg1 = tg2_tg1 * 1.D6 ! CONDEL(1) = IDINT(xtg2_tg1 ) ! CONDEL(2) = xtg2_tg1 - CONDEL(1) ! CONRAT = dtg2_tg1 ! !!!!!!!!!!!! Debug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! write(6,8)' CONDEL, CONRAT ', CONDEL, CONRAT ! write(6,8)' CONSENSUS ',CONSENSUS ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Call RANGE ( DATMC, SITEP, SITEV, SITEA, SUN, & ! & T0_T1, R1, R1dt, R1mag, R1magdt, & ! & R2, R2dt, R2mag, R2magdt, & ! & tr2_tr1, dtr2_tr1 ) ! !!!!!!!!!!!! Debug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !-------------------------------------------------------------------------- ! Check KTHED to determine if debug output is neccessary. IF ( KTHED .ne. 0 ) Then WRITE (6, 9100 ) 9100 FORMAT (1X, "Debug output for subroutine THERY." ) 8 FORMAT(A,4D25.16/(7X,5D25.16)) write(6,8)' tg2_tg1, dtg2_tg1 ',tg2_tg1, dtg2_tg1 write(6,8)' CONDEL, CONRAT ', CONDEL, CONRAT write(6,8)' CON_CNTRB ',CON_CNTRB write(6,8)' SUN_CNTRB ',SUN_CNTRB write(6,8)' CON_PART ',CON_PART write(6,8)' CONSENSUS ',CONSENSUS write(6,8)' Sunplus ',Sunplus Endif ! Return END !********************************************************************** SUBROUTINE SEKIDO ( DATMC, EARTH, EPBASE, SITEP, SITEV, & & SITEA, SUN, XMOON, STAR, STARdt, T0_T1, R1, R1dt, & & R1mag, R1magdt, R2, R2dt, R2mag, R2magdt, STAR12, & & STAR12dt, tg2_tg1, dtg2_tg1 ) IMPLICIT NONE ! ! Routine SEKIDO computes the theoretical values of delay and delay ! rate using the Sekido and Fukushima near-field model. ! ! Reference: M. Sekido and T. Fukushima, 'A VLBI delay model for ! radio sources at a finite distance', J. Geodesy, ! vol. 80, pp. 137-149, 2006. ! ! Input Variables: ! 1. DATMC(2,2) - The contributions to the delay and delay rate due ! to tropospheric refraction at each site. (s, s/s) ! 2. EARTH(3,3) - The solar system barycentric Earth position, ! velocity, and acceleration vectors. The first ! index runs over the vector components and the ! second runs over the time derivatives. ! (m, m/sec, m/sec**2) ! 3. EPBASE(3,2) - The J2000.0 baseline position and velocity ! vectors. (m, m/sec) ! 4. SITEP(3,2) - The J2000.0 geocentric position vectors of each ! observing site. (m) ! 5. SITEV(3,2) - The J2000.0 geocentric velocity vectors of each ! observing site. (m/sec) ! 6. SITEA(3,2) - The J2000.0 geocentric acceleration vectors of ! each observing site. (m/sec**2) ! 7. SUN(3,2) - The J2000.0 geocentric Sun position and velocity ! vectors. (m, m/sec) ! 8. XMOON(3,2) - The J2000.0 geocentric Moon position and velocity ! vectors. (m, m/sec) ! 9. STAR(3) - The J2000.0 source unit vector. (unitless) ! ! Output Variables: ! 1. tg2_tg1 - The geometric delay corrected for relativistic ! effects (but not atmospheric refraction) using ! the Consensus model. ! 2. dtg2_tg1 - The time derivative of the geometric delay ! corrected for relativistic effects using the ! Consensus model. ! 3. delta_t_grav - The total differential gravitational time delay, ! or "bending delay" from the Consensus model. ! 4. d_delta_t_grav- The time derivative of the total differential ! gravitational time delay, or "bending delay" from ! the Consensus model. ! 5. delta_t_grav_Sun-The Sun's contribution to the "bending delay" ! from the Consensus model. Now includes ! Sunplus(1). ! 6. d_delta_t_grav_Sun-The time derivative of the Sun's contribution ! to the "bending delay" from the Consensus model. ! Now includes Sunplus(2). ! 7. Con_part - Partial derivatives of the Consensus delays and ! rates with respect to Gamma. ! 8. Sunplus(2) - Higher order solar bending delay and rate ! contributions, as defined in IERS Conventions ! (1996), page 91, eqn. 14. (sec, sec/sec) It ! has also been added to the total solar bending ! delay and rate contibutions. ! 9. Bend_par(2) - The relativistic bending delay and rate partials ! with respect to Gamma. (sec, sec/sec) ! ! Common blocks used - ! INCLUDE 'cphys11.i' ! Variables 'from': ! 1. VLIGHT - The velocity of light in vacuum. (m/sec) ! 2. VLIGHT2 - The velocity of light squared. (m/sec)**2 ! 3. VLIGHT3 - The velocity of light cubed. (m/sec)**3 ! 4. GAMMA - The post Newtonian expansion parameter which ! affects light bending. (1.0 used here). (unitless) ! 5. GMSUN, GMMOON, GMEARTH, GMPLANET(7) - Gravitational constant ! times the masses of the Sun, Moon, Earth, and the ! other planets except Pluto. ! 6. REARTH - The equatorial radius of the Earth. (meters) ! INCLUDE 'csolsys11.i' ! Variables 'from': ! 1. SPLANET(3,2,7) - The J2000.0 Solar System Barycentric positions ! and velocities of all planets except the Earth ! and Pluto. (meters, meters/sec) The first index ! runs over X,Y, and Z, the second runs over ! position and velocity, and the third runs over ! the planets, where ! 1 = Mercury ! 2 = Venus ! 3 = Mars ! 4 = Jupiter ! 5 = Saturn ! 6 = Uranus ! 7 = Neptune ! ! 2. GPLANET(3,2,7) - The J2000.0 Geocentric positions and velocities ! of all planets except the Earth and Pluto. ! (meters, meters/sec) The first index runs over ! X,Y, and Z, the second runs over position and ! velocity, and the third runs over the planets, ! where ! 1 = Mercury ! 2 = Venus ! 3 = Mars ! 4 = Jupiter ! 5 = Saturn ! 6 = Uranus ! 7 = Neptune ! Real*8 PI, TWOPI, HALFPI, CONVD, CONVDS, CONVHS, SECDAY COMMON / CMATH / PI, TWOPI, HALFPI, CONVD, CONVDS, CONVHS, SECDAY ! INCLUDE 'ccon.i' ! Variables 'from': ! 1. KTHEC - The Theory routine flow control flag. ! (No longer has any function.) ! 2. KTHED - The theory routine debug control flag. ! 3. KRELC - The relativity module flow control flag. ! 0 --> Gravitational bending used. ! NOT 0 --> Gravitational bending not used. ! INCLUDE 'cobsn11.i' ! Variables from: ! 1. Nzero - Set to 1 or 2 if station 1 or 2 is at the geocenter, ! otherwise equals zero. For correlator usage. ! INCLUDE 'put2s.i' ! Variables to: ! 5. CON_PART(2) - The total Consensus delay and rate partials ! with respect to Gamma. (sec, sec/sec) ! 6. Sunplus(2) - Higher order solar bending delay and rate ! contributions, as defined in IERS Conventions ! (1996), page 91, eqn. 14. (sec, sec/sec) It ! has also been added to the total solar bending ! delay and rate contibutions. ! 7. Bend_par(2) - The relativistic bending delay and rate partials ! with respect to Gamma. (sec, sec/sec) ! INCLUDE 'd_input.i' ! Variables from: ! 1. SPpxyz(3) - ! 2. SPpvyz(3) - ! ! ! Program Specifications - REAL*8 DATMC(2,2), SUN(3,2), EARTH(3,3), EPBASE(3,2), SITEA(3,2), & & STAR(3), STARdt(3), SITEP(3,2), SITEV(3,2), XMOON(3,2), & & STAR12(3,2), STAR12dt(3,2), DOTP, VECMG Real*8 R1(3), R2(3), R1dt(3), R2dt(3), R1mag, R2mag, R1magdt, & & R2magdt, T0_T1 Real*8 XsubEarth(3), VsubEarth(3), AsubEarth(3) Real*8 x_sub1t1(3), x_sub2t1(3), w_sub1(3), w_sub2(3) Real*8 a_sub1(3), a_sub2(3) Real*8 unit_K(3), b_nought(3),db_nought(3), unitK1(3) Real*8 x_subSun(3), x_subMoon(3), v_subSun(3), v_subMoon(3) Real*8 XsubSun(3), XsubMoon(3), VsubSun(3), VsubMoon(3) Real*8 R_Earth_Sun(3), R_Earth_Moon(3) Real*8 V_Earth_Sun(3), V_Earth_Moon(3) Real*8 Xsub1(3), Xsub2(3), t_sub1, t_sub2 Real*8 dXsub1(3), dXsub2(3) Real*8 x1Sun(3), x1Moon(3), x1Planet(3,7) Real*8 del_t_Sun, del_t_Moon, del_t_Planet(7) Real*8 XSunt1J(3), XMoont1J(3), XPlant1J(3,7) Real*8 dXSunt1J(3), dXMoont1J(3), dXPlant1J(3,7) Real*8 R1Sunt1(3), R1Moont1(3), R1Plant1(3,7) Real*8 R2Sunt1(3), R2Moont1(3), R2Plant1(3,7) Real*8 dR1Sunt1(3), dR1Moont1(3), dR1Plant1(3,7) Real*8 dR2Sunt1(3), dR2Moont1(3), dR2Plant1(3,7) Real*8 delta_t_grav, delta_t_grav_Sun, delta_t_grav_Moon, & & delta_t_grav_Earth, delta_t_grav_Plan(7), & & delta_t_grav_Planets Real*8 d_delta_t_grav, d_delta_t_grav_Sun, d_delta_t_grav_Moon, & & d_delta_t_grav_Earth, d_delta_t_grav_Plan(7), & & d_delta_t_grav_Planets Real*8 delgrav_Sun,ddelgrav_Sun, delgrav_Moon, ddelgrav_Moon, & & delgrav_Plan(7), ddelgrav_Plan(7), delgrav, ddelgrav, & delgrav_Earth, ddelgrav_Earth Real*8 U, U_Sun, dU, absVEarth Real*8 C_Earth, C_Sun, C_Moon, C_Plan(7) Real*8 term_a, term_b, term_c, term_d, term_e, term_f, & & term_g, term_h, term_a1, term_b1, term_g1, term_h1 Real*8 dterm_a, dterm_b, dterm_c, dterm_d, dterm_e, dterm_f, & & dterm_g, dterm_h, dterm_a1, dterm_b1, dterm_g1, dterm_h1 Real*8 term1, term2a, term2b, term2c, term2d, term2, & & term2bcd, term3a, term3b, term3, term4, & & term123, vec_sum(3), tv2_tv1 Real*8 dterm1, dterm2a, dterm2b, dterm2c, dterm2d, dterm2, & & dterm2bcd, dterm3a, dterm3b, dterm3, dterm4, & & dterm123, dvec_sum(3), dtv2_tv1 Real*8 V_w1(3), V_w2(3), K_V_w1, K_V_w2 !! !** Real*8 tg2_tg1, dtg2_tg1, Con_part(2), Bend_par(2) Real*8 tg2_tg1, dtg2_tg1 Real*8 w2_w1(3), a2_a1(3), K_dotw2w1, dK_dotw2w1 ! Real*8 k_sub1(3), k_sub2(3) Real*8 KdotB, dKdotB, VdotB, dVdotB, KdotV, dKdotV Real*8 delta_t2_t1, d_delta_t2_t1, term10a, term10b, term10c, & & dterm10a, dterm10b, dterm10c, del_del_t(2,2) Real*8 x_delta_t2_t1, dx_delta_t2_t1 Real*8 vecmg1, vecmg2, vecmg0, dvecmg0 Real*8 N_hat(3), dN_hat(3), Vmag_S, CSun1, NplusK(3), V1, dV1 !! & V2, dV2 Real*8 R0sun(3), dR0sun(3), R1J, R2J, R10, R20, R0J, & & delTpNsun, R0moon(3), dR0moon(3), delTpNmoon, R0plan(3,7), & & dR0plan(3,7), delTpNPlan(7), delTpN, Vec3c(3), & & term3c, term4a, H, H1(3), H2, H3, dH1a(3), dH1b(3),dH1(3), & & dH2, dH3, dH, V2(3), dV2(3), dVec3c(3), dterm3c, dterm4a Real*8 dR1J, dR2J, dR10, dR20, dR0J, GMp, DelTpNPl, DelTpNPldt, & & delTpNsundt, DelTpNmoondt, DelTpNPlandt(7), delTpNdt Real*8 VsubEart0(3), AsubEart0(3) INTEGER*4 I, L, K ! ! Subroutines used: DOTP, IDINT, VECAD, VECMG ! ! Constants used - VLIGHT - Speed of light (m/s) ! VLIGHT2 - Speed of light squared ! VLIGHT3 - Speed of light cubed ! GMMOON - GM of the Moon. ! GMEARTH - GM of the Earth. ! GMSUN - GM of the Sun. ! GMPlanet(7) - GM's of all planets except Earth and Pluto ! ! Program variables - ! ! 1. XsubEarth(3) = Barycentric radius vector of the geocenter (meters). ! 2. VsubEarth(3) = Barycentric velocity vector of the geocenter (m/sec). ! 3. AsubEarth(3) = Barycentric acceleration of the geocenter (m/sec**2). ! ! 4. x_sub1t1(3) = Geocentric radius vector of receivers 1 and 2 at the ! 5. x_sub2t1(3) geocentric time t_sub1 (meters). ! ! 6. w_sub1(3) = Geocentric velocity of receiver 1 (m/sec). ! 7. w_sub2(3) = Geocentric velocity of receiver 2 (m/sec). ! ! 8. a_sub1(3) = Geocentric acceleration of receiver 1 (m/sec). ! 9. a_sub2(3) = Geocentric acceleration of receiver 2 (m/sec). ! ! 10. x_subSun(3) = Geocentric radius vector of the Sun (meters). ! 11. x_subMoon(3) = Geocentric radius vector of the Moon (meters). ! 12. GPLANET(l,1,k)= Geocentric radius vectors of the planets (meters). ! (l = X,Y,Z; k = planet #) ! ! 13. v_subSun(3) = Geocentric velocity vector of the Sun (m/sec). ! 14. v_subMoon(3) = Geocentric velocity vector of the Moon (m/sec). ! 15. GPLANET(l,2,k)= Geocentric velocity vector of the planets (m/sec). ! (l = X,Y,Z; k = planet #) ! ! 16. Xsub1(3) = Barycentric radius vector of receiver 1 (meters). ! 17. Xsub2(3) = Barycentric radius vector of receiver 2 (meters). ! 18. dXsub1(3) = Barycentric velocity vector of receiver 1 (m/sec). ! 19. dXsub2(3) = Barycentric velocity vector of receiver 2 (m/sec). ! ! 20. x1Sun(3) = Vector from receiver 1 to the Sun at time t_sub1 ! 21. x1Moon(3) = Vector from receiver 1 to the Moon at time t_sub1 ! 22. x1Planet(3,7) = Vector from receiver 1 to the planets at time t_sub1 ! ! 23. XsubSun(3) = Barycentric radius vector of the Sun/Moon/planets ! 24. XsubMoon(3) ! 25. SPLANET(l,1,k) ! ! 26. VsubSun(3) = Barycentric velocity vector of the Sun/Moon/planets ! 27. VsubMoon(3) ! 28. SPLANET(l,2,k) ! ! 29. XSunt1J(3) = SSBC radius vector to the Sun/Moon/planets ! 30. XMoont1J(3) at the time of closest approach. ! 31. XPlant1J(3,7) (meters) ! ! 32. R_Earth_Sun(3) = Vector from the Sun to the geocenter ! 33. R_Earth_Moon(3)= Vector from the Moon to the geocenter ! 34. -GPLANET(l,1,k)= Vector from planets to the geocenter ! ! 35. V_Earth_Sun(3) = Velocity Vector of the geocenter from the Sun ! ! 36. R1Sunt1(3) = Vector from the Sun/Moon/planets to receiver ! 37. R1Moont1(3) = 1 at the time of closest approach. ! 38. R1Plant1(3,7) = (meters) ! ! 39. R2Sunt1(3) = Vector from the Sun/Moon/planets to receiver ! 40. R2Moont1(3) = 2 at the time of closest approach. ! 41. R2Plant1(3,7) = (meters) ! ! 42. t_sub1 = Time of arrival of the signal at receivers 1 and 2. ! 43. t_sub2 [Note: not needed here so we don't actually define ! them.] ! ! 44. del_t_Sun = Difference between arrival time at receiver 1, t_sub1, ! 45. del_t_Moon and time of closest approach to Sun/Moon/planets ! 46. del_t_Planet(7) ( > 0 if closest approach is before arrival) ! ! 47. unit_K(3) = Barycentric unit vector to the source (in the absense ! of gravitational or aberrational bending). ! ! 48. k_sub1(3) = Aberrated unit vector from station 1 to the source. ! 49. k_sub2(3) = Aberrated unit vector from station 2 to the source. ! [Note: not computed here, see atmosphere module.] ! ! 50. b_nought(3) = A priori geocentric baseline vector at time t_sub1 ! (Baseline vector and it's velocity. Defined from ! site #1 to site #2 (M,M/S).) ! ! 51. U = Sun's potential ! ! 52. delgrav_Sun ! 53. delgrav_Moon = The differential gravitational time delay for ! 54. delgrav_Plan(7) the Sun/Moon/planets/Earth. ! 55. delgrav_Earth ! ! 56. delta_t_grav_Planets = Sum of bending delays for the planets ! ! 57. delta_t_grav = The total differential gravitational time delay, ! or "bending delay." ! ! 58. tv2_tv1 = The theoretical vacuum delay computed using the ! Consensus model. (sec) ! ! 59. tg2_tg1 = The total theoretical geometric delay (doesn't ! include tropospheric refraction effects which are ! usually added in Solve) computed using the Consensus ! model. (sec) ! ! 4.2.9 PROGRAMMER - DAVID GORDON 04/22/93 thru 08/02/93 - Written and debuged ! D. Gordon Nov. 1993 - Modified to use all planets ! except Pluto. ! D. Gordon Jan-Mar 1994 - Modified for axis offset ! correction. Not using Eubank's Step 10. ! D. Gordon 10.05.94 Many unused variable and much unused ! code removed. ! D. Gordon 96.01.30 Added section to compute partial ! derivatives of delay and rate w.r.t. Gamma. ! D. Gordon 96.02.09 Corrected typo in Step 2 debug printout. ! D. Gordon 98.08.18 Mods for geocenter station. ! D. Gordon 98.11.16 Added computation of higher order ! solar bending term, from IERS Conventions ! (1996), page 91, eqn. 14. This term is now ! added to the total solar bending term. It ! only becomes significant within about 2 ! degrees of the Sun. ! D. Gordon 98.12.16 Added computation of Bend_par(2), the ! partial derivatives of the gravitational ! bending contributions with respect to Gamma. ! ! Write(6,*) ' ' ! Write(6,*) ' SEKIDO: ' ! write(6,*) 'SEKIDO/R1mag,R1magdt ', R1mag,R1magdt ! write(6,*) 'SEKIDO/R2mag,R2magdt ', R2mag,R2magdt ! ! CONSEN program structure: ! ! Compute the theoretical delay and delay rate using Eubanks's Consensus ! Relativity model. !_____________________________________________________________________________ ! ! Copy CALC variables into variables with names which mimic the variables ! in the consensus paper. VLIGHT is used for the velocity of light because ! 'C' is a poor variable name that could be easily lost. ! Do l=1,3 ! XsubEarth(l) = EARTH(l,1) ! SSBC Earth position VsubEarth(l) = EARTH(l,2) ! ditto velocity AsubEarth(l) = EARTH(l,3) ! ditto acceleration ! VsubEart0(l) = 0.0D0 ! ditto velocity AsubEart0(l) = 0.0D0 ! ditto acceleration ! VsubEart0(l) = EARTH(l,2) ! ditto velocity ! AsubEart0(l) = EARTH(l,3) ! ditto acceleration ! x_sub1t1(l) = SITEP(l,1) ! Site 1 geocentric position w_sub1(l) = SITEV(l,1) ! ditto velocity a_sub1(l) = SITEA(l,1) ! ditto acceleration ! x_sub2t1(l) = SITEP(l,2) ! Site 2 geocentric position w_sub2(l) = SITEV(l,2) ! ditto velocity a_sub2(l) = SITEA(l,2) ! ditto acceleration ! unit_K(l) = STAR(l) ! J2000.0 pseudo source vector unitK1(l) = STAR12(l,1) ! x_subSun(l) = SUN(l,1) ! Geocentric Sun position v_subSun(l) = SUN(l,2) ! ditto velocity ! x_subMoon(l) = XMOON(l,1) ! Geocentric Moon position v_subMoon(l) = XMOON(l,2) ! ditto velocity ! b_nought(l) = -EPBASE(l,1) ! Baseline vector from site 1 to site 2 db_nought(l) = -EPBASE(l,2) ! Time derivative of baseline vector ! R_Earth_Sun(l) = -SUN(l,1) V_Earth_Sun(l) = -SUN(l,2) ! = -v_subSun(l) also R_Earth_Moon(l) = -XMOON(l,1) ! ! SSBC radius vectors of Sun and Moon: XsubSun(l) = XsubEarth(l) + x_subSun(l) XsubMoon(l) = XsubEarth(l) + x_subMoon(l) ! ! SSBC velocity vectors of Sun and Moon: VsubSun(l) = VsubEarth(l) + v_subSun(l) VsubMoon(l) = VsubEarth(l) + v_subMoon(l) ! Enddo ! ! Atomic time, TAI, at receiver #1: ! t_sub1 = ! Not needed here ! If (KTHED .ne. 0) Then write(6,'(/,15x," Debug output for subroutine CONSEN",/)') write(6,8)' GMSUN ', GMSUN write(6,8)' GMMOON ', GMMOON write(6,8)' GMEARTH', GMEARTH write(6,'("GM-Planets = ",3d25.16,/,5x,4d25.16)') GMPLANET write(6,'(/,"XsubEarth:",3D23.14)') XsubEarth write(6,'("VsubEarth:",3D23.14)') VsubEarth write(6,'("AsubEarth:",3D23.14)') AsubEarth write(6,'("x_sub1t1:",3D23.14)') x_sub1t1 write(6,'("w_sub1:",3D23.14)') w_sub1 write(6,'("a_sub1:",3D23.14)') a_sub1 write(6,'("x_sub2t1:",3D23.14)') x_sub2t1 write(6,'("w_sub2:",3D23.14)') w_sub2 write(6,'("a_sub2:",3D23.14)') a_sub2 write(6,'("unit_K:",3D23.14)') unit_K write(6,'("unitK1:",3D23.14)') unitK1 write(6,'("x_subSun:",3D23.14)') x_subSun write(6,'("x_subMoon:",3D23.14)') x_subMoon write(6,'("v_subSun:",3D23.14)') v_subSun write(6,'("v_subMoon:",3D23.14)') v_subMoon write(6,'("b_nought:",3D23.14)') b_nought write(6,'("db_nought:",3D23.14)') db_nought write(6,'("R_Earth_Sun:",3D23.14)') R_Earth_Sun write(6,'("V_Earth_Sun:",3D23.14)') V_Earth_Sun write(6,'("R_Earth_Moon:",3D23.14)') R_Earth_Moon write(6,'("XsubSun:",3D23.14)') XsubSun write(6,'("XsubMoon:",3D23.14)') XsubMoon write(6,'("VsubSun:",3D23.14)') VsubSun write(6,'("VsubMoon:",3D23.14)') VsubMoon Endif ! !****************************************************************************** ! Step 1: Estimate barycentric radius and velocity vectors for stations 1 ! and 2 at time t_sub1 (Equation 6 and time derivative of). ! Do l=1,3 Xsub1(l) = XsubEarth(l) + x_sub1t1(l) Xsub2(l) = XsubEarth(l) + x_sub2t1(l) dXsub1(l) = VsubEarth(l) + w_sub1(l) !derivative dXsub2(l) = VsubEarth(l) + w_sub2(l) !derivative Enddo ! If (KTHED .ne. 0) Then write(6,'(/,15x,"Step 1 dump:")') write(6,8)' Xsub1 ',Xsub1 write(6,8)' Xsub2 ',Xsub2 write(6,8)' dXsub1 ',dXsub1 write(6,8)' dXsub2 ',dXsub2 Endif ! !****************************************************************************** ! Step 2: Estimate the vectors from the Sun, the Moon, and each planet (except ! Earth and Pluto) to receiver 1. ! ! Eq. 5a - Find time of closet approach of the near-field object to the ! gravitating body. [Actually ........................ ! we find only how much earlier (or later) the ........ rays ! passed closest to the gravitating body. Then, if earlier, we ! extrapolate the body's position back to that earlier time. If ! not earlier we just keep the current position.] ! ! J2000.0 vector from receiver #1 to the Sun/Moon/planets: Do l=1,3 x1Sun(l) = XsubSun(l) - Xsub1(l) x1Moon(l) = XsubMoon(l) - Xsub1(l) Enddo Do k=1,7 do l=1,3 x1Planet(l,k) = SPLANET(l,1,k) - Xsub1(l) enddo Enddo ! ! more good stuff KdotB = DOTP(STAR,b_nought) dKdotB = DOTP(STAR,db_nought) + DOTP(STARdt,b_nought) ! Write(6,*) ' ' ! Write(6,*) ' SEKIDO: ' ! Write(6,*) ' KdotB,dKdotB ', KdotB,dKdotB ! ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! Sun: del_t_Sun = DOTP(unitK1,x1Sun) / Vlight ! ! write(6,8)' T0_T1, del_t_Sun(sec) ', T0_T1, del_t_Sun If(del_t_Sun .lt. 0.D0) del_t_Sun = 0.0D0 If(del_t_Sun .gt. T0_T1) del_t_Sun = T0_T1 ! write(6,8)' T0_T1, del_t_Sun(sec) ', T0_T1, del_t_Sun ! Do l=1,3 ! SSBC vector to Sun at time of closest approach and its time derivative: XSunt1J(l) = Xsubsun(l) - VsubSun(l)*del_t_Sun dXSunt1J(l) = Vsubsun(l) !Derivative (approx. - no acceleration) ! ! Vector from Sun to receiver 1 R1Sunt1(l) = Xsub1(l) - XSunt1J(l) dR1Sunt1(l) = dXsub1(l) - dXSunt1J(l) !Derivative ! ! Vector from Sun to receiver 2 R2Sunt1(l) = Xsub2(l) - VsubEarth(l)*KdotB/VLIGHT - XSunt1J(l) dR2Sunt1(l) = dXsub2(l) - AsubEarth(l)*KdotB/VLIGHT - & & VsubEarth(l)*dKdotB/VLIGHT - dXSunt1J(l) !Derivative ! ! Vector from Sun to near-field object R0sun(l) = R1Sunt1(l) + R1(l) dR0sun(l) = dR1Sunt1(l) + R1dt(l) ! Enddo ! write(6,8)' Xsubsun,VsubSun ', Xsubsun,VsubSun ! write(6,8)' Xsub1, XSunt1J ', Xsub1, XSunt1J ! write(6,8)' Xsub2, VsubEarth ', Xsub2, VsubEarth ! write(6,8)' KdotB,VLIGHT ', KdotB,VLIGHT ! write(6,8)' XSunt1J,dXSunt1J ', XSunt1J,dXSunt1J ! write(6,8)' R1Sunt1,dR1Sunt1 ', R1Sunt1,dR1Sunt1 ! write(6,8)' R2Sunt1,dR2Sunt1 ', R2Sunt1,dR2Sunt1 ! write(6,8)' R0sun,R1 ', R0sun,R1 ! R1J = VECMG (R1Sunt1) dR1J = DOTP(R1Sunt1,dR1Sunt1)/R1J R2J = VECMG (R2Sunt1) dR2J = DOTP(R2Sunt1,dR2Sunt1)/R2J R0J = VECMG (R0sun) dR0J = DOTP(R0Sun,dR0Sun)/R0J R10 = R1mag dR10 = R1magdt R20 = R2mag dR20 = R2magdt ! Write(6,*) ' Sun: ' ! Write(6,*) ' R1,R1Sunt1 ', R1, R1Sunt1 ! Write(6,*) ' R0sun ', R0sun ! Write(6,*) ' R2,R2Sunt1 ', R2, R2Sunt1 ! Write(6,*) ' R1J, dR1J ', R1J, dR1J ! Write(6,*) ' R2J, dR2J ', R2J, dR2J ! Write(6,*) ' R0J, dR0J ', R0J, dR0J ! Write(6,*) ' R10, dR10 ', R10, dR10 ! Write(6,*) ' R20, dR20 ', R20, dR20 ! Write(6,*) ' ' ! CALL TGRAV(GMSUN, R0J,R1J,R2J,R10,R20, dR0J,dR1J,dR2J, & & dR10,dR20, delgrav_Sun, ddelgrav_Sun) ! & dR10,dR20, DelTpNsun,DelTpNsundt) ! write(6,8)' TGRAV: delTpNsun ', delTpNsun, delTpNsundt ! write(6,8)'TGRAV: delgrav_Sun ', delgrav_Sun, ddelgrav_Sun ! ! Solar gravitational delay on near-field object: ! term_a = R1J + DOTP(STAR12(1,1),R1Sunt1) ! dterm_a = dR1J + DOTP(STAR12dt(1,1),R1Sunt1) + & ! & DOTP(STAR12(1,1),dR1Sunt1) ! ! term_a1 = R0J + DOTP(STAR12(1,2),R0Sun) ! dterm_a1 = dR0J + -DOTP(STAR12dt(1,2),R0Sun) + & ! & DOTP(STAR12(1,2),dR0Sun) ! ! term_b = R2J + DOTP(STAR12(1,2),R2Sunt1) ! dterm_b = dR2J + DOTP(STAR12dt(1,2),R2Sunt1) + & ! & DOTP(STAR12(1,2),dR2Sunt1) ! ! term_b1 = R0J + DOTP(STAR12(1,1),R0Sun) ! dterm_b1 = dR0J + DOTP(STAR12dt(1,1),R0Sun) + & ! & DOTP(STAR12(1,1),dR0Sun) ! ! delgrav_Sun = (2.D0*GMSUN/VLIGHT3) * & ! & DLOG ( (term_a*term_a1)/(term_b*term_b1)) ! ddelgrav_Sun = (2.D0*GMSUN/VLIGHT3) * (dterm_a/term_a + & ! & dterm_a1/term_a1 - dterm_b/term_b - dterm_b1/term_b1) ! write(6,*) ' Sun: ' ! write(6,*) ' GMSUN,VLIGHT3 ', GMSUN,VLIGHT3 ! write(6,*) ' term_a,dterm_a ', term_a,dterm_a ! write(6,*) ' term_a1,dterm_a1 ', term_a1,dterm_a1 ! write(6,*) ' term_b,dterm_b ', term_b,dterm_b ! write(6,*) ' term_b1,dterm_b1 ', term_b1,dterm_b1 ! write(6,8) 'Sekido: delgrav_Sun ', delgrav_Sun, ddelgrav_Sun ! ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! Moon: del_t_Moon = DOTP(unitK1,x1Moon) / Vlight ! If(del_t_Moon .lt. 0.D0) del_t_Moon = 0.0D0 If(del_t_Moon .gt. T0_T1) del_t_Moon = T0_T1 ! write(6,8)' del_t_Moon ',del_t_Moon ! Do l=1,3 ! SSBC vector to Moon at time of closest approach and its time derivative: XMoont1J(l) = XsubMoon(l) - VsubMoon(l)*del_t_Moon dXMoont1J(l) = VsubMoon(l) !Derivative (approx. - no acceleration) ! ! Vector from the Moon to receiver 1 R1Moont1(l) = Xsub1(l) - XMoont1J(l) dR1Moont1(l) = dXsub1(l) - dXMoont1J(l) !Derivative ! ! Vector from the Moon to receiver 2 R2Moont1(l) = Xsub2(l) - VsubEarth(l)* KdotB/VLIGHT - & & XMoont1J(l) dR2Moont1(l) = dXsub2(l) - AsubEarth(l)*KdotB/VLIGHT - & & VsubEarth(l)*dKdotB/VLIGHT - dXMoont1J(l) !Derivative ! ! Vector from Moon to near-field R0moon(l) = R1moont1(l) + R1(l) dR0moon(l) = dR1moont1(l) + R1dt(l) ! Enddo ! R1J = VECMG (R1Moont1) dR1J = DOTP(R1Moont1,dR1Moont1)/R1J R2J = VECMG (R2Moont1) dR2J = DOTP(R2Moont1,dR2Moont1)/R2J R0J = VECMG (R0Moon) dR0J = DOTP(R0Moon,dR0Moon)/R0J ! ! Write(6,*) ' Moon: ' ! Write(6,*) ' R1J, dR1J ', R1J, dR1J ! Write(6,*) ' R2J, dR2J ', R2J, dR2J ! Write(6,*) ' R0J, dR0J ', R0J, dR0J ! Write(6,*) ' ' ! ! CALL TGRAV(GMMOON,R0J,R1J,R2J,R10,R20, dR0J,dR1J,dR2J, & & dR10,dR20, delgrav_Moon, ddelgrav_Moon) ! & dR10,dR20, DelTpNmoon,DelTpNmoondt) ! write(6,8)'TGRAV: DelTpNmoon ', delTpNmoon, DelTpNmoondt ! write(6,8)'TGRAV: delgrav_Moon ', delgrav_Moon, ddelgrav_Moon ! ! Lunar gravitational delay on near-field object: ! term_a = R1J + DOTP(STAR12(1,1),R1Moont1) ! dterm_a = dR1J + DOTP(STAR12dt(1,1),R1Moont1) + & ! & DOTP(STAR12(1,1),dR1Moont1) ! ! term_a1 = R0J + DOTP(STAR12(1,2),R0Moon) ! dterm_a1 = dR0J + DOTP(STAR12dt(1,2),R0Moon) + & ! & DOTP(STAR12(1,2),dR0Moon) ! ! term_b = R2J + DOTP(STAR12(1,2),R2Moont1) ! dterm_b = dR2J + DOTP(STAR12dt(1,2),R2Moont1) + & ! & DOTP(STAR12(1,2),dR2Moont1) ! ! term_b1 = R0J + DOTP(STAR12(1,1),R0Moon) ! dterm_b1 = dR0J + DOTP(STAR12dt(1,1),R0Moon) + & ! & DOTP(STAR12(1,1),dR0Moon) ! ! delgrav_Moon = (2.D0*GMMOON/VLIGHT3) * & ! & DLOG ( (term_a*term_a1)/(term_b*term_b1) ) ! ddelgrav_Moon = (2.D0*GMMOON/VLIGHT3) * (dterm_a/term_a + & ! & dterm_a1/term_a1 - dterm_b/term_b - dterm_b1/term_b1) ! write(6,*) ' Moon: ' ! write(6,*) ' term_a,dterm_a ', term_a,dterm_a ! write(6,*) ' term_a1,dterm_a1 ', term_a1,dterm_a1 ! write(6,*) ' term_b,dterm_b ', term_b,dterm_b ! write(6,*) ' term_b1,dterm_b1 ', term_b1,dterm_b1 ! write(6,8) 'Sekido: delgrav_Moon ', delgrav_Moon, ddelgrav_Moon ! ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! Planets: Do k=1,7 !Planet loop del_t_Planet(k) = DOTP(unitK1,x1Planet(1,k)) / Vlight ! If(del_t_Planet(k) .lt. 0.D0) del_t_Planet(k) = 0.0D0 If(del_t_Planet(k) .gt. T0_T1) del_t_Planet(k) = T0_T1 ! write(6,'("del_t_Planet(",i1,") = ",d25.16)')k,del_t_Planet(k) ! Do l=1,3 ! SSBC vector to Planet at time of closest approach and its time derivative: XPlant1J(l,k) = SPLANET(l,1,k) - SPLANET(l,2,k)*del_t_Planet(k) dXPlant1J(l,k) = SPLANET(l,2,k) ! Derivative (approximate) ! ! Vector from Planet to receiver 1 R1Plant1(l,k) = Xsub1(l) - XPlant1J(l,k) dR1Plant1(l,k) = dXsub1(l) - dXPlant1J(l,k) !Derivative ! ! Vector from Planet to receiver 2 R2Plant1(l,k) = Xsub2(l) - VsubEarth(l)*KdotB/VLIGHT - & & XPlant1J(l,k) dR2Plant1(l,k) = dXsub2(l) - AsubEarth(l)*KdotB/VLIGHT - & & VsubEarth(l)*dKdotB/VLIGHT - dXPlant1J(l,k) !Derivative ! ! Vector from Planet k to near-field object R0plan(l,k) = R1plant1(l,k) + R1(l) dR0plan(l,k) = dR1plant1(l,k) + R1dt(l) ! Enddo ! R1J = VECMG (R1Plant1(1,k)) dR1J = DOTP(R1Plant1(1,k),dR1Plant1(1,k))/R1J R2J = VECMG (R2Plant1(1,k)) dR2J = DOTP(R2Plant1(1,k),dR2Plant1(1,k))/R2J R0J = VECMG (R0plan(1,k)) dR0J = DOTP(R0plan(1,k),dR0plan(1,k))/R0J GMp = GMPlanet(k) ! ! Write(6,*) ' Planet ', k ! Write(6,*) ' R1J, dR1J ', R1J, dR1J ! Write(6,*) ' R2J, dR2J ', R2J, dR2J ! Write(6,*) ' R0J, dR0J ', R0J, dR0J ! Write(6,*) ' GMp ', GMp ! Write(6,*) ' ' ! CALL TGRAV (GMp,R0J,R1J,R2J,R10,R20, dR0J,dR1J,dR2J, & & dR10,dR20, DelTpNPl,DelTpNPldt) ! DelTpNPlan(k) = DelTpNPl ! DelTpNPlandt(k) = DelTpNPldt ! write(6,*)' #, delTpNPLAN(k) ', k, delTpNPlan(k), DelTpNPlandt(k) delgrav_Plan(k) = DelTpNPl ddelgrav_Plan(k) = DelTpNPldt ! write(6,*)'TGRAV: delgrav_Plan ', k, delgrav_Plan(k), ddelgrav_Plan(k) ! ! Planet #k gravitational delay on near-field object: ! term_a = R1J + DOTP(STAR12(1,1),R1Plant1(1,k)) ! dterm_a = dR1J + DOTP(STAR12dt(1,1),R1Plant1(1,k)) + & ! & DOTP(STAR12(1,1),dR1Plant1(1,k)) ! ! term_a1 = R0J + DOTP(STAR12(1,2),R0Plan(1,k)) ! dterm_a1 = dR0J + DOTP(STAR12dt(1,2),R0Plan(1,k)) + & ! & DOTP(STAR12(1,2),dR0Plan(1,k)) ! ! term_b = R2J + DOTP(STAR12(1,2),R2Plant1(1,k)) ! dterm_b = dR2J + DOTP(STAR12dt(1,2),R2Plant1(1,k)) + & ! & DOTP(STAR12(1,2),dR2Plant1(1,k)) ! ! term_b1 = R0J + DOTP(STAR12(1,1),R0Plan(1,k)) ! dterm_b1 = dR0J + DOTP(STAR12dt(1,1),R0Plan(1,k)) + & ! & DOTP(STAR12(1,1),dR0Plan(1,k)) ! ! delgrav_Plan(k) = (2.D0*GMPlanet(k)/VLIGHT3) * & ! & DLOG ( (term_a*term_a1)/(term_b*term_b1) ) ! ddelgrav_Plan(k) = (2.D0*GMPlanet(k)/VLIGHT3) * (dterm_a/term_a & ! & + dterm_a1/term_a1 - dterm_b/term_b - dterm_b1/term_b1) ! write(6,*) ' Planet ', k ! write(6,*) ' term_a,dterm_a ', term_a,dterm_a ! write(6,*) ' term_a1,dterm_a1 ', term_a1,dterm_a1 ! write(6,*) ' term_b,dterm_b ', term_b,dterm_b ! write(6,*) ' term_b1,dterm_b1 ', term_b1,dterm_b1 ! write(6,*) ' k,delgrav_Plan ', k,delgrav_Plan(k), ddelgrav_Plan(k) ! write(6,*) ' ' ! Enddo !Planet loop ! !******************* Debug ******************** ! IF(KTHED .ne. 0) Then ! Debug ! write(6,'(/,15x,"Step 2 dump:")') ! write(6,8)' x1Sun ',x1Sun ! write(6,8)' x1Moon ',x1Moon ! Do k=1,7 ! write(6,'("x1Planet(",i1,") = ",3d25.16)') k,x1Planet(1,k), & ! & x1Planet(2,k),x1Planet(3,k) ! write(6,'("SPLANET(",i1,") = ",3d25.16)') k,SPLANET(1,1,k), & ! & SPLANET(2,1,k),SPLANET(3,1,k) ! Enddo ! ! write(6,8)' KdotB, dKdotB ', KdotB, dKdotB ! write(6,8)' del_t_Sun ',del_t_Sun ! write(6,8)' XSunt1J ',XSunt1J ! write(6,8)' dXSunt1J ',dXSunt1J ! write(6,8)' R1Sunt1 ',R1Sunt1 ! write(6,8)' dR1Sunt1 ',dR1Sunt1 ! write(6,8)' R2Sunt1 ',R2Sunt1 ! write(6,8)' dR2Sunt1 ',dR2Sunt1 ! write(6,8)' del_t_Moon ',del_t_Moon ! write(6,8)' XMoont1J ',XMoont1J ! write(6,8)' dXMoont1J ',dXMoont1J ! write(6,8)' R1Moont1 ',R1Moont1 ! write(6,8)' dR1Moont1 ',dR1Moont1 ! write(6,8)' R2Moont1 ',R2Moont1 ! write(6,8)' dR2Moont1 ',dR2Moont1 ! write(6,8)' Vmag_S ', Vmag_S ! write(6,8)' N_hat ', N_hat ! write(6,8)' dN_hat ', dN_hat ! ! Do k=1,7 ! write(6,'("del_t_Planet(",i1,") = ",d25.16)')k,del_t_Planet(k) ! write(6,'("XPlant1J(",i1,") = ",3d25.16)') k,XPlant1J(1,k), & ! & XPlant1J(2,k),XPlant1J(3,k) ! write(6,'("dXPlant1J(",i1,") = ",3d25.16)') k,dXPlant1J(1,k), & ! & dXPlant1J(2,k),dXPlant1J(3,k) ! write(6,'("R1Plant1(",i1,") = ",3d25.16)') k,R1Plant1(1,k), & ! & R1Plant1(2,k),R1Plant1(3,k) ! write(6,'("dR1Plant1(",i1,") = ",3d25.16)') k,dR1Plant1(1,k), & ! & dR1Plant1(2,k),dR1Plant1(3,k) ! write(6,'("R2Plant1(",i1,") = ",3d25.16)') k,R2Plant1(1,k), & ! & R2Plant1(2,k),R2Plant1(3,k) ! write(6,'("dR2Plant1(",i1,") = ",3d25.16)') k,dR2Plant1(1,k), & ! & dR2Plant1(2,k),dR2Plant1(3,k) ! Enddo ! ! Endif ! Debug ! !****************************************************************************** ! Step 4: Find the differential delay due to the Earth. ! C_Earth = (1.0D0 + gamma) * GMEarth/VLIGHT3 vecmg0 = VECMG(SPpxyz) dvecmg0 = DOTP(SPvxyz,SPpxyz)/vecmg0 ! IF(Nzero .ne. 1) THEN vecmg1 = VECMG(x_sub1t1) term_g = vecmg1 + DOTP(STAR12(1,1),x_sub1t1) dterm_g = DOTP(x_sub1t1,w_sub1)/vecmg1 + & & DOTP(STAR12(1,1),w_sub1) + DOTP(STAR12dt(1,1),x_sub1t1) term_g1 = vecmg0 + DOTP(STAR12(1,1),SPpxyz ) dterm_g1 = dvecmg0 + DOTP(STAR12(1,1),SPvxyz) & & + DOTP(STAR12dt(1,1),SPpxyz) ELSE vecmg1 = REARTH term_g = 2.0D0 * REARTH dterm_g = 0.0D0 term_g1 = 2.0D0 * vecmg0 dterm_g1 = 2.0D0 * dvecmg0 ENDIF ! IF(Nzero .ne. 2) THEN vecmg2 = VECMG(x_sub2t1) term_h = vecmg2 + DOTP(STAR12(1,2),x_sub2t1) dterm_h = DOTP(x_sub2t1,w_sub2)/vecmg2 + & & DOTP(STAR12(1,2),w_sub2) + DOTP(STAR12dt(1,2),x_sub2t1) term_h1 = vecmg0 + DOTP(STAR12(1,2),SPpxyz) dterm_h1 = dvecmg0 + DOTP(STAR12(1,2),SPvxyz) & & + DOTP(STAR12dt(1,2),SPpxyz) ELSE vecmg2 = REARTH term_h = 2.0D0 * REARTH dterm_h = 0.0D0 term_h1 = 2.0D0 * vecmg0 dterm_h1 = 2.0D0 * dvecmg0 ENDIF ! delgrav_Earth = C_Earth * DLOG( (term_g*term_h1) / & & (term_h*term_g1) ) ddelgrav_Earth = C_Earth * ( dterm_g/term_g + & & dterm_h1/term_h1 - dterm_h/term_h - dterm_g1/term_g1 ) ! ! write(6,8)' vecmg0, dvecmg0 ', vecmg0, dvecmg0 ! write(6,8)' vecmg1, vecmg2 ', vecmg1, vecmg2 ! write(6,8)' term_g, dterm_g ', term_g, dterm_g ! write(6,8)' term_g1,dterm_g1 ', term_g1, dterm_g1 ! write(6,8)' term_h, dterm_h ', term_h, dterm_h ! write(6,8)' term_h1,dterm_h1 ', term_h1, dterm_h1 ! write(6,8)'Sekido: delgrav_Earth, ddelgrav_Earth ', & ! & delgrav_Earth,ddelgrav_Earth ! !* IF(KTHED .ne. 0) Then ! Debug !* write(6,'(/,15x,"Step 4 dump:")') !* write(6,8)' C_Earth ',C_Earth !* write(6,8)' term_g, dterm_g ', term_g, dterm_g !* write(6,8)' term_h, dterm_h ', term_h, dterm_h !* write(6,8)' delgrav_Earth, ddelgrav_Earth ', & !* & delgrav_Earth,ddelgrav_Earth !* Endif ! Debug ! !****************************************************************************** ! Step 5: Add up all components to get the total ! differential gravitational delay. ! [Does not include the term for observations close to the Sun.] ! ! Add up the 10 gravitational components delgrav = delgrav_Sun + delgrav_Moon + delgrav_Earth ddelgrav = ddelgrav_Sun + ddelgrav_Moon + ddelgrav_Earth Do k=1,7 delgrav = delgrav + delgrav_Plan(k) ddelgrav = ddelgrav + ddelgrav_Plan(k) Enddo ! write(6,8)' Sekido: delgrav,ddelgrav ', delgrav,ddelgrav ! !* IF(KTHED .ne. 0) Then ! Debug !* write(6,'(/,15x,"Step 5 dump:")') !* write(6,8)' delgrav_Planets, ddelgrav_Planets ', & !* & delgrav_Planets, ddelgrav_Planets !* write(6,8)' delgrav, ddelgrav ', & !* & delgrav, ddelgrav !* Endif ! Debug ! !****************************************************************************** ! Step 6: Add the total differential gravitational delay to the rest of the ! a priori vacuum delay, Equation 9. ! ! Find Sun's potential: ! U_Sun = GMSUN/VLIGHT2 vecmg1 = VECMG(R_Earth_Sun) U = U_Sun/vecmg1 ! Derivative: dU = -U_Sun * Dotp(R_Earth_Sun,V_Earth_Sun) / vecmg1**3 ! ! Compute individual terms of Eqn. 20 ! term2a = KdotB/VLIGHT dterm2a = dKdotB/VLIGHT ! derivative ! Write(6,8) ' term2a,dterm2a ', term2a,dterm2a ! term2b = 1.D0 - ((1.D0 + gamma) * U) dterm2b = -(1.D0 + gamma) * dU ! derivative ! Write(6,8) ' term2b-1,dterm2b ', (term2b-1.D0),dterm2b ! absVEarth = VECMG(VsubEarth) term2c = (absVEarth)**2 / (2.D0*VLIGHT2) dterm2c = Dotp(VsubEarth,AsubEarth) / VLIGHT2 ! derivative !! term2c = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! dterm2c = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Write(6,8) ' term2c,dterm2c ', term2c,dterm2c ! term2d = DOTP(VsubEarth,w_sub2) / VLIGHT2 dterm2d = (DOTP(AsubEarth,w_sub2) + DOTP(VsubEarth,a_sub2)) & & / VLIGHT2 ! derivative !! term2d = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! dterm2d = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Write(6,8) ' term2d,dterm2d ', term2d,dterm2d ! ! Combine terms 2b,2c,2d term2bcd = term2b - term2c - term2d dterm2bcd = dterm2b - dterm2c - dterm2d ! derivative ! term2 = term2a * term2bcd dterm2 = term2a * dterm2bcd + dterm2a * term2bcd ! derivative ! Write(6,8) ' term2, dterm2 ', term2,dterm2 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! VdotB = DOTP(VsubEart0,b_nought) dVdotB = DOTP(AsubEart0,b_nought) + DOTP(VsubEart0,db_nought) !derivative ! !!!! KdotV = DOTP(unit_K,VsubEarth) !!!! dKdotV = DOTP(unit_K,AsubEarth) ! derivative ! term3a = VdotB/VLIGHT2 dterm3a = dVdotB/VLIGHT2 ! Write(6,8) ' term3a,dterm3a ', term3a,dterm3a ! Call VECAD (VsubEart0, w_sub2, V2 ) Call VECAD (AsubEart0, a_sub2, dV2 ) ! Write(6,8) ' V2, dV2 ', V2, dV2 term3b = DOTP(STAR12(1,2), V2)/VLIGHT ! term3b = DOTP(STAR12(1,2), w_sub2)/VLIGHT !!!!!! TEST !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! dterm3b = (DOTP(STAR12dt(1,2),V2) + DOTP(STAR12(1,2),dV2))/VLIGHT ! dterm3b = (DOTP(STAR12dt(1,2),w_sub2) + DOTP(STAR12(1,2),a_sub2))/VLIGHT !!!!! TEST !!!!!!! !?? term3b = DOTP(STAR12(1,2), w_sub2)/VLIGHT !?? dterm3b = (DOTP(STAR12dt(1,2),w_sub2) + DOTP(STAR12(1,2),a_sub2))/VLIGHT ! Write(6,8) ' STAR12(1,2) ', STAR12(1,2),STAR12(2,2),STAR12(3,2) ! Write(6,8) ' STAR12dt(1,2) ', STAR12dt(1,2),STAR12dt(2,2),STAR12dt(3,2) ! Write(6,8) ' term3b,dterm3b ', term3b,dterm3b ! DO I = 1,3 Vec3c(I) = VsubEart0(I) + 2.D0*w_sub2(I) dVec3c(I) = AsubEart0(I) + 2.D0*a_sub2(I) ENDDO term3c = DOTP(Vec3c,STAR) / (2.D0*VLIGHT) dterm3c = (DOTP(dVec3c,STAR) + DOTP(Vec3c,STARdt))/(2.D0*VLIGHT) ! Write(6,8) ' term3c,dterm3c ', term3c,dterm3c ! term3 = term3a * (1.D0 + term3b - term3c) dterm3 = dterm3a*(1.D0+term3b-term3c) + term3a*(dterm3b-dterm3c) !! term3 = term3a * (1.D0 ) !!!!!!!!!!!!!!!!! !! dterm3 = dterm3a*(1.D0) !!!!!!!!!!!!!!!!! !? term3 = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!! !? dterm3 = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!! ! Write(6,8) ' term3, dterm3 ', term3,dterm3 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! term4a = 1.D0 + DOTP(STAR12(1,2),V2)/VLIGHT ! term4a = 1.D0 + DOTP(STAR12(1,2),w_sub2)/VLIGHT dterm4a = (DOTP(STAR12dt(1,2),V2) + DOTP(STAR12(1,2),dV2))/VLIGHT ! dterm4a = (DOTP(STAR12dt(1,2),w_sub2) + DOTP(STAR12(1,2),w_sub2))/VLIGHT !?? term4a = 1.D0 + DOTP(STAR12(1,2),w_sub2)/VLIGHT !?? dterm4a = (DOTP(STAR12dt(1,2),w_sub2) + DOTP(STAR12(1,2),a_sub2))/VLIGHT ! Write(6,8) ' term4a,dterm4a ', term4a,dterm4a ! Call CROSP(V2,STAR12(1,2),H1) !?? Call CROSP(w_sub2,STAR12(1,2),H1) H2 = (VECMG(H1)/VLIGHT)**2 H3 = KdotB/(2.D0*R2mag) H = H2*H3 Call CROSP(dV2,STAR12(1,2),dH1a) Call CROSP(V2,STAR12dt(1,2),dH1b) !?? Call CROSP(a_sub2,STAR12(1,2),dH1a) !?? Call CROSP(w_sub2,STAR12dt(1,2),dH1b) Call VECAD (dH1a, dH1b, dH1) ! dH2 = 2.D0*(VECMG(H1)*(VECMG(dH1)))/(VLIGHT2) dH2 = (2.D0 * DOTP(H1,dH1)) / VLIGHT2 dH3 = (dKdotB/R2mag - KdotB*R2magdt/R2mag**2)/2.0D0 dH = dH2*H3 + H2*dH3 ! Write(6,8) ' H1, dH1 ', H1, dH1 ! Write(6,8) ' H2, dH2 ', H2, dH2 ! Write(6,8) ' H3, dH3 ', H3, dH3 ! Write(6,8) ' H, dH ', H, dH ! term4 = term4a * (1.D0 + H) dterm4 = dterm4a*(1.D0 + H) + term4a*dH !!? term4 = (1.D0 + H) !!!!!!!!!!!!!!!!!!!!!! !!? dterm4 = 0.D0 !!!!!!!!!!!!!!!!!!!!!! ! Write(6,8) ' term4, dterm4 ', term4,dterm4 ! tv2_tv1 = (-term2 - term3 + delgrav) / term4 dtv2_tv1 = (-dterm2 - dterm3 + ddelgrav) / term4 - & & (1.D0/term4**2) * (-term2 - term3 + delgrav) * dterm4 ! ! write(6,8) 'S: tv2_tv1, dtv2_tv1 ',tv2_tv1, dtv2_tv1 ! Write(6,*) ' No delTgrav (ps): ', (tv2_tv1 - delgrav/term4)*1.D12 ! ! IF(KTHED .ne. 0) Then ! Debug ! write(6,'(/,15x,"Step 6 dump:")') ! write(6,8)' U_Sun, U, dU ',U_Sun, U, dU ! write(6,8)' term1, dterm1 ',term1, dterm1 ! write(6,8)' term2a, dterm2a ',term2a, dterm2a ! write(6,8)' term2b, dterm2b ',term2b, dterm2b ! write(6,8)' absVEarth',absVEarth ! write(6,8)' term2c, dterm2c ',term2c, dterm2c ! write(6,8)' term2d, dterm2d ',term2d, dterm2d ! write(6,8)' term2bcd, dterm2bcd ',term2bcd, dterm2bcd ! write(6,8)' term2, dterm2 ',term2, dterm2 ! write(6,8)' VdotB, dVdotB ',VdotB, dVdotB ! write(6,8)' KdotV, dKdotV ',KdotV, dKdotV ! write(6,8)' term3a, dterm3a ',term3a, dterm3a ! write(6,8)' term3b, dterm3b ',term3b, dterm3b ! write(6,8)' term3, dterm3 ',term3, dterm3 ! write(6,8)' vec_sum, dvec_sum ',vec_sum, dvec_sum ! write(6,8)' term4, dterm4 ',term4, dterm4 ! write(6,8)' term123, dterm123 ',term123, dterm123 ! write(6,8)' tv2_tv1, dtv2_tv1 ',tv2_tv1, dtv2_tv1 ! Endif ! Debug ! !****************************************************************************** ! Step 7: Calculate the aberrated source vectors for use in the tropospheric ! propogation delay calculation. ! ! [Note: These calculations are in the atmosphere module, subroutine ! ATMG, and do not need to be done here. The aberrated source vectors are ! used in the atmosphere module to compute the topocentric azimuths and ! elevations, and in the axis offset module (along with refraction) to ! compute the vector axis offset.] ! ! call vecad(VsubEarth,w_sub1,V_w1) ! K_V_w1 = Dotp(unit_K,V_w1) ! ! call vecad(VsubEarth,w_sub2,V_w2) ! K_V_w2 = Dotp(unit_K,V_w2) ! ! Do l=1,3 ! k_sub1(l) = unit_K(l) + (VsubEarth(l) + w_sub1(l)) / VLIGHT ! . - unit_K(l) * K_V_w1 /VLIGHT ! k_sub2(l) = unit_K(l) + (VsubEarth(l) + w_sub2(l)) / VLIGHT ! . - unit_K(l) * K_V_w2 /VLIGHT ! Enddo ! ! IF(KTHED .ne. 0) Then ! Debug ! write(6,'(/,15x,"Step 7 dump:")') ! write(6,8)' V_w1, K_V_w1 ',V_w1, K_V_w1 ! write(6,8)' V_w2, K_V_w2 ',V_w2, K_V_w2 ! write(6,8)' k_sub1, k_sub2 ',k_sub1, k_sub2 ! Endif ! Debug ! !****************************************************************************** ! Step 8: Add the geometric part of the tropospheric propogation delay to the ! vacuum delay. ! [Geocenter station: DATMC(Nzero,1) and DATMC(Nzero,2) should ! already be zero.] ! call vecsb(w_sub2,w_sub1,w2_w1) call vecsb(a_sub2,a_sub1,a2_a1) K_dotw2w1 = DOTP(STAR,w2_w1) + DOTP(STARdt,b_nought) dK_dotw2w1 = DOTP(STARdt,w2_w1) + DOTP(STAR,a2_a1) + & & DOTP(STARdt,db_nought) !!! & + DOTP(STARdtdt,b_nought) ! ! Negative sign because station 1 atmosphere is negative. tg2_tg1 = tv2_tv1 - DATMC(1,1) * K_dotw2w1/VLIGHT dtg2_tg1 = dtv2_tv1 - DATMC(1,2) * K_dotw2w1/VLIGHT & & - DATMC(1,1) * dK_dotw2w1/VLIGHT ! write(6,8)' w2_w1, a2_a1 ', w2_w1, a2_a1 ! write(6,8)' DOTP(STAR,w2_w1)/c ', DOTP(STAR,w2_w1)/VLIGHT ! write(6,8)' DOTP(STARdt,b_nought)/c ', DOTP(STARdt,b_nought)/VLIGHT ! write(6,8)' K_dotw2w1, dK_dotw2w1 ', K_dotw2w1/VLIGHT, dK_dotw2w1/VLIGHT !* write(6,8)' DATMC ', DATMC ! write(6,8)' DATMC*K_dotw2w1/c ', DATMC(1,1)*K_dotw2w1/VLIGHT ! write(6,8)'S: tg2_tg1, dtg2_tg1 ', tg2_tg1, dtg2_tg1 !* write(6,8)' tg2_tg1,dtg2_tg1+atm:', tg2_tg1+DATMC(1,1)+DATMC(2,1), & !* & dtg2_tg1+DATMC(1,2)+DATMC(2,2) ! ! ! IF(KTHED .ne. 0) Then ! Debug ! write(6,'(/,15x,"Step 8 dump:")') ! write(6,8)' w2_w1, a2_a1 ',w2_w1, a2_a1 ! write(6,8)' K_dotw2w1, dK_dotw2w1 ', K_dotw2w1, dK_dotw2w1 ! write(6,8)' DATMC ', DATMC ! write(6,8)' tg2_tg1, dtg2_tg1 ', tg2_tg1, dtg2_tg1 ! Endif ! Debug ! !****************************************************************************** ! Step 9: Step 9 is to compute the total delay by adding in the "Best" ! estimate of the troposphere propogation delay for each site. ! Normally this will be done 'on the fly' in program SOLVE. ! However, we add the Niell atmosphere term here if the atmosphere ! flag, KATMC is 1. The default is not to add it. ! [Geocenter station: DATMC(Nzero,1) and DATMC(Nzero,2) should ! already be zero.] ! ! IF(KATMC.eq.1) then ! tg2_tg1 = tg2_tg1 + DATMC(1,1) + DATMC(2,1) ! dtg2_tg1 = dtg2_tg1 + DATMC(1,2) + DATMC(2,2) ! ! IF(KTHED .ne. 0) Then ! Debug ! write(6,'(/,15x,"Step 9 dump:")') ! write(6,8)'(KATMC=1:) tg2_tg1, dtg2_tg1 ', tg2_tg1, dtg2_tg1 ! Endif ! Debug ! ! Endif ! ! 8 FORMAT(A,4D25.16/(7X,5D25.16)) ! RETURN END ! !********************************************************************** SUBROUTINE TGRAV(GM, R0J,R1J,R2J,R10,R20, R0Jdt,R1Jdt,R2Jdt, & & R10dt,R20dt, DelTpN,DelTpNdt) IMPLICIT NONE ! ! Subroutine to compute the gravitational deflection for a solar ! system body (Sun, Moon, planets) on a near-field object. ! ! ! 1. GM - Gravitational constant times the mass of the gravitating ! body. Real*8 GM, R0J, R1J, R2J, R10, R20, R0Jdt, R1Jdt, R2Jdt, R10dt, & & R20dt, DelTpN,DelTpNdt Real*8 term1a, term1b, term2a, term2b, term1adt, term1bdt, & & term2adt, term2bdt, term12, term12dt !* INCLUDE 'cphys11.i' ! Variables 'from': ! 1. VLIGHT - The velocity of light in vacuum. (m/sec) ! 2. VLIGHT2 - The velocity of light squared. (m/sec)**2 ! 3. VLIGHT3 - The velocity of light cubed. (m/sec)**3 ! ! term1a = R2J + R0J + R20 term1adt = R2Jdt + R0Jdt + R20dt ! term1b = R2J + R0J - R20 term1bdt = R2Jdt + R0Jdt - R20dt ! term2a = R1J + R0J - R10 term2adt = R1Jdt + R0Jdt - R10dt ! term2b = R1J + R0J + R10 term2bdt = R1Jdt + R0Jdt + R10dt ! term12 = (term1a/term1b) * (term2a/term2b) term12dt = (term1adt/term1b) * (term2a/term2b) - & & (term1a*term1bdt)/(term1b**2) * (term2a/term2b) + & & (term1a/term1b) * (term2adt/term2b) - & & (term1a/term1b) * (term2a*term2bdt)/(term2b**2) ! ! DelTpN = (2.D0*GM/VLIGHT3) * DLOG ( ((R2J+R0J+R20)/(R2J+R0J-R20)) & ! * ((R1J+R0J+R10)/(R1J+R0J-R10)) ) ! DelTpN = (2.D0*GM/VLIGHT3) * DLOG (term12) DelTpNdt = (2.D0*GM/VLIGHT3) * (1.D0/term12)*term12dt ! RETURN END ! !********************************************************************** Subroutine DUEV (DATMC, SITEP, SITEV, SITEA, SUN, T0_T1, & & EARTH, XMOON, UTC, TT, TDB, TDBg, td2_td1, dtd2_td1) ! IMPLICIT NONE REAL*8 DATMC(2,2), SUN(3,2), EARTH(3,3), SITEA(3,2), SITEP(3,2), & & SITEV(3,2), DOTP, VECMG, UTC, TT, TDB, TDBg, XMOON(3,2) Real*8 T0_T1, td2_td1, dtd2_td1 Real*8 U, dU, U_Sun, Vecmg1, F1, VdotR1, VdotR2, VdotNFO, & & del_T0, gm1, RLT_01, R_01mag, R0sunT0mag, R1sunT1mag, & & T1_T0, ddel_T0, R01sunmag Real*8 VsubEarth(3), R_Earth_Sun(3), V_Earth_Sun(3), R1_T1(3), & & R2_T1(3), R0_T0(3), dR0_T0(3), R_01(3), delSun0(3), & & SUN_0(3), R0sunT0(3), R1sunT1(3), Rb_01(3), dR0(3), & & R1(3), EARTH_0(3), RsunT0(3), RsunT1(3), R01sun(3) Real*8 R_02(3), Site2(3), R2_T2(3), SunT2(3), EarthT2(3), & & RsunT2(3), R2sunT2(3), R02sun(3) Real*8 T2_T0, del_T2, R2_T2mag, R_02mag, R2sunT2mag, R02sunmag, & & RLT_02, ddel_T2 Real*8 T2_T1, VEsq, UE, B(3), VEdotB, VEdotV2 Real*8 delMoon0(3),MOON_0(3),RmoonT0(3),R0moonT0(3),R0moonT0mag, & & RmoonT1(3),R1moonT1(3),R1moonT1mag,R01moon(3),R01moonmag, & & RLT_01s,RLT_01m,RLT_01p,delpl(3),Spl_0(3),R0planT0(3), & & R0planT0mag,R1planT1(3),R1planT1mag,R01plan(3),R01planmag, & & RLTp(7) Real*8 MoonT2(3),RmoonT2(3),R2moonT2(3),R2moonT2mag,R02moon(3), & & R02moonmag,RLT_02s,RLT_02m,RLT_02p,Spl_2(3),R2planT2(3), & & R2planT2mag,R02plan(3),R02planmag Real*8 delEarth0(3), REarthT0(3), R0earthT0(3), R0earthT0mag, & & RearthT1(3), R1earthT1(3), R1earthT1mag, R01earth(3), & & R01earthmag, RLT_01e, RearthT2(3), R2earthT2(3), & & R2earthT2mag, R02earth(3), R02earthmag, RLT_02e, & & delearth2(3), delSun2(3), delmoon2(3) Integer*4 I, J, K ! INCLUDE 'd_input.i' ! Variables from: ! 1. SPpxyz(3) - Geocentric satellite position at T0 (meters). ! INCLUDE 'cphys11.i' ! Variables 'from': ! 1. VLIGHT - The velocity of light in vacuum. (m/sec) ! 2. VLIGHT2 - The velocity of light squared. (m/sec)**2 ! 3. VLIGHT3 - THE VELOCITY OF LIGHT IN VACUUM CUBED. ! (M**3/SEC**3) ! 4. GMSUN, GMMOON, GMEARTH, GMPLANET(7) - Gravitational constant ! times the masses of the Sun, Moon, Earth, and the ! other planets except Pluto. ! INCLUDE 'csolsys11.i' ! Variables 'from': ! 1. SPLANET(3,2,7) - The J2000.0 Solar System Barycentric positions ! and velocities of all planets except the Earth ! and Pluto. (meters, meters/sec) The first index ! runs over X,Y, and Z, the second runs over ! position and velocity, and the third runs over ! the planets, where ! 1 = Mercury ! 2 = Venus ! 3 = Mars ! 4 = Jupiter ! 5 = Saturn ! 6 = Uranus ! 7 = Neptune ! ! 2. GPLANET(3,2,7) - The J2000.0 Geocentric positions and velocities ! of all planets except the Earth and Pluto. ! (meters, meters/sec) The first index runs over ! X,Y, and Z, the second runs over position and ! velocity, and the third runs over the planets, ! where ! 1 = Mercury ! 2 = Venus ! 3 = Mars ! 4 = Jupiter ! 5 = Saturn ! 6 = Uranus ! 7 = Neptune ! ! 3. SUNb(3,2) - J2000 Barycentric Sun position and velocity. ! 4. MOONb(3,2) - J2000 Barycentric Moon position and velocity. ! ! Convert spacecraft and site vectors from TT to TDB Do I = 1, 3 VsubEarth(I) = EARTH(I,2) ! SSBC Earth Velocity R_Earth_Sun(I) = -SUN(I,1) ! Vector from Sun to Earth V_Earth_Sun(I) = -SUN(I,2) ! Velocity of Earth w.r.t. the Sun Enddo ! ! Sun's potential: U_Sun = GMSUN/VLIGHT2 vecmg1 = VECMG(R_Earth_Sun) U = U_Sun/vecmg1 ! Derivative: ! dU = -U_Sun * Dotp(R_Earth_Sun,V_Earth_Sun) / vecmg1**3 F1 = 1.D0 - U - 1.48082686741D-8 ! VdotR1 = DOTP(VsubEarth,SITEP(1,1)) ! VdotR2 = DOTP(VsubEarth,SITEP(1,2)) VdotNFO = DOTP(VsubEarth,SPpxyz) ! ! Barycentric TDB vectors: (Duev equation 4) Do I = 1,3 R1_T1(I) = SITEP(I,1)*F1 - 0.5D0*(VdotR1/VLIGHT2)*VsubEarth(I) & & + EARTH(I,1) R2_T1(I) = SITEP(I,2)*F1 - 0.5D0*(VdotR2/VLIGHT2)*VsubEarth(I) & & + EARTH(I,1) Enddo ! write(6,*) 'site1,R1_T1: ', SITEP(1,1),SITEP(2,1),SITEP(3,1),R1_T1 ! write(6,*) 'site2,R2_T1: ', SITEP(1,2),SITEP(2,2),SITEP(3,2),R2_T1 ! ! Initial estimate for (T1 - T0) CALL VECSB (SPpxyz, SITEP(1,1), R_01) R_01mag = VECMG(R_01) del_T0 = VECMG(R_01)/VLIGHT ! write(6,*) 'DUEV: T0_T1, Initial del_T0 ', T0_T1, del_T0 ! gm1 = 2.D0*GMSUN/VLIGHT2 ! ! Three iterations should be enough DO J = 1,3 Do I = 1,3 ! Barycentric at T1 - del_T0 R0_T0(I) = SPpxyz(I)*F1 - 0.5D0*(VdotNFO/VLIGHT2)*VsubEarth(I) & & + (EARTH(I,1) - del_T0*EARTH(I,2)) - & & SPvxyz(I)*(T0_T1 - del_T0) Enddo ! ! Barycentric vectors Call VECSB(R1_T1, R0_T0, R_01) ! equation 15 R_01mag = VECMG(R_01) !!! write(6,*) 'R_01,R_01mag: ', R_01,R_01mag ! CALL VECMU(SUNb(1,2),del_T0,delSun0) CALL VECSB(SUNb(1,1), delSun0, RsunT0) ! write(6,*) 'new RsunT0: ', RsunT0 CALL VECSB (R0_T0, RsunT0, R0sunT0) ! eqn 16, Sun to spacecraft at T0 R0sunT0mag = VECMG(R0sunT0) !!! write(6,*) 'R0sunT0,R0sunT0mag: ', R0sunT0,R0sunT0mag ! CALL VECMU(MOONb(1,2),del_T0,delMoon0) CALL VECSB(MOONb(1,1), delMoon0, RmoonT0) ! write(6,*) 'new RmoonT0: ', RmoonT0 CALL VECSB (R0_T0, RmoonT0, R0moonT0) ! eqn 16, Moon to spacecraft at T0 R0moonT0mag = VECMG(R0moonT0) ! CALL VECMU(EARTH(1,2),del_T0,delEarth0) CALL VECSB(EARTH(1,1), delMoon0, REarthT0) CALL VECSB (R0_T0, RearthT0, R0earthT0) ! eqn 16, Earth to spacecraft at T0 R0earthT0mag = VECMG(R0earthT0) ! write(6,*) 'RearthT0: ', RearthT0 ! write(6,*) 'R0earthT0: ', R0earthT0 ! write(6,*) 'R0earthT0mag: ', R0earthT0mag CALL VECEQ(EARTH(1,1),RearthT1) ! write(6,*) 'RearthT1: ', RearthT1 CALL VECSB(R1_T1, RearthT1, R1earthT1) ! eqn 16, Earth to site 1 at T1 R1earthT1mag = VECMG(R1earthT1) ! write(6,*) 'R1earthT1: ', R1earthT1 ! write(6,*) 'R1earthT1mag: ', R1earthT1mag CALL VECSB (R1earthT1,R0earthT0,R01earth) ! eqn 17, Earth R01earthmag = VECMG(R01earth) ! write(6,*) 'R01earth: ', R01earth ! write(6,*) 'R01earthmag: ', R01earthmag ! CALL VECEQ(SUNb(1,1),RsunT1) ! write(6,*) 'new RsunT1: ', RsunT1 CALL VECSB(R1_T1, RsunT1, R1sunT1) ! eqn 16, Sun to site 1 at T1 R1sunT1mag = VECMG(R1sunT1) !!! write(6,*) 'R1sunT1,R1sunT1mag: ', R1sunT1,R1sunT1mag ! CALL VECEQ(MOONb(1,1),RmoonT1) ! write(6,*) 'new RmoonT1: ', RmoonT1 CALL VECSB(R1_T1, RmoonT1, R1moonT1) ! eqn 16, Moon to site 1 at T1 R1moonT1mag = VECMG(R1moonT1) ! CALL VECSB (R1sunT1,R0sunT0,R01sun) ! eqn 17, Sun R01sunmag = VECMG(R01sun) !!! write(6,*) 'R01sun,R01sunmag: ', R01sun,R01sunmag CALL VECSB (R1moonT1,R0moonT0,R01moon) ! eqn 17, Moon R01moonmag = VECMG(R01moon) ! Sun: RLT_01s = (2.D0*GMSUN/VLIGHT3) * DLOG( (R0sunT0mag + R1sunT1mag & & + R01sunmag + gm1)/(R0sunT0mag + R1sunT1mag - R01sunmag + gm1) ) ! Moon: RLT_01m = (2.D0*GMMOON/VLIGHT3) * DLOG( (R0moonT0mag + R1moonT1mag & & + R01moonmag)/(R0moonT0mag + R1moonT1mag - R01moonmag) ) ! write(6,*) 'RLT_01s, RLT_01m ', RLT_01s,RLT_01m ! ! Earth: (Contribution is zero at the geocenter) ! RLT_01e = (2.D0*GMEarth/VLIGHT3) * DLOG( (R0earthT0mag + R1earthT1mag & ! & + R01earthmag)/(R0earthT0mag + R1earthT1mag - R01earthmag) ) ! write(6,*) 'RLT_01s ', RLT_01s ! write(6,*) 'RLT_01m ', RLT_01m ! write(6,*) 'RLT_01e ', RLT_01e ! ! Planets: RLT_01p = 0.0D0 ! write(6,*) 'J, del_T0', J,del_T0 Do K = 1,7 Call VECMU(SPLANET(1,2,K),del_T0,delpl) CALL VECSB(SPLANET(1,1,K), delpl, Spl_0) ! SSBC planet vector at T0 CALL VECSB(R0_T0, Spl_0, R0planT0) ! Planet to spacecraft at T0 R0planT0mag = VECMG(R0planT0) ! CALL VECSB(R1_T1,SPLANET(1,1,K), R1planT1) ! Planet to site 1 at T1 R1planT1mag = VECMG(R1planT1) ! CALL VECSB(R1planT1, R0planT0, R01plan) ! ~Planet-to-site1 R01planmag = VECMG(R01plan) ! RLTp(K) = (2.D0*GMPLANET(K)/VLIGHT3) * & & DLOG( (R0planT0mag + R1planT1mag + R01planmag) / & & (R0planT0mag + R1planT1mag - R01planmag) ) RLT_01p = RLT_01p + RLTp(K) Enddo ! RLT_01 = RLT_01s + RLT_01m + RLT_01p T1_T0 = R_01mag/VLIGHT + RLT_01 ddel_T0 = T1_T0 - del_T0 del_T0 = T1_T0 ! write(6,*) 'J,RLT_01, del_T0: ', J, RLT_01, del_T0 ! Enddo ! ! Now find (T2 - T0): ! Initial estimate for (T2 - T0) Call VECSB(R2_T1, R0_T0, R_02) T2_T0 = VECMG(R_02)/VLIGHT del_T2 = T1_T0 - T2_T0 ! difference from T1 ! write(6,*) 'R_02,T2_T0,del_T2: ', R_02,T2_T0,del_T2 ! DO J = 1,3 ! Three iterations ???? ! Barycentric TDB vectors: Do I = 1,3 Site2(I) = SITEP(I,2) - SITEV(I,2)*del_T2 ! Site 2 at T2 VdotR2 = DOTP(VsubEarth,Site2) R2_T2(I) = Site2(I)*F1 - 0.5D0*(VdotR2/VLIGHT2)*VsubEarth(I) & & + (EARTH(I,1) - del_T2*EARTH(I,2)) ! SSBC Site2 vector @ T2 Enddo R2_T2mag = VECMG(R2_T2) ! write(6,*) 'SITEP2,Site2:', SITEP(1,2),SITEP(2,2),SITEP(3,2),Site2 ! write(6,*) 'R2_T2,R2_T2mag:', R2_T2,R2_T2mag ! Call VECSB(R2_T2, R0_T0, R_02) ! eqn 15; Site2 to NFO R_02mag = VECMG(R_02) ! CALL VECMU(SUNb(1,2),del_T2,delSun2) CALL VECSB(SUNb(1,1), delSun2, RsunT2) CALL VECSB (R2_T2, RsunT2, R2sunT2) ! eqn 16, Sun to site 2 at T2 R2sunT2mag = VECMG(R2sunT2) ! write(6,*) 'new RsunT2: ', RsunT2 ! write(6,*) 'new R2sunT2mag: ', R2sunT2mag CALL VECSB (R2sunT2,R0sunT0,R02sun) ! eqn 17, Sun R02sunmag = VECMG(R02sun) ! write(6,*) 'new R02sun,R02sunmag: ', R02sun,R02sunmag ! CALL VECMU(MOONb(1,2),del_T2,delmoon2) CALL VECSB(MOONb(1,1), delmoon2, RmoonT2) CALL VECSB (R2_T2, RmoonT2, R2moonT2) ! eqn 16, Moon to site 2 at T2 R2moonT2mag = VECMG(R2moonT2) ! write(6,*) 'new R2moonT2: ', R2moonT2 ! write(6,*) 'new R2moonT2mag: ', R2moonT2mag CALL VECSB (R2moonT2,R0moonT0,R02moon) ! eqn 17, Moon R02moonmag = VECMG(R02moon) ! write(6,*) 'new R02moon,R02moonmag: ', R02moon,R02moonmag ! CALL VECMU(EARTH(1,2),del_T2,delearth2) CALL VECSB(EARTH(1,1), delearth2, RearthT2) !! write(6,*) 'new RearthT2: ', earthRT2 CALL VECSB (R2_T2, RearthT2, R2earthT2) ! eqn 16, Earth to site 2 at T2 R2earthT2mag = VECMG(R2earthT2) ! write(6,*) 'R2earthT2mag: ', R2earthT2mag CALL VECSB (R2earthT2,R0earthT0,R02earth) ! eqn 17, Earth R02earthmag = VECMG(R02earth) !! write(6,*) 'R02earth: ', R02earth ! write(6,*) 'R02earthmag: ', R02earthmag ! RLT_02s = (2.D0*GMSUN/VLIGHT3) * DLOG((R0sunT0mag + R2sunT2mag & & + R02sunmag + gm1)/(R0sunT0mag + R2sunT2mag - R02sunmag + gm1) ) RLT_02m = (2.D0*GMMOON/VLIGHT3) * DLOG((R0moonT0mag + R2moonT2mag & & + R02moonmag)/(R0moonT0mag + R2moonT2mag - R02moonmag) ) ! write(6,*) 'RLT_02s, RLT_02m', RLT_02s, RLT_02m RLT_02e = (2.D0*GMEarth/VLIGHT3) * DLOG( (R0earthT0mag + R2earthT2mag & & + R02earthmag)/(R0earthT0mag + R2earthT2mag - R02earthmag) ) ! write(6,*) 'RLT_02s ', RLT_02s ! write(6,*) 'RLT_02m ', RLT_02m ! write(6,*) 'RLT_02e ', RLT_02e ! ! Planets: RLT_02p = 0.0D0 Do K = 1,7 Call VECMU(SPLANET(1,2,K),del_T0,delpl) CALL VECSB(SPLANET(1,1,K), delpl, Spl_0) ! SSBC planet vector at T0 CALL VECSB(R0_T0, Spl_0, R0planT0) ! Planet to spacecraft at T0 R0planT0mag = VECMG(R0planT0) ! Call VECMU(SPLANET(1,2,K),del_T2,delpl) CALL VECSB(SPLANET(1,1,K), delpl, Spl_2) ! SSBC planet vector at T2 CALL VECSB(R2_T2,Spl_2, R2planT2) ! Planet to site 2 at T2 R2planT2mag = VECMG(R2planT2) ! CALL VECSB(R2planT2, R0planT0, R02plan) ! ~Planet-to-site2 R02planmag = VECMG(R02plan) ! RLTp(K) = (2.D0*GMPLANET(K)/VLIGHT3) * & & DLOG( (R0planT0mag + R2planT2mag + R02planmag) / & & (R0planT0mag + R2planT2mag - R02planmag) ) RLT_02p = RLT_02p + RLTp(K) Enddo ! write(6,*) 'RLT_02p ', RLT_02p ! RLT_02 = RLT_02s + RLT_02m + RLT_02e + RLT_02p ! write(6,*) 'RLT_02 ', RLT_02 T2_T0 = R_02mag/VLIGHT + RLT_02 ddel_T2 = T2_T0 - (T1_T0 - del_T2) del_T2 = T1_T0 - T2_T0 ! ! write (6,*) 'del_T2, T2-T1: ', del_T2, (T2_T0 - T1_T0) Enddo ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Convert to TT frame, Duev equation 20. ! T2_T1 = T2_T0 - T1_T0 ! Vacuum delay in TDB frame VEsq = DOTP(VsubEarth,VsubEarth) ! Square of Earth's SSBC velocity UE = GMSUN/vecmg1 ! Newtonian potential, Sun only Call VECSB(SITEP(1,2),SITEP(1,1),B) ! Geocentric baseline vector VEdotB = DOTP(VsubEarth,B) ! Earth SSB velocity * Baseline VEdotV2 = DOTP(VsubEarth,SITEV(1,2)) ! Earth SSB velocity * Site2 geocentric velocity ! td2_td1 = ( (T2_T1/(1.D0-1.48082686741D-8)) * & & ( 1.D0 - (VEsq/2.D0 + UE)/VLIGHT2 ) - VEdotB/VLIGHT2 ) / & & (1.D0 + VEdotV2/VLIGHT2) ! write (6,*) 'Duev: td2_td1: ', td2_td1 ! ! STOP ! RETURN END ! !********************************************************************** ! Subroutine RANGE (DATMC, SITEP, SITEV, SITEA, SUN, T0_T1, & & R1, R1dt, R1mag, R1magdt, R2, R2dt, R2mag, R2magdt, & & tr2_tr1, dtr2_tr1) IMPLICIT NONE ! REAL*8 DATMC(2,2), SUN(3,2), EARTH(3,3), SITEA(3,2), SITEP(3,2), & & SITEV(3,2), DOTP, VECMG Real*8 R1(3), R2(3), R1dt(3), R2dt(3), R1mag, R2mag, R1magdt, & & R2magdt, T0_T1, tr2_tr1, dtr2_tr1 Real*8 SunNFO(3), del2(3), Site2(3),x2x0(3), rsun2(3), vecx2x0(3) Real*8 G2sunV3, rsun0mag, rsun1mag, rsun2mag, x1x0mag, x2x0mag, & & xlog1, xlog2, t1t0, t2t0, t1t2, t1t2_old, G2EarthV3, & & xEx0mag, Site2mag, tgrav, tgravSun, tgravEarth Integer*4 J ! INCLUDE 'd_input.i' ! Variables from: ! 1. SPpxyz(3) - Geocentric satellite position (meters). ! INCLUDE 'cphys11.i' ! Variables 'from': ! 1. VLIGHT - The velocity of light in vacuum. (m/sec) ! 2. VLIGHT2 - The velocity of light squared. (m/sec)**2 ! 3. VLIGHT3 - THE VELOCITY OF LIGHT IN VACUUM CUBED. ! (M**3/SEC**3) ! 4. GMSUN, GMMOON, GMEARTH, GMPLANET(7) - Gravitational constant ! times the masses of the Sun, Moon, Earth, and the ! other planets except Pluto. ! ! From IERS Conventions (2010), section 11.2. ! ! 2*GMsun/VLIGHT3 G2sunV3 = 2.D0*GMSUN/VLIGHT3 Call VECSB(SPpxyz,SUN(1,1),SunNFO) rsun0mag = VECMG(SunNFO) rsun1mag = VECMG(SUN(1,1)) x1x0mag = VECMG(R1) xlog1 = rsun0mag + rsun1mag + x1x0mag xlog2 = rsun0mag + rsun1mag - x1x0mag tgravSun = G2sunV3 * DLOG(xlog1/xlog2) ! G2EarthV3 = 2.D0*GMEARTH/VLIGHT3 xEx0mag = VECMG(SPpxyz) xlog1 = xEx0mag + xEx0mag xlog2 = xEx0mag - xEx0mag ! tgravEarth = G2EarthV3 * DLOG(xlog1/xlog2) tgravEarth = 0.D0 ! zero at geocenter ! ! Delay from NFO to site1 (geocenter) ! t1t0 = R1mag/VLIGHT + tgravEarth + tgravSun t1t0 = R1mag/VLIGHT + tgravEarth ! write(6,*) 'tgravSun,tgravEarth ', tgravSun,tgravEarth ! t1t2 = 0.D0 Do J = 1,4 t1t2_old = t1t2 Call VECMU(SITEV(1,2),t1t2,del2) ! del2 = vector correction to site2 position Call VECSB(SITEP(1,2),del2,Site2) ! Site2 = site2 retarded position Call VECSB(Site2,SPpxyz,vecx2x0) x2x0mag = VECMG(vecx2x0) Call VECSB(Site2,SUN(1,1),rsun2) rsun2mag = VECMG(rsun2) ! xlog1 = rsun0mag + rsun2mag + x2x0mag xlog2 = rsun0mag + rsun2mag - x2x0mag tgravSun = G2sunV3 * DLOG(xlog1/xlog2) ! Site2mag = VECMG(Site2) xlog1 = xEx0mag + Site2mag + x2x0mag xlog2 = xEx0mag + Site2mag - x2x0mag tgravEarth = G2EarthV3 * DLOG(xlog1/xlog2) ! ! t2t0 = x2x0mag/VLIGHT + tgravEarth + tgravSun t2t0 = x2x0mag/VLIGHT + tgravEarth t1t2 = t1t0 - t2t0 ! write(6,*) 'J,t1t0,t2t0,t1t2,t1t2_old,tgravSun,tgravEarth ', J,t1t0,t2t0,t1t2,t1t2_old,tgravSun,tgravEarth Enddo ! write(6,*) 'tgravSun,tgravEarth ', tgravSun,tgravEarth ! tr2_tr1 = t2t0 - t1t0 ! Return End !********************************************************************** SUBROUTINE SEKIDOewns (Datmc_h_EWNS, Datmc_w_EWNS, EARTH, EPBASE, & & SITEP, SITEV, SITEA, SUN, XMOON, T0_T1 ) IMPLICIT NONE ! ! Routine SEKIDO computes the theoretical values of delay and delay ! rate using the Sekido and Fukushima near-field model. ! ! Reference: M. Sekido and T. Fukushima, 'A VLBI delay model for ! radio sources at a finite distance', J. Geodesy, ! vol. 80, pp. 137-149, 2006. ! ! Input Variables: ! 1. DATMC(2,2) - The contributions to the delay and delay rate due ! to tropospheric refraction at each site. (s, s/s) ! 2. EARTH(3,3) - The solar system barycentric Earth position, ! velocity, and acceleration vectors. The first ! index runs over the vector components and the ! second runs over the time derivatives. ! (m, m/sec, m/sec**2) ! 3. EPBASE(3,2) - The J2000.0 baseline position and velocity ! vectors. (m, m/sec) ! 4. SITEP(3,2) - The J2000.0 geocentric position vectors of each ! observing site. (m) ! 5. SITEV(3,2) - The J2000.0 geocentric velocity vectors of each ! observing site. (m/sec) ! 6. SITEA(3,2) - The J2000.0 geocentric acceleration vectors of ! each observing site. (m/sec**2) ! 7. SUN(3,2) - The J2000.0 geocentric Sun position and velocity ! vectors. (m, m/sec) ! 8. XMOON(3,2) - The J2000.0 geocentric Moon position and velocity ! vectors. (m, m/sec) ! 9. STAR(3) - The J2000.0 source unit vector. (unitless) ! ! Output Variables: ! 1. tg2_tg1ewns(4) - The geometric delay corrected for relativistic ! effects (but not atmospheric refraction) using ! the Consensus model. ! 2. dtg2_tg1 - The time derivative of the geometric delay ! corrected for relativistic effects using the ! Consensus model. ! 3. delta_t_grav - The total differential gravitational time delay, ! or "bending delay" from the Consensus model. ! 4. d_delta_t_grav- The time derivative of the total differential ! gravitational time delay, or "bending delay" from ! the Consensus model. ! 5. delta_t_grav_Sun-The Sun's contribution to the "bending delay" ! from the Consensus model. Now includes ! Sunplus(1). ! 6. d_delta_t_grav_Sun-The time derivative of the Sun's contribution ! to the "bending delay" from the Consensus model. ! Now includes Sunplus(2). ! 7. Con_part - Partial derivatives of the Consensus delays and ! rates with respect to Gamma. ! 8. Sunplus(2) - Higher order solar bending delay and rate ! contributions, as defined in IERS Conventions ! (1996), page 91, eqn. 14. (sec, sec/sec) It ! has also been added to the total solar bending ! delay and rate contibutions. ! 9. Bend_par(2) - The relativistic bending delay and rate partials ! with respect to Gamma. (sec, sec/sec) ! ! Common blocks used - ! INCLUDE 'cphys11.i' ! Variables 'from': ! 1. VLIGHT - The velocity of light in vacuum. (m/sec) ! 2. VLIGHT2 - The velocity of light squared. (m/sec)**2 ! 3. VLIGHT3 - The velocity of light cubed. (m/sec)**3 ! 4. GAMMA - The post Newtonian expansion parameter which ! affects light bending. (1.0 used here). (unitless) ! 5. GMSUN, GMMOON, GMEARTH, GMPLANET(7) - Gravitational constant ! times the masses of the Sun, Moon, Earth, and the ! other planets except Pluto. ! 6. REARTH - The equatorial radius of the Earth. (meters) ! INCLUDE 'csolsys11.i' ! Variables 'from': ! 1. SPLANET(3,2,7) - The J2000.0 Solar System Barycentric positions ! and velocities of all planets except the Earth ! and Pluto. (meters, meters/sec) The first index ! runs over X,Y, and Z, the second runs over ! position and velocity, and the third runs over ! the planets, where ! 1 = Mercury ! 2 = Venus ! 3 = Mars ! 4 = Jupiter ! 5 = Saturn ! 6 = Uranus ! 7 = Neptune ! ! 2. GPLANET(3,2,7) - The J2000.0 Geocentric positions and velocities ! of all planets except the Earth and Pluto. ! (meters, meters/sec) The first index runs over ! X,Y, and Z, the second runs over position and ! velocity, and the third runs over the planets, ! where ! 1 = Mercury ! 2 = Venus ! 3 = Mars ! 4 = Jupiter ! 5 = Saturn ! 6 = Uranus ! 7 = Neptune ! Real*8 PI, TWOPI, HALFPI, CONVD, CONVDS, CONVHS, SECDAY COMMON / CMATH / PI, TWOPI, HALFPI, CONVD, CONVDS, CONVHS, SECDAY ! INCLUDE 'ccon.i' ! Variables 'from': ! 1. KTHEC - The Theory routine flow control flag. ! (No longer has any function.) ! 2. KTHED - The theory routine debug control flag. ! 3. KRELC - The relativity module flow control flag. ! 0 --> Gravitational bending used. ! NOT 0 --> Gravitational bending not used. ! INCLUDE 'cobsn11.i' ! Variables from: ! 1. Nzero - Set to 1 or 2 if station 1 or 2 is at the geocenter, ! otherwise equals zero. For correlator usage. ! INCLUDE 'put2s.i' ! Variables from: ! INCLUDE 'd_input.i' ! Variables from: ! 1. SPpxyz(3) - ! 2. SPpvyz(3) - ! Real*8 SpEWNS(3,4), R1EWNS(3,4), R2EWNS(3,4), R1EWNSmag(4), & & R2EWNSmag(4), STAR12EWNS(3,4,2), STARewns(3,4) Common / NFewns/ SpEWNS, R1EWNS, R2EWNS, R1EWNSmag, R2EWNSmag, & & STAR12EWNS, STARewns ! ! ! Program Specifications - REAL*8 SUN(3,2), EARTH(3,3), EPBASE(3,2), SITEA(3,2), SITEP(3,2), & & SITEV(3,2), XMOON(3,2), DOTP, VECMG, Datmc_h_EWNS(4), & & Datmc_w_EWNS(4), T0_T1 Real*8 XsubEarth(3), VsubEarth(3), AsubEarth(3) Real*8 x_sub1t1(3), x_sub2t1(3), w_sub1(3), w_sub2(3) ! Real*8 a_sub1(3), a_sub2(3) Real*8 unit_K(3), b_nought(3),db_nought(3), unitK1(3) Real*8 x_subSun(3), x_subMoon(3), v_subSun(3), v_subMoon(3) Real*8 XsubSun(3), XsubMoon(3), VsubSun(3), VsubMoon(3) Real*8 R_Earth_Sun(3), R_Earth_Moon(3) Real*8 V_Earth_Sun(3), V_Earth_Moon(3) Real*8 Xsub1(3), Xsub2(3), t_sub1, t_sub2 Real*8 dXsub1(3), dXsub2(3) Real*8 x1Sun(3), x1Moon(3), x1Planet(3,7) Real*8 del_t_Sun, del_t_Moon, del_t_Planet(7) Real*8 XSunt1J(3), XMoont1J(3), XPlant1J(3,7) ! Real*8 dXSunt1J(3), dXMoont1J(3), dXPlant1J(3,7) Real*8 R1Sunt1(3), R1Moont1(3), R1Plant1(3,7) Real*8 R2Sunt1(3), R2Moont1(3), R2Plant1(3,7) ! Real*8 dR1Sunt1(3), dR1Moont1(3), dR1Plant1(3,7) ! Real*8 dR2Sunt1(3), dR2Moont1(3), dR2Plant1(3,7) Real*8 delta_t_grav, delta_t_grav_Sun, delta_t_grav_Moon, & & delta_t_grav_Earth, delta_t_grav_Plan(7), & & delta_t_grav_Planets ! Real*8 d_delta_t_grav, d_delta_t_grav_Sun, d_delta_t_grav_Moon, & ! & d_delta_t_grav_Earth, d_delta_t_grav_Plan(7), & ! & d_delta_t_grav_Planets Real*8 delgrav_Sun,ddelgrav_Sun, delgrav_Moon, ddelgrav_Moon, & & delgrav_Plan(7), ddelgrav_Plan(7), delgrav, ddelgrav, & delgrav_Earth, ddelgrav_Earth Real*8 U, U_Sun, dU, absVEarth Real*8 C_Earth, C_Sun, C_Moon, C_Plan(7) Real*8 term_a, term_b, term_c, term_d, term_e, term_f, & & term_g, term_h, term_a1, term_b1, term_g1, term_h1 ! Real*8 dterm_a, dterm_b, dterm_c, dterm_d, dterm_e, dterm_f, & ! & dterm_g, dterm_h, dterm_a1, dterm_b1, dterm_g1, dterm_h1 Real*8 term1, term2a, term2b, term2c, term2d, term2, & & term2bcd, term3a, term3b, term3, term4, & & term123, vec_sum(3), tv2_tv1 ! Real*8 dterm1, dterm2a, dterm2b, dterm2c, dterm2d, dterm2, & ! & dterm2bcd, dterm3a, dterm3b, dterm3, dterm4, & ! & dterm123, dvec_sum(3), dtv2_tv1 Real*8 V_w1(3), V_w2(3), K_V_w1, K_V_w2 !! Real*8 tg2_tg1ewns(4) Real*8 w2_w1(3), a2_a1(3), K_dotw2w1, dK_dotw2w1 ! Real*8 k_sub1(3), k_sub2(3) Real*8 KdotB, dKdotB, VdotB, dVdotB, KdotV, dKdotV Real*8 delta_t2_t1, d_delta_t2_t1, term10a, term10b, term10c, & & dterm10a, dterm10b, dterm10c, del_del_t(2,2) Real*8 x_delta_t2_t1, dx_delta_t2_t1 Real*8 vecmg1, vecmg2, vecmg0, dvecmg0 Real*8 N_hat(3), dN_hat(3), Vmag_S, CSun1, NplusK(3), V1, dV1 !! & V2, dV2 Real*8 R0sun(3), dR0sun(3), R1J, R2J, R10, R20, R0J, & & delTpNsun, R0moon(3), dR0moon(3), delTpNmoon, R0plan(3,7), & & dR0plan(3,7), delTpNPlan(7), delTpN, Vec3c(3), & & term3c, term4a, H, H1(3), H2, H3, dH1a(3), dH1b(3),dH1(3), & & dH2, dH3, dH, V2(3), dV2(3), dVec3c(3), dterm3c, dterm4a Real*8 dR1J, dR2J, dR10, dR20, dR0J, GMp, DelTpNPl, DelTpNPldt, & & delTpNsundt, DelTpNmoondt, DelTpNPlandt(7), delTpNdt Real*8 VsubEart0(3), AsubEart0(3) INTEGER*4 I, L, K, N ! ! Program variables - ! ! 1. XsubEarth(3) = Barycentric radius vector of the geocenter (meters). ! 2. VsubEarth(3) = Barycentric velocity vector of the geocenter (m/sec). ! 3. AsubEarth(3) = Barycentric acceleration of the geocenter (m/sec**2). ! ! 4. x_sub1t1(3) = Geocentric radius vector of receivers 1 and 2 at the ! 5. x_sub2t1(3) geocentric time t_sub1 (meters). ! ! 6. w_sub1(3) = Geocentric velocity of receiver 1 (m/sec). ! 7. w_sub2(3) = Geocentric velocity of receiver 2 (m/sec). ! ! 8. a_sub1(3) = Geocentric acceleration of receiver 1 (m/sec). ! 9. a_sub2(3) = Geocentric acceleration of receiver 2 (m/sec). ! ! 10. x_subSun(3) = Geocentric radius vector of the Sun (meters). ! 11. x_subMoon(3) = Geocentric radius vector of the Moon (meters). ! 12. GPLANET(l,1,k)= Geocentric radius vectors of the planets (meters). ! (l = X,Y,Z; k = planet #) ! ! 13. v_subSun(3) = Geocentric velocity vector of the Sun (m/sec). ! 14. v_subMoon(3) = Geocentric velocity vector of the Moon (m/sec). ! 15. GPLANET(l,2,k)= Geocentric velocity vector of the planets (m/sec). ! (l = X,Y,Z; k = planet #) ! ! 16. Xsub1(3) = Barycentric radius vector of receiver 1 (meters). ! 17. Xsub2(3) = Barycentric radius vector of receiver 2 (meters). ! 18. dXsub1(3) = Barycentric velocity vector of receiver 1 (m/sec). ! 19. dXsub2(3) = Barycentric velocity vector of receiver 2 (m/sec). ! ! 20. x1Sun(3) = Vector from receiver 1 to the Sun at time t_sub1 ! 21. x1Moon(3) = Vector from receiver 1 to the Moon at time t_sub1 ! 22. x1Planet(3,7) = Vector from receiver 1 to the planets at time t_sub1 ! ! 23. XsubSun(3) = Barycentric radius vector of the Sun/Moon/planets ! 24. XsubMoon(3) ! 25. SPLANET(l,1,k) ! ! 26. VsubSun(3) = Barycentric velocity vector of the Sun/Moon/planets ! 27. VsubMoon(3) ! 28. SPLANET(l,2,k) ! ! 29. XSunt1J(3) = SSBC radius vector to the Sun/Moon/planets ! 30. XMoont1J(3) at the time of closest approach. ! 31. XPlant1J(3,7) (meters) ! ! 32. R_Earth_Sun(3) = Vector from the Sun to the geocenter ! 33. R_Earth_Moon(3)= Vector from the Moon to the geocenter ! 34. -GPLANET(l,1,k)= Vector from planets to the geocenter ! ! 35. V_Earth_Sun(3) = Velocity Vector of the geocenter from the Sun ! ! 36. R1Sunt1(3) = Vector from the Sun/Moon/planets to receiver ! 37. R1Moont1(3) = 1 at the time of closest approach. ! 38. R1Plant1(3,7) = (meters) ! ! 39. R2Sunt1(3) = Vector from the Sun/Moon/planets to receiver ! 40. R2Moont1(3) = 2 at the time of closest approach. ! 41. R2Plant1(3,7) = (meters) ! ! 42. t_sub1 = Time of arrival of the signal at receivers 1 and 2. ! 43. t_sub2 [Note: not needed here so we don't actually define ! them.] ! ! 44. del_t_Sun = Difference between arrival time at receiver 1, t_sub1, ! 45. del_t_Moon and time of closest approach to Sun/Moon/planets ! 46. del_t_Planet(7) ( > 0 if closest approach is before arrival) ! ! 47. unit_K(3) = Barycentric unit vector to the source (in the absense ! of gravitational or aberrational bending). ! ! 48. k_sub1(3) = Aberrated unit vector from station 1 to the source. ! 49. k_sub2(3) = Aberrated unit vector from station 2 to the source. ! [Note: not computed here, see atmosphere module.] ! ! 50. b_nought(3) = A priori geocentric baseline vector at time t_sub1 ! (Baseline vector and it's velocity. Defined from ! site #1 to site #2 (M,M/S).) ! ! 51. U = Sun's potential ! ! 52. delgrav_Sun ! 53. delgrav_Moon = The differential gravitational time delay for ! 54. delgrav_Plan(7) the Sun/Moon/planets/Earth. ! 55. delgrav_Earth ! ! 56. delta_t_grav_Planets = Sum of bending delays for the planets ! ! 57. delta_t_grav = The total differential gravitational time delay, ! or "bending delay." ! ! 58. tv2_tv1 = The theoretical vacuum delay computed using the ! Consensus model. (sec) ! ! 59. tg2_tg1 = The total theoretical geometric delay (doesn't ! include tropospheric refraction effects which are ! usually added in Solve) computed using the Consensus ! model. (sec) ! ! 4.2.9 PROGRAMMER - DAVID GORDON 04/22/93 thru 08/02/93 - Written and debuged ! D. Gordon Nov. 1993 - Modified to use all planets ! except Pluto. ! D. Gordon Jan-Mar 1994 - Modified for axis offset ! correction. Not using Eubank's Step 10. ! D. Gordon 10.05.94 Many unused variable and much unused ! code removed. ! D. Gordon 96.01.30 Added section to compute partial ! derivatives of delay and rate w.r.t. Gamma. ! D. Gordon 96.02.09 Corrected typo in Step 2 debug printout. ! D. Gordon 98.08.18 Mods for geocenter station. ! D. Gordon 98.11.16 Added computation of higher order ! solar bending term, from IERS Conventions ! (1996), page 91, eqn. 14. This term is now ! added to the total solar bending term. It ! only becomes significant within about 2 ! degrees of the Sun. ! D. Gordon 98.12.16 Added computation of Bend_par(2), the ! partial derivatives of the gravitational ! bending contributions with respect to Gamma. ! ! Write(6,*) ' ' ! Write(6,*) ' SEKIDO: ' ! !_____________________________________________________________________________ ! ! Copy CALC variables into variables with names which mimic the variables ! in the consensus paper. VLIGHT is used for the velocity of light because ! 'C' is a poor variable name that could be easily lost. ! Do l=1,3 ! XsubEarth(l) = EARTH(l,1) ! SSBC Earth position VsubEarth(l) = EARTH(l,2) ! ditto velocity AsubEarth(l) = EARTH(l,3) ! ditto acceleration ! VsubEart0(l) = 0.0D0 ! modified velocity AsubEart0(l) = 0.0D0 ! modified acceleration ! x_sub1t1(l) = SITEP(l,1) ! Site 1 geocentric position w_sub1(l) = SITEV(l,1) ! ditto velocity !d a_sub1(l) = SITEA(l,1) ! ditto acceleration ! x_sub2t1(l) = SITEP(l,2) ! Site 2 geocentric position w_sub2(l) = SITEV(l,2) ! ditto velocity !d a_sub2(l) = SITEA(l,2) ! ditto acceleration ! ! unit_K(l) = STAR(l) ! J2000.0 pseudo source vector ! unitK1(l) = STAR12(l,1) ! x_subSun(l) = SUN(l,1) ! Geocentric Sun position v_subSun(l) = SUN(l,2) ! ditto velocity ! x_subMoon(l) = XMOON(l,1) ! Geocentric Moon position v_subMoon(l) = XMOON(l,2) ! ditto velocity ! b_nought(l) = -EPBASE(l,1) ! Baseline vector from site 1 to site 2 db_nought(l) = -EPBASE(l,2) ! Time derivative of baseline vector ! R_Earth_Sun(l) = -SUN(l,1) V_Earth_Sun(l) = -SUN(l,2) ! = -v_subSun(l) also R_Earth_Moon(l) = -XMOON(l,1) ! ! SSBC radius vectors of Sun and Moon: XsubSun(l) = XsubEarth(l) + x_subSun(l) XsubMoon(l) = XsubEarth(l) + x_subMoon(l) ! ! SSBC velocity vectors of Sun and Moon: VsubSun(l) = VsubEarth(l) + v_subSun(l) VsubMoon(l) = VsubEarth(l) + v_subMoon(l) ! Enddo ! !****************************************************************************** ! Step 1: Estimate barycentric radius and velocity vectors for stations 1 ! and 2 at time t_sub1 ! Do l=1,3 Xsub1(l) = XsubEarth(l) + x_sub1t1(l) Xsub2(l) = XsubEarth(l) + x_sub2t1(l) !d dXsub1(l) = VsubEarth(l) + w_sub1(l) !derivative !d dXsub2(l) = VsubEarth(l) + w_sub2(l) !derivative Enddo ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! Start East, West, North, South offsets loop Do N = 1, 4 ! do l = 1,3 unitK1(l) = STAR12EWNS(l,N,1) enddo ! !****************************************************************************** ! Step 2: Estimate the vectors from the Sun, the Moon, and each planet (except ! Earth and Pluto) to receiver 1. ! ! Find time of closet approach of the near-field object to the ! gravitating body. [Actually ........................ ! we find only how much earlier (or later) the ........ rays ! passed closest to the gravitating body. Then, if earlier, we ! extrapolate the body's position back to that earlier time. If ! not earlier we just keep the current position.] ! ! J2000.0 vector from receiver #1 to the Sun/Moon/planets: Do l=1,3 x1Sun(l) = XsubSun(l) - Xsub1(l) x1Moon(l) = XsubMoon(l) - Xsub1(l) Enddo Do k=1,7 do l=1,3 x1Planet(l,k) = SPLANET(l,1,k) - Xsub1(l) enddo Enddo ! ! more good stuff KdotB = DOTP(STARewns(1,N),b_nought) !d dKdotB = DOTP(STAR,db_nought) + DOTP(STARdt,b_nought) ! ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! Sun: del_t_Sun = DOTP(unitK1,x1Sun) / Vlight ! ! write(6,8)' T0_T1, del_t_Sun(sec) ', T0_T1, del_t_Sun If(del_t_Sun .lt. 0.D0) del_t_Sun = 0.0D0 If(del_t_Sun .gt. T0_T1) del_t_Sun = T0_T1 ! write(6,8)' T0_T1, del_t_Sun(sec) ', T0_T1, del_t_Sun ! Do l=1,3 ! SSBC vector to Sun at time of closest approach and its time derivative: XSunt1J(l) = Xsubsun(l) - VsubSun(l)*del_t_Sun !d dXSunt1J(l) = Vsubsun(l) !Derivative (approx. - no acceleration) ! ! Vector from Sun to receiver 1 R1Sunt1(l) = Xsub1(l) - XSunt1J(l) !d dR1Sunt1(l) = dXsub1(l) - dXSunt1J(l) !Derivative ! ! Vector from Sun to receiver 2 R2Sunt1(l) = Xsub2(l) - VsubEarth(l)*KdotB/VLIGHT - XSunt1J(l) !d dR2Sunt1(l) = dXsub2(l) - AsubEarth(l)*KdotB/VLIGHT - & !d & VsubEarth(l)*dKdotB/VLIGHT - dXSunt1J(l) !Derivative ! ! Vector from Sun to near-field object !!! R0sun(l) = R1Sunt1(l) + R1(l) R0sun(l) = R1Sunt1(l) + R1EWNS(l,N) !d dR0sun(l) = dR1Sunt1(l) + R1dt(l) ! Enddo ! write(6,8)' Xsubsun,VsubSun ', Xsubsun,VsubSun ! write(6,8)' Xsub1, XSunt1J ', Xsub1, XSunt1J ! write(6,8)' Xsub2, VsubEarth ', Xsub2, VsubEarth ! write(6,8)' KdotB,VLIGHT ', KdotB,VLIGHT ! write(6,8)' XSunt1J,dXSunt1J ', XSunt1J,dXSunt1J ! write(6,8)' R1Sunt1,dR1Sunt1 ', R1Sunt1,dR1Sunt1 ! write(6,8)' R2Sunt1,dR2Sunt1 ', R2Sunt1,dR2Sunt1 ! write(6,8)' R0sun,R1 ', R0sun,R1 ! R1J = VECMG (R1Sunt1) !d dR1J = DOTP(R1Sunt1,dR1Sunt1)/R1J R2J = VECMG (R2Sunt1) !d dR2J = DOTP(R2Sunt1,dR2Sunt1)/R2J R0J = VECMG (R0sun) !d dR0J = DOTP(R0Sun,dR0Sun)/R0J !!!!! R10 = R1mag R10 = R1EWNSmag(N) !d dR10 = R1magdt !!!!! R20 = R2mag R20 = R2EWNSmag(N) !d dR20 = R2magdt ! Write(6,*) ' Sun: ' ! Write(6,*) ' R1,R1Sunt1 ', R1, R1Sunt1 ! Write(6,*) ' R0sun ', R0sun ! Write(6,*) ' R2,R2Sunt1 ', R2, R2Sunt1 ! Write(6,*) ' R1J, dR1J ', R1J, dR1J ! Write(6,*) ' R2J, dR2J ', R2J, dR2J ! Write(6,*) ' R0J, dR0J ', R0J, dR0J ! Write(6,*) ' R10, dR10 ', R10, dR10 ! Write(6,*) ' R20, dR20 ', R20, dR20 ! Write(6,*) ' ' ! CALL TGRAV1(GMSUN, R0J,R1J,R2J,R10,R20, delgrav_Sun) ! ! ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! Moon: del_t_Moon = DOTP(unitK1,x1Moon) / Vlight ! If(del_t_Moon .lt. 0.D0) del_t_Moon = 0.0D0 If(del_t_Moon .gt. T0_T1) del_t_Moon = T0_T1 ! write(6,8)' del_t_Moon ',del_t_Moon ! Do l=1,3 ! SSBC vector to Moon at time of closest approach and its time derivative: XMoont1J(l) = XsubMoon(l) - VsubMoon(l)*del_t_Moon !d dXMoont1J(l) = VsubMoon(l) !Derivative (approx. - no acceleration) ! ! Vector from the Moon to receiver 1 R1Moont1(l) = Xsub1(l) - XMoont1J(l) !d dR1Moont1(l) = dXsub1(l) - dXMoont1J(l) !Derivative ! ! Vector from the Moon to receiver 2 R2Moont1(l) = Xsub2(l) - VsubEarth(l)* KdotB/VLIGHT - & & XMoont1J(l) !d dR2Moont1(l) = dXsub2(l) - AsubEarth(l)*KdotB/VLIGHT - & !d & VsubEarth(l)*dKdotB/VLIGHT - dXMoont1J(l) !Derivative ! ! Vector from Moon to near-field R0moon(l) = R1moont1(l) + R1EWNS(l,N) !d dR0moon(l) = dR1moont1(l) + R1dt(l) ! Enddo ! R1J = VECMG (R1Moont1) !d dR1J = DOTP(R1Moont1,dR1Moont1)/R1J R2J = VECMG (R2Moont1) !d dR2J = DOTP(R2Moont1,dR2Moont1)/R2J R0J = VECMG (R0Moon) !d dR0J = DOTP(R0Moon,dR0Moon)/R0J ! ! CALL TGRAV1(GMMOON,R0J,R1J,R2J,R10,R20, delgrav_Moon) ! write(6,8)'TGRAV: delgrav_Moon ', delgrav_Moon, ddelgrav_Moon ! ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! Planets: Do k=1,7 ! Planet loop del_t_Planet(k) = DOTP(unitK1,x1Planet(1,k)) / Vlight ! If(del_t_Planet(k) .lt. 0.D0) del_t_Planet(k) = 0.0D0 If(del_t_Planet(k) .gt. T0_T1) del_t_Planet(k) = T0_T1 ! write(6,'("del_t_Planet(",i1,") = ",d25.16)')k,del_t_Planet(k) ! Do l=1,3 ! SSBC vector to Planet at time of closest approach and its time derivative: XPlant1J(l,k) = SPLANET(l,1,k) - SPLANET(l,2,k)*del_t_Planet(k) !d dXPlant1J(l,k) = SPLANET(l,2,k) ! Derivative (approximate) ! ! Vector from Planet to receiver 1 R1Plant1(l,k) = Xsub1(l) - XPlant1J(l,k) !d dR1Plant1(l,k) = dXsub1(l) - dXPlant1J(l,k) !Derivative ! ! Vector from Planet to receiver 2 R2Plant1(l,k) = Xsub2(l) - VsubEarth(l)*KdotB/VLIGHT - & & XPlant1J(l,k) !d dR2Plant1(l,k) = dXsub2(l) - AsubEarth(l)*KdotB/VLIGHT - & !d & VsubEarth(l)*dKdotB/VLIGHT - dXPlant1J(l,k) !Derivative ! ! Vector from Planet k to near-field object R0plan(l,k) = R1plant1(l,k) + R1EWNS(l,N) !d dR0plan(l,k) = dR1plant1(l,k) + R1dt(l) ! Enddo ! R1J = VECMG (R1Plant1(1,k)) !d dR1J = DOTP(R1Plant1(1,k),dR1Plant1(1,k))/R1J R2J = VECMG (R2Plant1(1,k)) !d dR2J = DOTP(R2Plant1(1,k),dR2Plant1(1,k))/R2J R0J = VECMG (R0plan(1,k)) !d dR0J = DOTP(R0plan(1,k),dR0plan(1,k))/R0J GMp = GMPlanet(k) ! ! Write(6,*) ' Planet ', k ! Write(6,*) ' R1J, dR1J ', R1J, dR1J ! Write(6,*) ' R2J, dR2J ', R2J, dR2J ! Write(6,*) ' R0J, dR0J ', R0J, dR0J ! Write(6,*) ' GMp ', GMp ! Write(6,*) ' ' ! CALL TGRAV1 (GMp,R0J,R1J,R2J,R10,R20, DelTpNPl) ! write(6,*)' #, delTpNPLAN(k) ', k, delTpNPlan(k) delgrav_Plan(k) = DelTpNPl !d ddelgrav_Plan(k) = DelTpNPldt ! write(6,*)'TGRAV: delgrav_Plan ', k, delgrav_Plan(k), ddelgrav_Plan(k) ! Enddo ! Planet loop ! !****************************************************************************** ! Step 4: Find the differential delay due to the Earth. ! C_Earth = (1.0D0 + gamma) * GMEarth/VLIGHT3 vecmg0 = VECMG(SpEWNS(1,N)) !d dvecmg0 = DOTP(SPvxyz,SPpxyz)/vecmg0 ! IF(Nzero .ne. 1) THEN vecmg1 = VECMG(x_sub1t1) term_g = vecmg1 + DOTP(STAR12EWNS(1,N,1),x_sub1t1) !d dterm_g = DOTP(x_sub1t1,w_sub1)/vecmg1 + & !d & DOTP(STAR12(1,1),w_sub1) + DOTP(STAR12dt(1,1),x_sub1t1) term_g1 = vecmg0 + DOTP(STAR12EWNS(1,N,1),SpEWNS(1,N)) !d dterm_g1 = dvecmg0 + DOTP(STAR12(1,1),SPvxyz) & !d & + DOTP(STAR12dt(1,1),SPpxyz) ELSE vecmg1 = REARTH term_g = 2.0D0 * REARTH !d dterm_g = 0.0D0 term_g1 = 2.0D0 * vecmg0 !d dterm_g1 = 2.0D0 * dvecmg0 ENDIF ! IF(Nzero .ne. 2) THEN vecmg2 = VECMG(x_sub2t1) term_h = vecmg2 + DOTP(STAR12EWNS(1,N,2),x_sub2t1) !d dterm_h = DOTP(x_sub2t1,w_sub2)/vecmg2 + & !d & DOTP(STAR12(1,2),w_sub2) + DOTP(STAR12dt(1,2),x_sub2t1) term_h1 = vecmg0 + DOTP(STAR12EWNS(1,N,2),SpEWNS(1,N)) !d dterm_h1 = dvecmg0 + DOTP(STAR12(1,2),SPvxyz) & !d & + DOTP(STAR12dt(1,2),SPpxyz) ELSE vecmg2 = REARTH term_h = 2.0D0 * REARTH !d dterm_h = 0.0D0 term_h1 = 2.0D0 * vecmg0 !d dterm_h1 = 2.0D0 * dvecmg0 ENDIF ! delgrav_Earth = C_Earth * DLOG( (term_g*term_h1) / & & (term_h*term_g1) ) !d ddelgrav_Earth = C_Earth * ( dterm_g/term_g + & !d & dterm_h1/term_h1 - dterm_h/term_h - dterm_g1/term_g1 ) ! !****************************************************************************** ! Step 5: Add up all components to get the total ! differential gravitational delay. ! [Does not include the term for observations close to the Sun.] ! ! Add up the 10 gravitational components delgrav = delgrav_Sun + delgrav_Moon + delgrav_Earth !d ddelgrav = ddelgrav_Sun + ddelgrav_Moon + ddelgrav_Earth Do k=1,7 delgrav = delgrav + delgrav_Plan(k) !d ddelgrav = ddelgrav + ddelgrav_Plan(k) Enddo ! write(6,8)' Sekido: delgrav,ddelgrav ', delgrav,ddelgrav ! !****************************************************************************** ! Step 6: Add the total differential gravitational delay to the rest of the ! a priori vacuum delay, Equation 9. ! ! Find Sun's potential: ! U_Sun = GMSUN/VLIGHT2 vecmg1 = VECMG(R_Earth_Sun) U = U_Sun/vecmg1 ! Derivative: !d dU = -U_Sun * Dotp(R_Earth_Sun,V_Earth_Sun) / vecmg1**3 ! ! Compute individual terms of Eqn. 20 ! term2a = KdotB/VLIGHT !d dterm2a = dKdotB/VLIGHT ! derivative ! Write(6,8) ' term2a,dterm2a ', term2a,dterm2a ! term2b = 1.D0 - ((1.D0 + gamma) * U) !d dterm2b = -(1.D0 + gamma) * dU ! derivative ! Write(6,8) ' term2b-1,dterm2b ', (term2b-1.D0),dterm2b ! absVEarth = VECMG(VsubEarth) term2c = (absVEarth)**2 / (2.D0*VLIGHT2) !d dterm2c = Dotp(VsubEarth,AsubEarth) / VLIGHT2 ! derivative !! term2c = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! dterm2c = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Write(6,8) ' term2c,dterm2c ', term2c,dterm2c ! term2d = DOTP(VsubEarth,w_sub2) / VLIGHT2 !d dterm2d = (DOTP(AsubEarth,w_sub2) + DOTP(VsubEarth,a_sub2)) & !d & / VLIGHT2 ! derivative !! term2d = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! dterm2d = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Write(6,8) ' term2d,dterm2d ', term2d,dterm2d ! ! Combine terms 2b,2c,2d term2bcd = term2b - term2c - term2d !d dterm2bcd = dterm2b - dterm2c - dterm2d ! derivative ! term2 = term2a * term2bcd !d dterm2 = term2a * dterm2bcd + dterm2a * term2bcd ! derivative ! Write(6,8) ' term2, dterm2 ', term2,dterm2 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! VdotB = DOTP(VsubEart0,b_nought) !d dVdotB = DOTP(AsubEarth,b_nought) + DOTP(VsubEarth,db_nought) !derivative ! term3a = VdotB/VLIGHT2 !d dterm3a = dVdotB/VLIGHT2 ! Write(6,8) ' term3a,dterm3a ', term3a,dterm3a ! Call VECAD (VsubEart0, w_sub2, V2 ) !d Call VECAD (AsubEarth, a_sub2, dV2 ) ! Write(6,8) ' V2, dV2 ', V2, dV2 term3b = DOTP(STAR12EWNS(1,N,2), V2)/VLIGHT !d dterm3b = (DOTP(STAR12dt(1,2),V2) + DOTP(STAR12(1,2),dV2))/VLIGHT ! Write(6,8) ' STAR12(1,2) ', STAR12(1,2),STAR12(2,2),STAR12(3,2) ! Write(6,8) ' STAR12dt(1,2) ', STAR12dt(1,2),STAR12dt(2,2),STAR12dt(3,2) ! Write(6,8) ' term3b,dterm3b ', term3b,dterm3b ! DO I = 1,3 Vec3c(I) = VsubEart0(I) + 2.D0*w_sub2(I) !d dVec3c(I) = AsubEarth(I) + 2.D0*a_sub2(I) ENDDO term3c = DOTP(Vec3c,STARewns(1,N)) / (2.D0*VLIGHT) !d dterm3c = (DOTP(dVec3c,STAR) + DOTP(Vec3c,STARdt))/(2.D0*VLIGHT) ! Write(6,8) ' term3c,dterm3c ', term3c,dterm3c ! term3 = term3a * (1.D0 + term3b - term3c) !d dterm3 = dterm3a*(1.D0+term3b-term3c) + term3a*(dterm3b-dterm3c) !! term3 = term3a * (1.D0 ) !!!!!!!!!!!!!!!!! !! dterm3 = dterm3a*(1.D0) !!!!!!!!!!!!!!!!! !? term3 = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!! !? dterm3 = 0.D0 !!!!!!!!!!!!!!!!!!!!!!!!! ! Write(6,8) ' term3, dterm3 ', term3,dterm3 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! term4a = 1.D0 + DOTP(STAR12EWNS(1,N,2),V2)/VLIGHT term4a = 1.D0 + term3b !d dterm4a = (DOTP(STAR12dt(1,2),V2) + DOTP(STAR12(1,2),dV2))/VLIGHT ! Write(6,8) ' term4a,dterm4a ', term4a,dterm4a ! Call CROSP(V2, STAR12EWNS(1,N,2), H1) H2 = (VECMG(H1)/VLIGHT)**2 H3 = KdotB/(2.D0*R2EWNSmag(N)) H = H2*H3 !d Call CROSP(dV2,STAR12(1,2),dH1a) !d Call CROSP(V2,STAR12dt(1,2),dH1b) !d Call VECAD (dH1a, dH1b, dH1) ! dH2 = 2.D0*(VECMG(H1)*(VECMG(dH1)))/(VLIGHT2) !d dH2 = (2.D0 * DOTP(H1,dH1)) / VLIGHT2 !d dH3 = (dKdotB/R2mag - KdotB*R2magdt/R2mag**2)/2.0D0 !d dH = dH2*H3 + H2*dH3 ! Write(6,8) ' H1, dH1 ', H1, dH1 ! Write(6,8) ' H2, dH2 ', H2, dH2 ! Write(6,8) ' H3, dH3 ', H3, dH3 ! Write(6,8) ' H, dH ', H, dH ! term4 = term4a * (1.D0 + H) !d dterm4 = dterm4a*(1.D0 + H) + term4a*dH !!? term4 = (1.D0 + H) !!!!!!!!!!!!!!!!!!!!!! !!? dterm4 = 0.D0 !!!!!!!!!!!!!!!!!!!!!! ! Write(6,8) ' term4, dterm4 ', term4,dterm4 ! tv2_tv1 = (-term2 - term3 + delgrav) / term4 !d dtv2_tv1 = (-dterm2 - dterm3 + ddelgrav) / term4 - & !d & (1.D0/term4**2) * (-term2 - term3 + delgrav) * dterm4 ! ! write(6,8) 'S: tv2_tv1, dtv2_tv1 ',tv2_tv1, dtv2_tv1 ! !****************************************************************************** ! Step 8: Add the geometric part of the tropospheric propogation delay to the ! vacuum delay. ! [Geocenter station: DATMC(Nzero,1) and DATMC(Nzero,2) should ! already be zero.] ! call vecsb(w_sub2,w_sub1,w2_w1) !d call vecsb(a_sub2,a_sub1,a2_a1) !??? K_dotw2w1 = DOTP(STARewns(1,N),w2_w1) + DOTP(STARdt,b_nought) K_dotw2w1 = DOTP(STARewns(1,N),w2_w1) !d dK_dotw2w1 = DOTP(STARdt,w2_w1) + DOTP(STAR,a2_a1) + & !d & DOTP(STARdt,db_nought) !!! & + DOTP(STARdtdt,b_nought) ! ! Negative sign because station 1 atmosphere is negative. !?? tg2_tg1ewns(N) = tv2_tv1 - Datmc_h_EWNS(N) * K_dotw2w1/VLIGHT tg2_tg1ewns(N) = tv2_tv1 !d dtg2_tg1 = dtv2_tv1 - DATMC(1,2) * K_dotw2w1/VLIGHT & !d & - DATMC(1,1) * dK_dotw2w1/VLIGHT ! write(6,8)' w2_w1, a2_a1 ', w2_w1, a2_a1 ! If (UVW .eq. 'exact ') Then tg2_tg1ewns(N) = tg2_tg1ewns(N) + Datmc_h_EWNS(N) + Datmc_w_EWNS(N) ! write(6,'("I,tg2_tg1(psec):",I2,F18.2)') I, tg2_tg1(I)*1.D12 Endif ! Enddo ! East, West, North, South offsets loop !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! 8 FORMAT(A,4D25.16/(7X,5D25.16)) ! write (6,*) ' Sekido model ' ! write (6,*) ' tg2_tg1east ', tg2_tg1ewns(1) ! write (6,*) ' tg2_tg1west ', tg2_tg1ewns(2) ! write (6,*) ' tg2_tg1north ', tg2_tg1ewns(3) ! write (6,*) ' tg2_tg1south ', tg2_tg1ewns(4) ! U_V(1) = VLIGHT*(tg2_tg1ewns(1) - tg2_tg1ewns(2))/0.0002D0 U_V(2) = VLIGHT*(tg2_tg1ewns(3) - tg2_tg1ewns(4))/0.0002D0 ! write(6,*) ' U (m, sec) ', U_V(1), U_V(1)/VLIGHT ! write(6,*) ' V (m, sec) ', U_V(2), U_V(2)/VLIGHT ! RETURN END ! !********************************************************************** SUBROUTINE TGRAV1(GM, R0J,R1J,R2J,R10,R20, DelTpN) IMPLICIT NONE ! ! Subroutine to compute the gravitational deflection for a solar ! system body (Sun, Moon, planets) on a near-field object. ! Same as TGRAV but without time derivatives. ! ! ! 1. GM - Gravitational constant times the mass of the gravitating ! body. Real*8 GM, R0J, R1J, R2J, R10, R20, DelTpN Real*8 term1a, term1b, term2a, term2b, term1adt, term1bdt, & & term2adt, term2bdt, term12, term12dt !* INCLUDE 'cphys11.i' ! Variables 'from': ! 1. VLIGHT - The velocity of light in vacuum. (m/sec) ! 2. VLIGHT2 - The velocity of light squared. (m/sec)**2 ! 3. VLIGHT3 - The velocity of light cubed. (m/sec)**3 ! ! term1a = R2J + R0J + R20 ! term1b = R2J + R0J - R20 ! term2a = R1J + R0J - R10 ! term2b = R1J + R0J + R10 ! term12 = (term1a/term1b) * (term2a/term2b) ! DelTpN = (2.D0*GM/VLIGHT3) * DLOG (term12) ! RETURN END ! !********************************************************************** Subroutine DUEVewns (Datmc_h_EWNS, Datmc_w_EWNS, SITEP, SITEV, & & SITEA, SUN, T0_T1, EARTH, XMOON, UTC, TT, TDB, TDBg) ! IMPLICIT NONE REAL*8 SUN(3,2), EARTH(3,3), SITEA(3,2), SITEP(3,2), & & SITEV(3,2), DOTP, VECMG, UTC, TT, TDB, TDBg, XMOON(3,2), & & Datmc_h_EWNS(4), Datmc_w_EWNS(4) Real*8 td2_td1(4), dtd2_td1(4), T0_T1 Real*8 U, dU, U_Sun, Vecmg1, F1, VdotR1, VdotR2, VdotNFO, & & del_T0, gm1, RLT_01, R_01mag, R0sunT0mag, R1sunT1mag, & & T1_T0, ddel_T0, R01sunmag Real*8 VsubEarth(3), R_Earth_Sun(3), V_Earth_Sun(3), R1_T1(3), & & R2_T1(3), R0_T0(3), dR0_T0(3), R_01(3), delSun0(3), & & SUN_0(3), R0sunT0(3), R1sunT1(3), Rb_01(3), dR0(3), & & R1(3), EARTH_0(3), RsunT0(3), RsunT1(3), R01sun(3) Real*8 R_02(3), Site2(3), R2_T2(3), SunT2(3), EarthT2(3), & & RsunT2(3), R2sunT2(3), R02sun(3) Real*8 T2_T0, del_T2, R2_T2mag, R_02mag, R2sunT2mag, R02sunmag, & & RLT_02, ddel_T2 Real*8 T2_T1, VEsq, UE, B(3), VEdotB, VEdotV2 Real*8 delMoon0(3),MOON_0(3),RmoonT0(3),R0moonT0(3),R0moonT0mag, & & RmoonT1(3),R1moonT1(3),R1moonT1mag,R01moon(3),R01moonmag, & & RLT_01s,RLT_01m,RLT_01p,delpl(3),Spl_0(3),R0planT0(3), & & R0planT0mag,R1planT1(3),R1planT1mag,R01plan(3),R01planmag, & & RLTp(7) Real*8 MoonT2(3),RmoonT2(3),R2moonT2(3),R2moonT2mag,R02moon(3), & & R02moonmag,RLT_02s,RLT_02m,RLT_02p,Spl_2(3),R2planT2(3), & & R2planT2mag,R02plan(3),R02planmag Real*8 delEarth0(3), REarthT0(3), R0earthT0(3), R0earthT0mag, & & RearthT1(3), R1earthT1(3), R1earthT1mag, R01earth(3), & & R01earthmag, RLT_01e, RearthT2(3), R2earthT2(3), & & R2earthT2mag, R02earth(3), R02earthmag, RLT_02e, & & delearth2(3), delSun2(3), delmoon2(3) Integer*4 I, J, K, N ! INCLUDE 'd_input.i' ! Variables from: ! 1. SPpxyz(3) - Geocentric satellite position at T0 (meters). ! Real*8 SpEWNS(3,4), R1EWNS(3,4), R2EWNS(3,4), R1EWNSmag(4), & & R2EWNSmag(4), STAR12EWNS(3,4,2), STARewns(3,4) Common / NFewns/ SpEWNS, R1EWNS, R2EWNS, R1EWNSmag, R2EWNSmag, & & STAR12EWNS, STARewns ! Variables from: ! 1. SpEWNS(3,4) - Geocentric satellite East, West, North, and ! South offset positions at T0 (meters). ! INCLUDE 'put2s.i' ! Variables to: ! INCLUDE 'cphys11.i' ! Variables 'from': ! 1. VLIGHT - The velocity of light in vacuum. (m/sec) ! 2. VLIGHT2 - The velocity of light squared. (m/sec)**2 ! 3. VLIGHT3 - THE VELOCITY OF LIGHT IN VACUUM CUBED. ! (M**3/SEC**3) ! 4. GMSUN, GMMOON, GMEARTH, GMPLANET(7) - Gravitational constant ! times the masses of the Sun, Moon, Earth, and the ! other planets except Pluto. ! INCLUDE 'csolsys11.i' ! Variables 'from': ! 1. SPLANET(3,2,7) - The J2000.0 Solar System Barycentric positions ! and velocities of all planets except the Earth ! and Pluto. (meters, meters/sec) The first index ! runs over X,Y, and Z, the second runs over ! position and velocity, and the third runs over ! the planets, where ! 1 = Mercury ! 2 = Venus ! 3 = Mars ! 4 = Jupiter ! 5 = Saturn ! 6 = Uranus ! 7 = Neptune ! ! 2. GPLANET(3,2,7) - The J2000.0 Geocentric positions and velocities ! of all planets except the Earth and Pluto. ! (meters, meters/sec) The first index runs over ! X,Y, and Z, the second runs over position and ! velocity, and the third runs over the planets, ! where ! 1 = Mercury ! 2 = Venus ! 3 = Mars ! 4 = Jupiter ! 5 = Saturn ! 6 = Uranus ! 7 = Neptune ! ! 3. SUNb(3,2) - J2000 Barycentric Sun position and velocity. ! 4. MOONb(3,2) - J2000 Barycentric Moon position and velocity. ! ! Convert spacecraft and site vectors from TT to TDB Do I = 1, 3 VsubEarth(I) = EARTH(I,2) ! SSBC Earth Velocity R_Earth_Sun(I) = -SUN(I,1) ! Vector from Sun to Earth V_Earth_Sun(I) = -SUN(I,2) ! Velocity of Earth w.r.t. the Sun Enddo ! ! Sun's potential: U_Sun = GMSUN/VLIGHT2 vecmg1 = VECMG(R_Earth_Sun) U = U_Sun/vecmg1 ! Derivative: ! dU = -U_Sun * Dotp(R_Earth_Sun,V_Earth_Sun) / vecmg1**3 F1 = 1.D0 - U - 1.48082686741D-8 ! VdotR1 = DOTP(VsubEarth,SITEP(1,1)) ! VdotR2 = DOTP(VsubEarth,SITEP(1,2)) ! ! Barycentric TDB vectors: (Duev equation 4) Do I = 1,3 R1_T1(I) = SITEP(I,1)*F1 - 0.5D0*(VdotR1/VLIGHT2)*VsubEarth(I) & & + EARTH(I,1) R2_T1(I) = SITEP(I,2)*F1 - 0.5D0*(VdotR2/VLIGHT2)*VsubEarth(I) & & + EARTH(I,1) Enddo ! write(6,*) 'site1,R1_T1: ', SITEP(1,1),SITEP(2,1),SITEP(3,1),R1_T1 ! write(6,*) 'site2,R2_T1: ', SITEP(1,2),SITEP(2,2),SITEP(3,2),R2_T1 ! ! Start the East, West, North, South offsets loop DO N = 1, 4 ! E-W-N-S loop ! VdotNFO = DOTP(VsubEarth,SpEWNS(1,N)) ! ! Initial estimate for (T1 - T0) !RR CALL VECSB (SPpxyz, SITEP(1,1), R_01) CALL VECSB (SpEWNS(1,N), SITEP(1,1), R_01) R_01mag = VECMG(R_01) del_T0 = VECMG(R_01)/VLIGHT ! write(6,*) 'DUEV: T0_T1, Initial del_T0 ', T0_T1, del_T0 ! gm1 = 2.D0*GMSUN/VLIGHT2 ! ! Three iterations should be enough DO J = 1,3 Do I = 1,3 ! Barycentric at T1 - del_T0 R0_T0(I) = SpEWNS(I,N)*F1 - 0.5D0*(VdotNFO/VLIGHT2)*VsubEarth(I) & & + (EARTH(I,1) - del_T0*EARTH(I,2)) - & & SPvxyz(I)*(T0_T1 - del_T0) Enddo ! ! Barycentric vectors Call VECSB(R1_T1, R0_T0, R_01) ! equation 15 R_01mag = VECMG(R_01) !!! write(6,*) 'R_01,R_01mag: ', R_01,R_01mag ! CALL VECMU(SUNb(1,2),del_T0,delSun0) CALL VECSB(SUNb(1,1), delSun0, RsunT0) ! write(6,*) 'new RsunT0: ', RsunT0 CALL VECSB (R0_T0, RsunT0, R0sunT0) ! eqn 16, Sun to spacecraft at T0 R0sunT0mag = VECMG(R0sunT0) !!! write(6,*) 'R0sunT0,R0sunT0mag: ', R0sunT0,R0sunT0mag ! CALL VECMU(MOONb(1,2),del_T0,delMoon0) CALL VECSB(MOONb(1,1), delMoon0, RmoonT0) ! write(6,*) 'new RmoonT0: ', RmoonT0 CALL VECSB (R0_T0, RmoonT0, R0moonT0) ! eqn 16, Moon to spacecraft at T0 R0moonT0mag = VECMG(R0moonT0) ! CALL VECMU(EARTH(1,2),del_T0,delEarth0) CALL VECSB(EARTH(1,1), delMoon0, REarthT0) CALL VECSB (R0_T0, RearthT0, R0earthT0) ! eqn 16, Earth to spacecraft at T0 R0earthT0mag = VECMG(R0earthT0) ! write(6,*) 'RearthT0: ', RearthT0 ! write(6,*) 'R0earthT0: ', R0earthT0 ! write(6,*) 'R0earthT0mag: ', R0earthT0mag CALL VECEQ(EARTH(1,1),RearthT1) ! write(6,*) 'RearthT1: ', RearthT1 CALL VECSB(R1_T1, RearthT1, R1earthT1) ! eqn 16, Earth to site 1 at T1 R1earthT1mag = VECMG(R1earthT1) ! write(6,*) 'R1earthT1: ', R1earthT1 ! write(6,*) 'R1earthT1mag: ', R1earthT1mag CALL VECSB (R1earthT1,R0earthT0,R01earth) ! eqn 17, Earth R01earthmag = VECMG(R01earth) ! write(6,*) 'R01earth: ', R01earth ! write(6,*) 'R01earthmag: ', R01earthmag ! CALL VECEQ(SUNb(1,1),RsunT1) ! write(6,*) 'new RsunT1: ', RsunT1 CALL VECSB(R1_T1, RsunT1, R1sunT1) ! eqn 16, Sun to site 1 at T1 R1sunT1mag = VECMG(R1sunT1) !!! write(6,*) 'R1sunT1,R1sunT1mag: ', R1sunT1,R1sunT1mag ! CALL VECEQ(MOONb(1,1),RmoonT1) ! write(6,*) 'new RmoonT1: ', RmoonT1 CALL VECSB(R1_T1, RmoonT1, R1moonT1) ! eqn 16, Moon to site 1 at T1 R1moonT1mag = VECMG(R1moonT1) ! CALL VECSB (R1sunT1,R0sunT0,R01sun) ! eqn 17, Sun R01sunmag = VECMG(R01sun) !!! write(6,*) 'R01sun,R01sunmag: ', R01sun,R01sunmag CALL VECSB (R1moonT1,R0moonT0,R01moon) ! eqn 17, Moon R01moonmag = VECMG(R01moon) ! Sun: RLT_01s = (2.D0*GMSUN/VLIGHT3) * DLOG( (R0sunT0mag + R1sunT1mag & & + R01sunmag + gm1)/(R0sunT0mag + R1sunT1mag - R01sunmag + gm1) ) ! Moon: RLT_01m = (2.D0*GMMOON/VLIGHT3) * DLOG( (R0moonT0mag + R1moonT1mag & & + R01moonmag)/(R0moonT0mag + R1moonT1mag - R01moonmag) ) ! write(6,*) 'RLT_01s, RLT_01m ', RLT_01s,RLT_01m ! ! Earth: ! RLT_01e = (2.D0*GMEarth/VLIGHT3) * DLOG( (R0earthT0mag + R1earthT1mag & ! & + R01earthmag)/(R0earthT0mag + R1earthT1mag - R01earthmag) ) ! write(6,*) 'RLT_01s ', RLT_01s ! write(6,*) 'RLT_01m ', RLT_01m ! write(6,*) 'RLT_01e ', RLT_01e ! ! Planets: RLT_01p = 0.0D0 ! write(6,*) 'J, del_T0', J,del_T0 Do K = 1,7 Call VECMU(SPLANET(1,2,K),del_T0,delpl) CALL VECSB(SPLANET(1,1,K), delpl, Spl_0) ! SSBC planet vector at T0 CALL VECSB(R0_T0, Spl_0, R0planT0) ! Planet to spacecraft at T0 R0planT0mag = VECMG(R0planT0) ! CALL VECSB(R1_T1,SPLANET(1,1,K), R1planT1) ! Planet to site 1 at T1 R1planT1mag = VECMG(R1planT1) ! CALL VECSB(R1planT1, R0planT0, R01plan) ! ~Planet-to-site1 R01planmag = VECMG(R01plan) ! RLTp(K) = (2.D0*GMPLANET(K)/VLIGHT3) * & & DLOG( (R0planT0mag + R1planT1mag + R01planmag) / & & (R0planT0mag + R1planT1mag - R01planmag) ) RLT_01p = RLT_01p + RLTp(K) Enddo ! RLT_01 = RLT_01s + RLT_01m + RLT_01p T1_T0 = R_01mag/VLIGHT + RLT_01 ddel_T0 = T1_T0 - del_T0 del_T0 = T1_T0 ! write(6,*) 'J,RLT_01, del_T0: ', J, RLT_01, del_T0 ! Enddo ! ! Now find (T2 - T0): ! Initial estimate for (T2 - T0) Call VECSB(R2_T1, R0_T0, R_02) T2_T0 = VECMG(R_02)/VLIGHT del_T2 = T1_T0 - T2_T0 ! difference from T1 ! write(6,*) 'R_02,T2_T0,del_T2: ', R_02,T2_T0,del_T2 ! DO J = 1,3 ! Three iterations ! Barycentric TDB vectors: Do I = 1,3 Site2(I) = SITEP(I,2) - SITEV(I,2)*del_T2 ! Site 2 at T2 VdotR2 = DOTP(VsubEarth,Site2) R2_T2(I) = Site2(I)*F1 - 0.5D0*(VdotR2/VLIGHT2)*VsubEarth(I) & & + (EARTH(I,1) - del_T2*EARTH(I,2)) ! SSBC Site2 vector @ T2 Enddo R2_T2mag = VECMG(R2_T2) ! write(6,*) 'SITEP2,Site2:', SITEP(1,2),SITEP(2,2),SITEP(3,2),Site2 ! write(6,*) 'R2_T2,R2_T2mag:', R2_T2,R2_T2mag ! Call VECSB(R2_T2, R0_T0, R_02) ! eqn 15; Site2 to NFO R_02mag = VECMG(R_02) ! CALL VECMU(SUNb(1,2),del_T2,delSun2) CALL VECSB(SUNb(1,1), delSun2, RsunT2) CALL VECSB (R2_T2, RsunT2, R2sunT2) ! eqn 16, Sun to site 2 at T2 R2sunT2mag = VECMG(R2sunT2) ! write(6,*) 'new RsunT2: ', RsunT2 ! write(6,*) 'new R2sunT2mag: ', R2sunT2mag CALL VECSB (R2sunT2,R0sunT0,R02sun) ! eqn 17, Sun R02sunmag = VECMG(R02sun) ! write(6,*) 'new R02sun,R02sunmag: ', R02sun,R02sunmag ! CALL VECMU(MOONb(1,2),del_T2,delmoon2) CALL VECSB(MOONb(1,1), delmoon2, RmoonT2) CALL VECSB (R2_T2, RmoonT2, R2moonT2) ! eqn 16, Moon to site 2 at T2 R2moonT2mag = VECMG(R2moonT2) ! write(6,*) 'new R2moonT2: ', R2moonT2 ! write(6,*) 'new R2moonT2mag: ', R2moonT2mag CALL VECSB (R2moonT2,R0moonT0,R02moon) ! eqn 17, Moon R02moonmag = VECMG(R02moon) ! write(6,*) 'new R02moon,R02moonmag: ', R02moon,R02moonmag !! CALL VECMU(EARTH(1,2),del_T2,delearth2) CALL VECSB(EARTH(1,1), delearth2, RearthT2) !! write(6,*) 'new RearthT2: ', earthRT2 CALL VECSB (R2_T2, RearthT2, R2earthT2) ! eqn 16, Earth to site 2 at T2 R2earthT2mag = VECMG(R2earthT2) ! write(6,*) 'R2earthT2mag: ', R2earthT2mag CALL VECSB (R2earthT2,R0earthT0,R02earth) ! eqn 17, Earth R02earthmag = VECMG(R02earth) !! write(6,*) 'R02earth: ', R02earth ! write(6,*) 'R02earthmag: ', R02earthmag ! RLT_02s = (2.D0*GMSUN/VLIGHT3) * DLOG((R0sunT0mag + R2sunT2mag & & + R02sunmag + gm1)/(R0sunT0mag + R2sunT2mag - R02sunmag + gm1) ) RLT_02m = (2.D0*GMMOON/VLIGHT3) * DLOG((R0moonT0mag + R2moonT2mag & & + R02moonmag)/(R0moonT0mag + R2moonT2mag - R02moonmag) ) ! write(6,*) 'RLT_02s, RLT_02m', RLT_02s, RLT_02m RLT_02e = (2.D0*GMEarth/VLIGHT3) * DLOG( (R0earthT0mag + R2earthT2mag & & + R02earthmag)/(R0earthT0mag + R2earthT2mag - R02earthmag) ) ! write(6,*) 'RLT_02s ', RLT_02s ! write(6,*) 'RLT_02m ', RLT_02m ! write(6,*) 'RLT_02e ', RLT_02e ! write(6,*) ' ' ! ! Planets: RLT_02p = 0.0D0 Do K = 1,7 Call VECMU(SPLANET(1,2,K),del_T0,delpl) CALL VECSB(SPLANET(1,1,K), delpl, Spl_0) ! SSBC planet vector at T0 CALL VECSB(R0_T0, Spl_0, R0planT0) ! Planet to spacecraft at T0 R0planT0mag = VECMG(R0planT0) ! Call VECMU(SPLANET(1,2,K),del_T2,delpl) CALL VECSB(SPLANET(1,1,K), delpl, Spl_2) ! SSBC planet vector at T2 CALL VECSB(R2_T2,Spl_2, R2planT2) ! Planet to site 2 at T2 R2planT2mag = VECMG(R2planT2) ! CALL VECSB(R2planT2, R0planT0, R02plan) ! ~Planet-to-site2 R02planmag = VECMG(R02plan) ! RLTp(K) = (2.D0*GMPLANET(K)/VLIGHT3) * & & DLOG( (R0planT0mag + R2planT2mag + R02planmag) / & & (R0planT0mag + R2planT2mag - R02planmag) ) RLT_02p = RLT_02p + RLTp(K) Enddo ! write(6,*) 'RLT_02p ', RLT_02p ! RLT_02 = RLT_02s + RLT_02m + RLT_02e + RLT_02p ! write(6,*) 'RLT_02 ', RLT_02 T2_T0 = R_02mag/VLIGHT + RLT_02 ddel_T2 = T2_T0 - (T1_T0 - del_T2) del_T2 = T1_T0 - T2_T0 ! ! write (6,*) 'del_T2, T2-T1: ', del_T2, (T2_T0 - T1_T0) Enddo ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Convert to TT frame, Duev equation 20. ! T2_T1 = T2_T0 - T1_T0 ! Vacuum delay in TDB frame VEsq = DOTP(VsubEarth,VsubEarth) ! Square of Earth's SSBC velocity UE = GMSUN/vecmg1 ! Newtonian potential, Sun only Call VECSB(SITEP(1,2),SITEP(1,1),B) ! Geocentric baseline vector VEdotB = DOTP(VsubEarth,B) ! Earth SSB velocity * Baseline VEdotV2 = DOTP(VsubEarth,SITEV(1,2)) ! Earth SSB velocity * Site2 geocentric velocity ! td2_td1(N) = (( (T2_T1/(1.D0-1.48082686741D-8)) * & & (( 1.D0 - (VEsq/2.D0 + UE)/VLIGHT2 )) - VEdotB/VLIGHT2 ) / & & (1.D0 + VEdotV2/VLIGHT2)) ! If (UVW .eq. 'exact ') Then td2_td1(N) = td2_td1(N) + Datmc_h_EWNS(N) + Datmc_w_EWNS(N) ! write(6,'("I,td2_td1(psec):",I2,F18.2)') I, td2_td1(I)*1.D12 Endif ! write (6,*) 'DUEVewns/new: N, td2_td1(ps): ', N, td2_td1(N)*1.D12 ! ENDDO ! E-W-N-S loop ! U_V(1) = VLIGHT*(td2_td1(1) - td2_td1(2))/0.0002D0 U_V(2) = VLIGHT*(td2_td1(3) - td2_td1(4))/0.0002D0 ! write(6,*) 'Duev New: U (m, psec) ', U_V(1), U_V(1)*1.D12/VLIGHT ! write(6,*) 'Duev New: V (m, psec) ', U_V(2), U_V(2)*1.D12/VLIGHT ! RETURN END ! !********************************************************************** Subroutine RANGEewns (Datmc_h_EWNS, Datmc_w_EWNS, SITEP, SITEV, & & SITEA, SUN, T0_T1, & & R1, R1dt, R1mag, R1magdt, R2, R2dt, R2mag, R2magdt) ! & tr2_tr1, dtr2_tr1) ! IMPLICIT NONE ! REAL*8 SUN(3,2), EARTH(3,3), SITEA(3,2), SITEP(3,2), & & SITEV(3,2), DOTP, VECMG, Datmc_h_EWNS(4), Datmc_w_EWNS(4) Real*8 R1(3), R2(3), R1dt(3), R2dt(3), R1mag, R2mag, R1magdt, & & R2magdt, T0_T1, tr2_tr1(4), dtr2_tr1(4) Real*8 SunNFO(3), del2(3), Site2(3),x2x0(3), rsun2(3), vecx2x0(3) Real*8 G2sunV3, rsun0mag, rsun1mag, rsun2mag, x1x0mag, x2x0mag, & & xlog1, xlog2, t1t0, t2t0, t1t2, t1t2_old, G2EarthV3, & & xEx0mag, Site2mag, tgrav, tgravSun, tgravEarth Integer*4 J, N ! INCLUDE 'd_input.i' ! Variables from: ! INCLUDE 'cphys11.i' ! Variables 'from': ! 1. VLIGHT - The velocity of light in vacuum. (m/sec) ! 2. VLIGHT2 - The velocity of light squared. (m/sec)**2 ! 3. VLIGHT3 - THE VELOCITY OF LIGHT IN VACUUM CUBED. ! (M**3/SEC**3) ! 4. GMSUN, GMMOON, GMEARTH, GMPLANET(7) - Gravitational constant ! times the masses of the Sun, Moon, Earth, and the ! other planets except Pluto. ! INCLUDE 'put2s.i' ! Variables from: ! Real*8 SpEWNS(3,4), R1EWNS(3,4), R2EWNS(3,4), R1EWNSmag(4), & & R2EWNSmag(4), STAR12EWNS(3,4,2), STARewns(3,4) Common / NFewns/ SpEWNS, R1EWNS, R2EWNS, R1EWNSmag, R2EWNSmag, & & STAR12EWNS, STARewns ! DO N = 1, 4 ! loop over East, West, North, South offsets ! G2sunV3 = 2.D0*GMSUN/VLIGHT3 Call VECSB(SpEWNS(1,N),SUN(1,1),SunNFO) rsun0mag = VECMG(SunNFO) rsun1mag = VECMG(SUN(1,1)) x1x0mag = VECMG(R1) xlog1 = rsun0mag + rsun1mag + x1x0mag xlog2 = rsun0mag + rsun1mag - x1x0mag tgravSun = G2sunV3 * DLOG(xlog1/xlog2) ! G2EarthV3 = 2.D0*GMEARTH/VLIGHT3 xEx0mag = VECMG(SpEWNS(1,N)) ! xlog1 = xEx0mag + xEx0mag ! xlog2 = xEx0mag - xEx0mag ! tgravEarth = G2EarthV3 * DLOG(xlog1/xlog2) tgravEarth = 0.D0 ! zero at geocenter ! ! Delay from NFO to site1 (geocenter) t1t0 = R1mag/VLIGHT + tgravEarth ! t1t2 = 0.D0 Do J = 1,4 t1t2_old = t1t2 Call VECMU(SITEV(1,2),t1t2,del2) ! del2 = vector correction to site2 position Call VECSB(SITEP(1,2),del2,Site2) ! Site2 = site2 retarded position Call VECSB(Site2,SpEWNS(1,N),vecx2x0) x2x0mag = VECMG(vecx2x0) Call VECSB(Site2,SUN(1,1),rsun2) rsun2mag = VECMG(rsun2) ! xlog1 = rsun0mag + rsun2mag + x2x0mag xlog2 = rsun0mag + rsun2mag - x2x0mag tgravSun = G2sunV3 * DLOG(xlog1/xlog2) ! Site2mag = VECMG(Site2) xlog1 = xEx0mag + Site2mag + x2x0mag xlog2 = xEx0mag + Site2mag - x2x0mag tgravEarth = G2EarthV3 * DLOG(xlog1/xlog2) ! t2t0 = x2x0mag/VLIGHT + tgravSun + tgravEarth t2t0 = x2x0mag/VLIGHT + tgravEarth ! t2t0 = x2x0mag/VLIGHT t1t2 = t1t0 - t2t0 Enddo ! tr2_tr1(N) = t2t0 - t1t0 ! If (UVW .eq. 'exact ') Then tr2_tr1(N) = tr2_tr1(N) + Datmc_h_EWNS(N) + Datmc_w_EWNS(N) Endif ! ENDDO ! end loop over East, West, North, South offsets ! ! write (6,*) ' Ranging model ' ! write (6,*) ' tr2_tr1east ', tr2_tr1(1) ! write (6,*) ' tr2_tr1west ', tr2_tr1(2) ! write (6,*) ' tr2_tr1north ', tr2_tr1(3) ! write (6,*) ' tr2_tr1south ', tr2_tr1(4) ! U_V(1) = VLIGHT*(tr2_tr1(1) - tr2_tr1(2))/0.0002D0 U_V(2) = VLIGHT*(tr2_tr1(3) - tr2_tr1(4))/0.0002D0 ! write(6,*) ' U (m, psec) ', U_V(1), U_V(1)*1.D12/VLIGHT ! write(6,*) ' V (m, psec) ', U_V(2), U_V(2)*1.D12/VLIGHT ! Return End !**********************************************************************