! riseset.f90
! Written by Nicholas Moe, completed 23 April 2007
! Computes the sunrise and sunset times for an observer, given the location and date.
! Follows algorithms and procedures given by Peter Duffett-Smith, Practical astronomy
! with your calculator, Second Edition. Cambridge, 1981.

! Convention: West longitudes are positive

SUBROUTINE sun_coords(day, month, year, lambda_sun, alpha_sun, delta_sun)
  ! This subroutine takes the day, month, and year, and gives the ecliptic longitude
  ! of the sun, and converts this coordinate to equatorial coordinates. It follows
  ! the pseudocode given by Duffett-Smith, section 43.
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  REAL(KIND=double), INTENT(IN) :: day         ! The day of the month
  INTEGER, INTENT(IN) :: month    ! The month of the year (expressed as an integer)
  INTEGER, INTENT(IN) :: year     ! The year
  REAL(KIND=double), INTENT(OUT) :: lambda_sun  ! Ecliptic Longitude of the sun (degrees)
  REAL(KIND=double), INTENT(OUT) :: alpha_sun  ! The right ascension in decimal form
  REAL(KIND=double), INTENT(OUT) :: delta_sun  ! The declination in decimal form
  REAL(KIND=double), EXTERNAL :: julian_date
  REAL(KIND=double) :: D          ! Number of days since 1980.0
  REAL(KIND=double) :: N          ! N = 360 / 365.2422 * D (in degrees); intermediate
                                  ! step for getting M, the mean anomaly
  REAL(KIND=double) :: M          ! Mean anomaly (degrees)
  REAL(KIND=double), PARAMETER :: e_g = 278.83354_double   ! Sun's ecliptic longitude at epoch 1980.0 (degrees)
  REAL(KIND=double), PARAMETER :: w_g = 282.596403_double  ! Sun's ecliptic longitude of perigree (degrees)
  REAL(KIND=double), PARAMETER :: pi = 3.141592654D0
  REAL(KIND=double) :: E           ! Eccentric anomaly
  REAL(KIND=double), PARAMETER :: epsilon = (1._double/(10._double**6))  ! Required accuracy
  REAL(KIND=double), PARAMETER :: ecc = 0.016718_double  ! Eccentricity of orbit
  REAL(KIND=double), EXTERNAL :: kepler_sol
  REAL(KIND=double) :: nu          ! True anomaly (degrees)
  REAL(KIND=double) :: beta_sun = 0.D0 ! Ecliptic latitude (decimal degrees)

!jc  Nicely commented header.

  ! Find the number of days since the beginning of the year in question.
  ! To simplify things (like the concern for leap years), we will first convert
  ! the date in question to the Julian date. Then we will calculate the number of days elapsed
  ! since 1980.0, which is January 0.0 1980. January 0.0 is midnight between 30 Dec and 31 Jan
  ! of the previous year.
  D = julian_date(day, month, year) - julian_date(0._double, 1, 1980)
  !PRINT *, 'D', D
  N = (360._double / 365.2422_double) * D
  !PRINT *, 'N', N
  
  ! Add or subtract 360 degrees until 0 <= N <= 360 degrees
!jc  You could probably do this more cleanly using MOD ...
  IF (N > 360._double) THEN
    DO
      N = N - 360._double
      IF (N < 360._double) EXIT
    END DO
  ELSE IF (N < 0._double) THEN
    DO
      N = N + 360._double
      IF (N > 0._double) EXIT
    END DO
  END IF

  !PRINT *, 'N', N
   
  M = N + e_g - w_g
  
  IF (M < 0._double) M = M + 360._double
  
  !PRINT *, 'M', M
  
  ! Convert M to radians:
  M = M * (pi / 180._double)
  
  !PRINT *, 'M (radians)', M
  
  E = kepler_sol(M)
  
  !PRINT *, 'E', E
  
  nu = 2. * (180._double/pi) * ATAN( SQRT((1._double + ecc)/(1._double - ecc)) * TAN(E/2._double))
  
  !PRINT *, 'nu', nu
  
  lambda_sun = nu + w_g
  
  IF (lambda_sun > 360._double) lambda_sun = lambda_sun - 360._double
  IF (lambda_sun < 0._double) lambda_sun = lambda_sun + 360._double
  
