Updates to PositionalAstronomy.f90

	Updates to PositionalAstronomy.f90

PositionalAstronomy.f90

@@ -374,4 +374,104 @@ SUBROUTINE HELIOBARY_CENTRIC_CORRECTION(MN_RA,MN_DEC,EQX,TMJD,CORHV,CORHZ,&
    
 END SUBROUTINE HELIOBARY_CENTRIC_CORRECTION
 
+! -------------------------------------------------------------
+! -------------------------------------------------------------
+
+SUBROUTINE HELIOBARY_CENTRIC_VELOCITY(MN_RA,MN_DEC,EQX,TMJD,DIRHV,DIRBV,HAFLG)
+
+   ! --------
+   ! Purpose:
+   ! --------
+   ! Given the mean RA and DEC of a target and a date of observation
+   ! in MJD per specified Equinox, determine the velocity in that direction
+   ! (in km/s)  due to the Earths motion around the sun. Note the Earths
+   ! velocity around the Sun is roughly 30km/s
+   
+   ! -----------
+   ! Description:
+   ! -----------
+   ! Derivation of the correction velocity is textbook as taken from the 
+   ! slalib guide and uses the following
+   ! salib routines:
+   ! SLAEVP      ~ Calculate the Earths position velocity and velocity vector
+   !               relative to the sun (DPH and DPV) and the Earth Sun
+   !               barycentre (DPB and DVB) from given date/time. Note units
+   !               are in AU and AU/s
+   !
+   ! SLAEVP      ~ Same as SLAEVP but listed as deriving higher accuracy.Note
+   !               units are in AU and AU/d so return velocity values need
+   !               to be divided by 86400
+   !
+   ! SLADCS2C    ~ Calculate the direction cosines (unit direction vector)
+   !               of a target specified in RA and Dec.
+   !
+   ! SLADVDV     ~ Perform a double precision scalar product of two vectors
+   !
+   ! Correction z value is simply (correction velocity)/(speed of light)
+   !
+   
+   ! ----------
+   ! Arguments:
+   ! ----------
+   DOUBLE PRECISION, INTENT(IN)  :: MN_RA   ! Mean Right Ascension (deg)
+   DOUBLE PRECISION, INTENT(IN)  :: MN_DEC  ! Mean Declination     (deg)
+   DOUBLE PRECISION, INTENT(IN)  :: EQX     ! Chosen Refernce Equinox
+   DOUBLE PRECISION, INTENT(IN)  :: TMJD    ! Time in MJD
+   DOUBLE PRECISION, INTENT(OUT) :: DIRHV   ! Return helio-centric velocity
+   DOUBLE PRECISION, INTENT(OUT) :: DIRBV   ! Return bary- centric velocity
+   LOGICAL,          INTENT(IN)  :: HAFLG   ! Flag to use high accuracy
+                                            ! slalib routines
+   ! ----------------
+   ! Local Variables:
+   ! ----------------
+   DOUBLE PRECISION ::DVB(3), DPB(3), DVH(3), DPH(3), UNIT_V(3)
+   DOUBLE PRECISION ::VEL_INAU
+   DOUBLE PRECISION ::VEL_INKM
+   
+   ! ----------
+   ! Constants:
+   ! ----------
+   DOUBLE PRECISION, PARAMETER :: AU_dp  = 149.597870E6 ! Km
+   DOUBLE PRECISION, PARAMETER ::  c_dp  =        3.0E5 ! Km/s (= 3.0E8 m/s)
+
+
+   ! Get the position vector and the velocity vector of the Earth relative
+   ! to the sun (DPH and DVH) at the specified date time (TMJD wrt EQX)
+   ! Note return vectors are in units of AU and AU/s
+   IF (.NOT.HAFLG) THEN
+      ! Use the standard SLAEVP slalib routine
+      CALL SLAEVP ( TMJD, EQX, DVB, DPB, DVH, DPH )
+   ELSE
+      ! Use the "high accuracy" SLAEPV slalib routine that returns velocites
+      ! in AU/day
+      CALL SLAEPV ( TMJD,      DPH, DVH, DPB, DVB )
+      DVH=DVH/86400;DVB=DVB/86400;
+   ENDIF
+
+   ! Get the unit direction vector from the Earth to the specified target at
+   ! spherical coordinates in radians (RA, DEC)
+   ! Note: due to the lack of parallax, this direction vector is constant
+   ! from the Earth at any point along it's path of orbit around the Sun.
+   CALL SLADCS2C ( MN_RA/DR2D, MN_DEC/DR2D, UNIT_V )
+   
+   ! Velocity of Earth towards target is simply the scalar product of the
+   ! velocity vector with the unit direction vector. Note this velocity is in
+   ! AU/s
+   VEL_INAU = SLADVDV ( UNIT_V, DVH ) 
+   
+   ! Convert to velocity in km/s
+   VEL_INKM = VEL_INAU * AU_dp
+   
+   ! Thus required velocity (in km/s) value to return
+   DIRHV = VEL_INKM 
+   
+   
+   ! And the same calculations for the bary-centric corrections
+   VEL_INAU = SLADVDV ( UNIT_V, DVB ) 
+   VEL_INKM = VEL_INAU * AU_dp
+   DIRBV = VEL_INKM 
+   
+   
+END SUBROUTINE HELIOBARY_CENTRIC_VELOCITY
+
 END MODULE POSITIONAL_ASTRONOMY