!   PRINT *, 'lambda_sun', lambda_sun
!   PRINT *, 'beta_sun', beta_sun
  
  CALL convert_to_ra_dec(lambda_sun, beta_sun, alpha_sun, delta_sun)
  
!   PRINT *, 'lambda_sun', lambda_sun
!   PRINT *, 'R.A.', alpha_sun
!   PRINT *, 'Dec', delta_sun
  
  RETURN
END SUBROUTINE sun_coords

FUNCTION julian_date(day, month, year)
  ! Purpose: given the day, month, and year, output the Julian date. Follows pseudocode
  ! given by Duffett-Smith, section 4.
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  REAL(KIND=double) :: julian_date
  REAL(KIND=double), INTENT(IN) :: day         ! The day of the month, which may be in decimal form
  INTEGER, INTENT(IN) :: month    ! The month of the year (expressed as an integer)
  INTEGER, INTENT(IN) :: year     ! The year
  INTEGER :: year_prime           ! Altered year (y')
  INTEGER :: month_prime          ! Altered month (m")
  INTEGER :: A                    ! Integer part of (y'/100)
  INTEGER :: B                    ! B = 2 - A + integer part of (A/4), or 0
  INTEGER :: C                    ! C = integer part of (365.25 * y')
  INTEGER :: D                    ! D = integer part of (30.6001 * (m' + 1))
  
  IF (month == 1) THEN
    year_prime = year - 1
    month_prime = month + 12
  ELSE IF (month == 2) THEN
    year_prime = year - 1
    month_prime = month + 12
  ELSE
    year_prime = year
    month_prime = month
  END IF
  
  ! Do different things depending on whether the date is later than 1582 October 15 or not  
  IF (year > 1582) THEN
    A = INT(REAL(year_prime)/100.)
    B = 2 - A + INT(REAL(A)/4.)
  ELSE IF (year < 1582) THEN
    B = 0
  ELSE IF (month < 10) THEN ! If the program gets to this statement, then the year is 1582.
    B = 0
  ELSE IF (month > 10) THEN ! Now it must either be 1582 Oct, Nov, or Dec
    A = INT(REAL(year_prime)/100.)    ! Do this if it is 
    B = 2 - A + INT(REAL(A)/4.)
  ELSE IF (day > 15._double) THEN ! Now it must be some day in October 1582
    A = INT(REAL(year_prime)/100.)
    B = 2 - A + INT(REAL(A)/4.)
  ELSE
    B = 0
  END IF
  
  C = INT(365.25 * year_prime)
  
  D = INT(30.6001 * (month_prime + 1))
  
  julian_date = B + C + D + day + 1720994.5_double
  
  RETURN
END FUNCTION julian_date

FUNCTION kepler_sol(M)
  ! This function finds a solution to the Kepler equation, E - e sin E = M, using
  ! a routine described in Duffett-Smith, p. 85.
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  REAL(KIND=double) :: kepler_sol
  REAL(KIND=double), INTENT(IN) :: M
  REAL(KIND=double) :: E               ! Eccentric anomaly
  REAL(KIND=double), PARAMETER :: epsilon = (1._double/(10._double**6))  ! Required accuracy
  REAL(KIND=double), PARAMETER :: ecc = 0.016718_double  ! Eccentricity of orbit
  REAL(KIND=double) :: delta           ! Dummy variable
  REAL(KIND=double) :: delta_E         ! Dummy variable
  
  E = M
  DO
    delta = E - (ecc * SIN(E)) - M
    IF (delta <= epsilon) EXIT
    delta_E = delta / (1._double - (ecc * COS(E)))
    E = E - delta_E
  END DO
  
  kepler_sol = E
  RETURN
END FUNCTION kepler_sol

SUBROUTINE convert_to_ra_dec(lambda_deg, beta_deg, ra, dec)
  ! This subroutine takes lambda (ecliptic longitude) and beta (ecliptic latitude) in
  ! degrees (in decimal form) and converts it to right ascension in decimal hours and
  ! declination in decimal degrees. It follows pseudocode given in Duffett-Smith,
  ! section 27.
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  REAL(KIND=double), INTENT(IN) :: lambda_deg ! Ecliptic longitude (decimal degrees)
  REAL(KIND=double), INTENT(IN) :: beta_deg   ! Ecliptic latitude (decimal degrees)
  REAL(KIND=double) :: lambda, beta           ! Radian versions of the above two
  REAL(KIND=double), INTENT(OUT) :: ra        ! Right ascension (decimal hours)
  REAL(KIND=double), INTENT(OUT) :: dec       ! Declination (decimal hours)
  REAL(KIND=double), PARAMETER :: pi = 3.141592654D0
  REAL(KIND=double), PARAMETER :: pi180 = (pi / 180._double) ! Multiply angle variable in degrees by this
                                                             ! to convert to radians
  REAL(KIND=double) :: delta                  ! Dummy variable for declination (in degrees)
  REAL(KIND=double) :: epsilon = 23.441884_double ! Obliquity of the ecliptic at epoch 1980.0 (degrees)
  REAL(KIND=double) :: y                      ! Dummy variable
  REAL(KIND=double) :: x                      ! Dummy variable
  REAL(KIND=double) :: alpha                  ! R.A. in degree form
  
  ! Convert lambda and beta to radians
  lambda = lambda_deg * pi180
  beta = beta_deg * pi180
  epsilon = 23.441884_double * pi180
  
  delta = ASIN(SIN(beta) * COS(epsilon) + COS(beta) * SIN(epsilon) * SIN(lambda)) * (180._double / pi)
  y = SIN(lambda) * COS(epsilon) - TAN(beta) * SIN(epsilon)
  x = COS(lambda)
  
!  PRINT *, 'y/x', y/x
  
  alpha = ATAN2(y, x) * (180._double / pi) ! ATAN2 gives the correct quadrant for (y/x)
!  PRINT *, 'alpha', alpha
  
  ra = alpha / 15._double
  dec = delta
  
  RETURN
END SUBROUTINE convert_to_ra_dec

SUBROUTINE LST_to_GST(LST, longitude, GST)
  ! This subroutine converts local sidereal time (LST) to Greenwich sidereal time.
  ! It follows pseudocode given in Duffett-Smith, section 15.
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  REAL(KIND=double), INTENT(IN) :: LST    ! Local sidereal time in decimal hours
  REAL(KIND=double), INTENT(IN) :: longitude    ! Local longitude, in decimal degrees
  REAL(KIND=double), INTENT(OUT) :: GST   ! Greenwich sidereal time in decimal hours
  REAL(KIND=double) :: longitude_h              ! Local longitude, in decimal hours
  
  longitude_h = longitude / 15.D0
  GST = LST + longitude_h
END SUBROUTINE LST_to_GST

SUBROUTINE GST_to_GMT(GST, day, month, year, GMT)
  ! This subroutine converts GST to GMT. It follows pseudocode given in Duffett-Smith,
  ! section 13.
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  REAL(KIND=double), INTENT(IN) :: GST
  REAL(KIND=double), INTENT(IN) :: day
  REAL(KIND=double), INTENT(IN) :: month
  REAL(KIND=double), INTENT(IN) :: year
  REAL(KIND=double), INTENT(OUT) :: GMT
  REAL(KIND=double), EXTERNAL :: constant_B
  REAL(KIND=double), EXTERNAL :: julian_date
  REAL(KIND=double), PARAMETER :: A = 0.0657098D0
  REAL(KIND=double), PARAMETER :: C = 1.002738D0
  REAL(KIND=double), PARAMETER :: D = 0.997270D0
  REAL(KIND=double) :: T0           ! Dummy variable
  
  T0 = ((julian_date(day, month, year) - julian_date(0.D0, 1, year)) * A) - constant_B(year)
  !PRINT *, 'T0 before convert', T0
  DO
    IF (T0 >= 0.D0) EXIT  
    T0 = T0 + 24.D0
  END DO
  !PRINT *, 'T0 after convert', T0
  
  IF ((GST - T0) < 24.D0) THEN
    GMT = (GST - T0 + 24.D0) * D
  ELSE
    GMT = (GST - T0) * D
  END IF
  
END SUBROUTINE GST_to_GMT

FUNCTION constant_B(year)
  ! Purpose: given the year, calculate constant B. This is used for the conversion
  ! from GST to GMT. Follows a routine given by Duffett-Smith in section 12.
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  INTEGER, INTENT(IN) :: year
  REAL(KIND=double) :: constant_B
  REAL(KIND=double), EXTERNAL :: julian_date
  REAL(KIND=double) :: JD     ! Julian date on 0.0 January of the year
  REAL(KIND=double) :: S      ! Dummy variable
  REAL(KIND=double) :: T      ! Dummy variable
  REAL(KIND=double) :: R      ! Dummy variable
  REAL(KIND=double) :: U      ! Dummy variable
  
  JD = julian_date(0.D0, 1, year)
  S = JD - 2415020.D0
  T = S / 36525.D0
  R = 6.6460656D0 + (2400.051262D0 * T) + (0.00002581D0 * (T*T))
  U = R - (24.D0 * REAL( (year - 1900),double) )
  constant_B = 24.D0 - U
  
END FUNCTION constant_B
  
SUBROUTINE rough_riseset(alpha, delta, phi, longitude, day, month, year, LSTr, LSTs)
  ! This subroutine calculates the local rising and setting times of an celestial object
  ! based on its RA (alpha), declination (delta), and the observer's geographical latitude (phi). It follows pseudocode given in Duffett-Smith, section 32.
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  REAL(KIND=double), INTENT(IN) :: alpha      ! Right Ascension (decimal hours)
  REAL(KIND=double), INTENT(IN) :: delta      ! Declination (decimal degrees)
  REAL(KIND=double), INTENT(IN) :: phi        ! Geographical latitude (decimal degrees)
  REAL(KIND=double), INTENT(IN) :: longitude  ! Geographical longitude (decimal degrees)
  REAL(KIND=double), INTENT(IN) :: day
  REAL(KIND=double), INTENT(IN) :: month
  REAL(KIND=double), INTENT(IN) :: year
  REAL(KIND=double) :: GMTr
  REAL(KIND=double) :: GMTs
  REAL(KIND=double), PARAMETER :: pi = 3.141592654D0
  REAL(KIND=double), PARAMETER :: pi180 = (pi / 180.D0) ! Multiply angle variable in degrees by this
  REAL(KIND=double) :: Ar                     ! Azimuth at rising
  REAL(KIND=double) :: As                     ! Azimuth at setting
  REAL(KIND=double) :: H                      ! Dummy variable (Hours)
  REAL(KIND=double), INTENT(OUT)  :: LSTr                   ! Local sidereal time for rising
  REAL(KIND=double), INTENT(OUT)  :: LSTs                   ! Local sidereal time for setting
  REAL(KIND=double) :: GSTr                   ! Greenwich sidereal time for rising
  REAL(KIND=double) :: GSTs                   ! Greenwich sidereal time for setting
  
  !PRINT *, 'delta', delta
  !PRINT *, 'phi', phi
  
  ! Find azimuth of object at rising and setting
  Ar = ACOS( SIN(delta * pi180) / COS(phi * pi180) ) / pi180    ! Result should be in degrees
  As = 360.D0 - Ar
  
  !PRINT *, 'Ar', Ar
  !PRINT *, 'As', As
  
  ! Find the times of rising and setting
  H = ( 1.D0 / 15.D0 ) * ACOS ( -TAN(phi * pi180) * TAN(delta * pi180) ) / pi180 
  
  !PRINT *, 'H', H
  LSTr = 24.D0 + alpha - H
  IF ( LSTr > 24.D0 ) LSTr = LSTr - 24.D0
  
  LSTs = alpha + H
  IF (LSTs > 24.D0 ) LSTs = LSTs - 24.D0
  
  !PRINT *, 'LSTr', LSTr
  !PRINT *, 'LSTs', LSTs
  
  ! Convert these local sidereal times to GST
  CALL LST_to_GST(LSTr, longitude, GSTr)
  CALL LST_to_GST(LSTs, longitude, GSTs)
  
  !PRINT *, 'GSTr', GSTr
  !PRINT *, 'GSTs', GSTs
 
  ! Convert GST to GMT
  CALL GST_to_GMT(GSTr, day, month, year, GMTr)
  CALL GST_to_GMT(GSTs, day, month, year, GMTs)
  
  !PRINT *, 'GMTr', GMTr
  !PRINT *, 'GMTs', GMTs

END SUBROUTINE rough_riseset
  
SUBROUTINE parallax_correction(delta, phi, angular_diameter, horiz_parallax, &
                               atm_refraction, delta_t)
  ! This subroutine calculates the parallax correction, delta_t, to rise and set times.
  ! It follows calculations given in section 32 of Duffett-Smith.
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  REAL(KIND=double), INTENT(IN) :: delta  ! declination (decimal degrees)
  REAL(KIND=double), INTENT(IN) :: phi    ! latitude (decimal degrees)
  REAL(KIND=double), INTENT(IN) :: angular_diameter ! decimal degrees
  REAL(KIND=double), INTENT(IN) :: horiz_parallax   ! decimal degrees
  REAL(KIND=double), INTENT(IN) :: atm_refraction   !decimal degrees
  REAL(KIND=double), INTENT(out) :: delta_t ! time correction (seconds)
  REAL(KIND=double) :: psi    ! Angle between true path and horizon (decimal degrees)
  REAL(KIND=double) :: x      ! 
  REAL(KIND=double), PARAMETER :: pi = 3.141592654D0
  REAL(KIND=double), PARAMETER :: pi180 = (pi / 180.D0) ! Multiply angle variable in degrees by this
  REAL(KIND=double) :: Delta_A   ! Ang. distance between apparent path and true path,
                                            ! parallel to the horizon
  REAL(KIND=double) :: Ar, As               ! Azimuth at rising, setting (decimal degrees)
  REAL(KIND=double) :: A_prime_r, A_prime_s ! Apparent azimuths of rising and setting
  REAL(KIND=double) :: y                    ! see p. 56 of Duffett-Smith (dedimal degrees)
  
  psi = ACOS( SIN(phi * pi180) / COS(delta * pi180)) / pi180
  !PRINT *, 'psi', psi
  
  x = (angular_diameter / 2) + horiz_parallax + atm_refraction
  !PRINT *, 'x', x
  
  Delta_A = ASIN( TAN(x * pi180) / TAN(psi * pi180) ) / pi180
  Ar = ACOS( SIN(delta * pi180) / COS(phi * pi180) ) / pi180    ! Result should be in degrees
  As = 360.D0 - Ar
  A_prime_r = Ar - Delta_A
  
  y = ASIN( SIN(x * pi180) / SIN(psi * pi180) ) / pi180
  !PRINT *, 'y', y
  
  delta_t = (240.D0 * y) / COS(delta * pi180) ! result is in seconds

END SUBROUTINE parallax_correction

SUBROUTINE GMT_to_local(GMT, time_zone, local_time)
  ! This subroutine converts GMT to local time. It follows pseudocode given in
  ! Duffett-Smith, section 10.
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  INTEGER, INTENT(IN) :: time_zone
  REAL(KIND=double), INTENT(IN) :: GMT
  REAL(KIND=double), INTENT(OUT) :: local_time

  local_time = GMT + REAL(time_zone,double)
    
  IF (local_time > 24.D0) local_time = local_time - 24.D0
  IF (local_time < 0.D0) local_time = local_time + 24.D0
  
END SUBROUTINE GMT_to_local

PROGRAM riseset
  IMPLICIT NONE
  INTEGER, PARAMETER :: double = 8
  REAL(KIND=double) :: day
  INTEGER :: day_in
  INTEGER :: month
  INTEGER :: year
  INTEGER :: time_zone ! (Local time) - GMT
!   REAL(KIND=double) :: local_rise_sun, local_set_sun    ! Local rise and set times for sun
!   REAL(KIND=double) :: local_rise_civil_sun, local_set_civil_sun ! Local civil rise and set
  REAL(KIND=double) :: phi! Local latitude
  REAL :: phi_in, phi_min      ! Input latitude, and minutes latitude
  REAL(KIND=double) :: longitude ! Local longitude
  REAL :: longitude_in, longitude_min ! Input local longitude
  REAL(KIND=double) :: julian_date  ! Declare function
  REAL(KIND=double) :: alpha_1  ! The right ascension in hours, decimal form
  REAL(KIND=double) :: delta_1  ! The declination in degrees, decimal form
  REAL(KIND=double) :: lambda_sun
  REAL(KIND=double) :: angular_diameter_sun = 0.533D0 ! degrees
  REAL(KIND=double):: horiz_parallax_sun = (8.79D0 / 60.D0) / 60.D0
  REAL(KIND=double):: atm_refraction_sun = (34.D0 / 60.D0)
  REAL(KIND=double) :: beta_sun = 0.D0
  REAL(KIND=double) :: lambda_2   ! lambda_2 = lambda_sun + 0.985647 degrees (amount sun moves in 24 hours)
  REAL(KIND=double) :: alpha_2, delta_2   ! R.A. and Dec. calculated from lambda_2
  REAL(KIND=double) :: ST1r, ST1s, ST2r, ST2s         ! Rise and set times of sun
  REAL(KIND=double) :: GMTr, GMTs                     ! Rise and set times of sun (GMT)
  REAL(KIND=double), PARAMETER :: pi = 3.141592654D0
  REAL(KIND=double), PARAMETER :: pi180 = (pi / 180.D0) ! Multiply angle variable in degrees by this
  REAL(KIND=double) :: Tr, Ts   ! Actual sidereal time of rising or setting
  REAL(KIND=double) :: delta_prime    ! delta_prime = (delta_1 + delta_2) / 2
  REAL(KIND=double) :: delta_t ! time correction (seconds at first, then hours)
  
  PRINT *, ''
  PRINT *, 'This program calculates the rise and set times of the sun'
  PRINT *, 'for a given date and location.'
  PRINT *, ''
  PRINT *, 'Please enter the date. First, enter the day.'
  READ *, day_in
  PRINT *, ''
  PRINT *, 'Please enter the month as an integer.'
  READ *, month
  PRINT *, 'Please enter the year as an integer.'
  READ *, year
  PRINT *, ''
  PRINT *, 'Please enter the latitude. If you wish to enter it in decimal'
  PRINT *, 'degrees, enter it at this time. If you wish to enter it in'
  PRINT *, 'Degree-minute form, only enter the degree portion at this time.'
  PRINT *, 'Note: Northern latitudes are positive.'
  READ *, phi_in
  PRINT *, ''
  PRINT *, 'Please enter the minutes latitude. If you entered decimal degrees'
  PRINT *, 'at the previous prompt, enter zero.'
  READ *, phi_min
  PRINT *, ''
  PRINT *, 'Please enter the longitude. If you wish to enter it in decimal'
  PRINT *, 'degrees, enter it at this time. If you wish to enter it in'
  PRINT *, 'Degree-minute form, only enter the degree portion at this time.'
  PRINT *, 'Note: Western latitudes are positive.'
  !jc  Should be longitude above as well as below a few lines
  READ *, longitude_in
  PRINT *, ''
  PRINT *, 'Please enter the minutes latitude. If you entered decimal degrees'
  PRINT *, 'at the previous prompt, enter zero.'
  READ *, longitude_min
  PRINT *, ''
  PRINT *, 'Please enter your time zone as an integer. The time zone is the'
  PRINT *, 'number of hours that local time differs from GMT. For example,'
  PRINT *, 'CST is -6.'
  READ *, time_zone
   
  day = REAL(day_in,double)
  phi = REAL(phi_in + (phi_min / 60.),double)
  longitude = REAL(longitude_in + (longitude_min / 60.),double)
  
  PRINT *, ''
  PRINT *, 'You entered:'
  PRINT *
  PRINT *, 'Day:  ', day_in
  PRINT *, 'Month:', month
  PRINT *, 'Year: ', year
  PRINT *, ''
  PRINT *, 'Latitude:', phi_in, 'degrees,', phi_min, 'minutes'
  PRINT *, 'Latitude (decimal form):', phi
  PRINT *, ''
  PRINT *, 'Longitude:', longitude_in, 'degrees,', longitude_min, 'minutes'
  PRINT *, 'Longitude (decimal form):', longitude
  PRINT *, ''
  
  ! For testing only
  !lambda_sun = 163.778867D0
  !alpha_1 = 11.003687D0
  !delta_1 = 6.380389D0
  
  ! Find the ecliptic and equatorial coordinates of the Sun
  CALL sun_coords(day, month, year, lambda_sun, alpha_1, delta_1)
  
  !PRINT *, 'lambda_sun', lambda_sun
  !PRINT *, 'alpha_1', alpha_1
  !PRINT *, 'delta_1', delta_1
  
  ! Find the ecliptic longitude of the sun 24 hours later
  lambda_2 = lambda_sun + 0.985647D0
  
  !PRINT *, 'lambda_2', lambda_2
  
  !PRINT *, 'Test of convert_to_ra_dec'
  
  ! Find the ecliptic and equatorial coordinates of the Sun 24 hours later
  CALL convert_to_ra_dec(lambda_2, beta_sun, alpha_2, delta_2)
  
  !PRINT *, 'alpha_2', alpha_2
  !PRINT *, 'delta_2', delta_2
  
  ! Calculate the local sidereal times of rising and setting appropriate to those coordinates
  CALL rough_riseset(alpha_1, delta_1, phi, longitude, day, month, year, ST1r, ST1s)
  CALL rough_riseset(alpha_2, delta_2, phi, longitude, day, month, year, ST2r, ST2s)
  !PRINT *, 'ST1r', ST1r
  !PRINT *, 'ST1s', ST1s
  !PRINT *, 'ST2r', ST2r
  !PRINT *, 'ST2s', ST2s
  
  ! Calculate Tr and Ts
  Tr = (24.07D0 * ST1r) / (24.07D0 + ST1r - ST2r)
  Ts = (24.07D0 * ST1s) / (24.07D0 + ST1s - ST2s)
  
  !PRINT *, 'Tr', Tr
  !PRINT *, 'Ts', Ts
  
  ! Calculate the corection, delta_t, due to parallax, refraction, and the Sun's finite diameter.
  delta_prime = (delta_1 + delta_2) / 2
  !PRINT *, 'delta_prime', delta_prime
  
  CALL parallax_correction(delta_prime, phi, angular_diameter_sun, horiz_parallax_sun, &
                               atm_refraction_sun, delta_t)
  !PRINT *, 'delta_t', delta_t
  
  ! Convert delta_t from seconds to hours
  delta_t = delta_t / 3600.D0
  
  ! Apply the parallax correction
  Tr = Tr - delta_t
  Ts = Ts + delta_t
  
  ! Tr and Ts are LST times. To convert to local civil times, first convert to GST:
  !PRINT *, 'Tr', Tr
  !PRINT *, 'Ts', Ts
  CALL LST_to_GST(Tr, longitude, Tr)
  CALL LST_to_GST(Ts, longitude, Ts)
  
  ! Tr and Ts are now GST times; convert to GMT times:
  CALL GST_to_GMT(Tr, day, month, year, Tr)
  CALL GST_to_GMT(Ts, day, month, year, Ts)
  !PRINT *, 'Tr', Tr
  !PRINT *, 'Ts', Ts
  
  ! Apply the time zone offset to Tr and Ts
  CALL GMT_to_local(Tr, time_zone, Tr)
  CALL GMT_to_local(Ts, time_zone, Ts)
    
  PRINT *, 'Sunrise:', INT(Tr), 'hours', & 
  REAL((Tr - REAL(INT(Tr),double)) * 60.D0), 'minutes'
  
  PRINT *, 'Sunset: ', INT(Ts), 'hours', & 
  REAL((Ts - REAL(INT(Ts),double)) * 60.D0), 'minutes'
  
END PROGRAM riseset

!jc  Nicely written program - there is a lot here. The style is good,
!jc  as are the commments.