my code:
```
"""
Sunrise and Sunset Calculator

This script calculates the sunrise and sunset times for a given location (latitude and longitude) and date.
The calculation accounts for the date's day of the year and timezone adjustments.
It uses astronomical calculations derived from the NOAA solar calculator and,
relies on custom implementations of trigonometric functions (sin, cos, tan, asin, acos, atan) based on Taylor series expansions.
Most importantly, it uses no external libraries, making it a self-contained solution.

Resources:
- Trigonometric function implementations based on Taylor series:
  - https://en.wikipedia.org/wiki/Taylor_series#Trigonometric_functions
- Sunrise and sunset calculation algorithm:
  - https://en.wikipedia.org/wiki/Sunrise_equation
  - https://www.aa.quae.nl/en/reken/zonpositie.html
  - https://gml.noaa.gov/grad/solcalc/solareqns.PDF
"""

# Define constants
PI = 3.14159265358979323846264338327950288419716939937510  # Pi to high precision
EPSILON = 1e-15  # Threshold for term convergence in series expansions
MAX_TERMS = 10

def radians(degrees):
    """
    Convert degrees to radians.

    Parameters:
        degrees (float): Angle in degrees.

    Returns:
        float: Angle in radians.
    """
    return degrees * PI / 180

def degrees(radians):
    """
    Convert radians to degrees.

    Parameters:
        radians (float): Angle in radians.

    Returns:
        float: Angle in degrees.
    """
    return radians * 180 / PI

def sin(x):
    """
    Calculate the sine of x (in degrees) using Taylor series expansion.

    Parameters:
        x (float): The angle in degrees.

    Returns:
        float: The sine of the angle.
    """
    # Convert degrees to radians
    x = radians(x)
    
    # Range reduction to improve convergence: bring x within (-π, π]
    x = x % (2 * PI)
    if x > PI:
        x -= (2 * PI)
    
    # Initialize variables for Taylor series
    term = x  # First term of the series (x^1 / 1!)
    sin_x = term
    n = 1  # Term counter
    
    while n < MAX_TERMS:
        # Compute the next term using the recurrence relation
        term *= -x * x / ((2 * n) * (2 * n + 1))
        sin_x += term
        
        # Check if the current term is small enough to stop
        if abs(term) < EPSILON:
            break
        n += 1
    
    return sin_x

def cos(x):
    """
    Calculate the cosine of x (in degrees) using Taylor series expansion.

    Parameters:
        x (float): The angle in degrees.

    Returns:
        float: The cosine of the angle.
    """
    # Convert degrees to radians
    x = radians(x)
    
    # Range reduction to improve convergence: bring x within (-π, π]
    x = x % (2 * PI)
    if x > PI:
        x -= (2 * PI)
    
    # Initialize variables for Taylor series
    term = 1.0  # First term of the series (x^0 / 0!)
    cos_x = term
    n = 1  # Term counter
    
    while n < MAX_TERMS:
        # Compute the next term using the recurrence relation
        term *= -x * x / ((2 * n - 1) * (2 * n))
        cos_x += term
        
        # Check if the current term is small enough to stop
        if abs(term) < EPSILON:
            break
        n += 1
    
    return cos_x

def tan(x):
    """
    Calculate the tangent of x (in degrees) using sine and cosine functions.

    Parameters:
        x (float): The angle in degrees.

    Returns:
        float: The tangent of the angle.

    Raises:
        ValueError: If cosine of x is zero (tan is undefined).
    """
    cosine = cos(x)
    if abs(cosine) < EPSILON:
        raise ValueError("Tangent undefined for angles where cosine is zero.")
    return sin(x) / cosine

def asin(x):
    """
    Calculate the arcsine of x using the standard Taylor series expansion.

    Parameters:
        x (float): The value for which to compute arcsin. Must be in [-1, 1].

    Returns:
        float: The arcsine of x in degrees.

    Raises:
        ValueError: If x is outside the domain [-1, 1].
    """
    if x < -1.0 or x > 1.0:
        raise ValueError("asin(x) is undefined for |x| > 1.")

    # Use the standard Taylor series expansion for asin(x)
    term = x
    result = term
    n = 1

    while n < MAX_TERMS:
        term *= x * x * (2 * n - 1) / (2 * n)
        current_term = term / (2 * n + 1)
        result += current_term

        if abs(current_term) < EPSILON:
            break
        n += 1

    return degrees(result)

def acos(x):
    """
    Calculate the arccosine of x using the identity acos(x) = π/2 - asin(x).

    Parameters:
        x (float): The value for which to compute arccos. Must be in [-1, 1].

    Returns:
        float: The arccosine of x in degrees.

    Raises:
        ValueError: If x is outside the domain [-1, 1].
    """
    return 90 - asin(x)

def atan(x):
    """
    Calculate the arctangent of x using Taylor series expansion.

    Parameters:
        x (float): The value for which to compute arctan.

    Returns:
        float: The arctangent of x in degrees.
    """
    # Handle special case where x is 0
    if x == 0.0:
        return 0.0

    # Determine the sign of x
    sign = 1
    if x < 0:
        sign = -1
        x = -x  # Work with positive x for simplicity

    if x <= 1.0:
        # Use Taylor series expansion for |x| <= 1
        term = x  # First term of the series
        result = term
        n = 1  # Term counter

        while n < MAX_TERMS:
            term *= -x * x
            current_term = term / (2 * n + 1)
            result += current_term

            # Check if the current term is small enough to stop
            if abs(current_term) < EPSILON:
                break
            n += 1
    else:
        # Use identity arctan(x) = π/2 - arctan(1/x) for x > 1
        inverse_x = 1.0 / x
        term = inverse_x  # First term of the series
        atan_inverse = term
        n = 1  # Term counter

        while n < MAX_TERMS:
            term *= -inverse_x * inverse_x
            current_term = term / (2 * n + 1)
            atan_inverse += current_term

            # Check if the current term is small enough to stop
            if abs(current_term) < EPSILON:
                break
            n += 1

        result = (PI / 2) - atan_inverse

    # Convert the result from radians to degrees and apply the original sign
    return degrees(sign * result)

def mod(a, b):
    """
    Compute the modulus of a and b.

    Parameters:
        a (float): Dividend.
        b (float): Divisor.

    Returns:
        float: The modulus of a by b.
    """
    return a % b

def is_leap_year(year):
    """
    Determine if a given year is a leap year.

    Parameters:
        year (int): The year to check.

    Returns:
        int: 1 if leap year, 0 otherwise.
    """
    return int(year % 4 == 0 and (year % 100 != 0 or year % 400 == 0))

def day_of_year(year, month, day):
    """
    Calculate the day of the year for a given date.

    Parameters:
        year (int): The year of the date.
        month (int): The month of the date.
        day (int): The day of the date.

    Returns:
        int: The day of the year (1-366).
    """
    # Days in each month, accounting for leap years in February
    month_days = [31, 28 + is_leap_year(year), 31, 30, 31, 30,
                  31, 31, 30, 31, 30, 31]
    
    if day > month_days[month - 1]:
        raise ValueError("Invalid day for the given month.")
    
    return sum(month_days[:month - 1]) + day

def adjust_quadrant(RA, L):
    """
    Adjust the right ascension (RA) to the correct quadrant based on the Sun's true longitude (L).

    Parameters:
        RA (float): Right ascension in degrees.
        L (float): Sun's true longitude in degrees.

    Returns:
        float: Adjusted right ascension in degrees.
    """
    # Ensure both RA and L are within 0-360 degrees
    RA = mod(RA, 360)
    L = mod(L, 360)
    
    # Determine the current quadrants
    L_quadrant  = (int(L / 90)) * 90
    RA_quadrant = (int(RA / 90)) * 90
    
    # Adjust RA to be in the same quadrant as L
    RA = RA + (L_quadrant - RA_quadrant)
    return mod(RA, 360)

def format_time(time_in_hours, seconds=False):
    """
    Format time from decimal hours to HH:MM:SS AM/PM format.

    Parameters:
        time_in_hours (float): Time in decimal hours.

    Returns:
        str: Formatted time string.
    """        
    hour = int(time_in_hours)
    remainder = time_in_hours - hour
    minute = int(remainder * 60)
    remainder = remainder * 60 - minute
    second = int(round(remainder * 60))

    if second >= 60:
        second -= 60
        minute += 1
    if minute >= 60:
        minute -= 60
        hour += 1
    hour %= 24

    # Determine AM or PM
    meridiem = "AM" if hour < 12 else "PM"
    # Convert to 12-hour format
    hour_display = hour % 12
    if hour_display == 0:
        hour_display = 12

    if seconds:
        return f"{hour_display}:{minute:02d}:{second:02d} {meridiem}"
    else:
        if second >= 30:
            minute += 1
            if minute >= 60:
                minute -= 60
                hour += 1
                hour %= 24
                meridiem = "AM" if hour < 12 else "PM"
                hour_display = hour % 12
                if hour_display == 0:
                    hour_display = 12
        return f"{hour_display}:{minute:02d} {meridiem}"

def calculate_sunrise_sunset(latitude, longitude, year, month, day, timezone_offset=0):
    """
    Calculate sunrise and sunset times for a given location and date.

    Parameters:
        latitude (float): Latitude in decimal degrees (positive for north, negative for south).
        longitude (float): Longitude in decimal degrees (positive for east, negative for west).
        year (int): Year of the date.
        month (int): Month of the date.
        day (int): Day of the date.
        timezone_offset (float): Offset from UTC in hours.

    Returns:
        tuple: Sunrise and sunset times in decimal hours (0-24), adjusted for timezone.
               Returns (None, None) if the sun does not rise or set on the given date/location.
    """
    N = day_of_year(year, month, day)  # Day of the year
    
    # Convert longitude to hour value and calculate approximate time
    lng_hour = longitude / 15.0
    
    # Calculate approximate times for sunrise and sunset
    t_rise = N + ((6 - lng_hour) / 24)
    t_set = N + ((18 - lng_hour) / 24)
    
    # Calculate the Sun's mean anomaly
    M_rise = (0.9856 * t_rise) - 3.289
    M_set = (0.9856 * t_set) - 3.289
    
    # Calculate the Sun's true longitude
    # L = M + (1.916 * sin(M)) + (0.020 * sin(2*M)) + 282.634
    L_rise = mod(M_rise + (1.916 * sin(M_rise)) + (0.020 * sin(2 * M_rise)) + 282.634, 360)
    L_set = mod(M_set + (1.916 * sin(M_set)) + (0.020 * sin(2 * M_set)) + 282.634, 360)
    
    # Calculate the Sun's right ascension
    RA_rise = mod(atan(0.91764 * tan(L_rise)), 360)
    RA_set = mod(atan(0.91764 * tan(L_set)), 360)
    
    # Adjust the right ascension to the correct quadrant
    RA_rise = adjust_quadrant(RA_rise, L_rise)
    RA_set = adjust_quadrant(RA_set, L_set)
    
    # Convert RA from degrees to hours
    RA_rise /= 15.0
    RA_set /= 15.0
    
    # Calculate the Sun's declination
    # sin_dec = 0.39782 * sin(L)
    sin_dec_rise = 0.39782 * sin(L_rise)
    cos_dec_rise = cos(asin(sin_dec_rise))
    
    sin_dec_set = 0.39782 * sin(L_set)
    cos_dec_set = cos(asin(sin_dec_set))
    
    zenith = 90.833333  # Official zenith for sunrise/sunset
    
    # Calculate the Sun's local hour angle
    # cos_h = (cos(zenith) - sin_dec * sin(latitude)) / (cos_dec * cos(latitude))
    cos_h_rise = (cos(zenith) - (sin_dec_rise * sin(latitude))) / (cos_dec_rise * cos(latitude))
    cos_h_set = (cos(zenith) - (sin_dec_set * sin(latitude))) / (cos_dec_set * cos(latitude))
    
    # Check for polar day/night
    if cos_h_rise > 1 or cos_h_set > 1:
        # The sun never rises or sets on this location and date (Arctic/Antarctic)
        return None, None
    
    # Calculate the local hour angle in degrees
    H_rise = (360 - acos(cos_h_rise)) / 15.0  # Sunrise hour angle
    H_set = acos(cos_h_set) / 15.0            # Sunset hour angle
    
    # Calculate local mean time of rising and setting
    # T = H + RA - (0.06571 * t) - 6.622
    T_rise = H_rise + RA_rise - (0.06571 * t_rise) - 6.622
    T_set = H_set + RA_set - (0.06571 * t_set) - 6.622
    
    # Convert to UTC
    UT_rise = mod(T_rise - lng_hour, 24)
    UT_set = mod(T_set - lng_hour, 24)
    
    # Adjust for timezone (no Daylight Saving Time considerations as of now)
    sunrise_time = mod(UT_rise + timezone_offset, 24)
    sunset_time = mod(UT_set + timezone_offset, 24)
    
    return sunrise_time, sunset_time

if __name__ == "__main__":
    # Location coordinates
    latitude  = 13.0326372542501    # Latitude in degrees (positive for north)
    longitude = 77.5755860908572    # Longitude in degrees (positive for east)
    
    # Date and timezone information
    day, month, year = 31, 12, 2024
    timezone_offset = 0  # UTC offset in hours (e.g., IST is UTC+5.5)
    
    # Calculate sunrise and sunset times
    sunrise, sunset = calculate_sunrise_sunset(latitude, longitude, year, month, day, timezone_offset)
    
    if sunrise is not None and sunset is not None:
        # Format and display the results
        print(f"Sunrise: {format_time(sunrise, seconds=True)}")
        print(f"Sunset: {format_time(sunset, seconds=True)}")
    else:
        # Handle cases where sunrise or sunset does not occur
        print("Sunrise and sunset times not available for this location on the given date.")
```

Astral
```
from astral import LocationInfo
from astral.sun import sun
from datetime import datetime, timedelta
import pytz

def round_time(dt):
    if dt.second >= 30:
        return dt.replace(second=0, microsecond=0) + timedelta(minutes=1)
    return dt.replace(second=0, microsecond=0)

def calculate_sunrise_sunset_astral(latitude, longitude, timezone, date=None):
    if date is None:
        date = datetime.now()
    timezone = pytz.timezone(timezone)
    location = LocationInfo(latitude=latitude, longitude=longitude, timezone=timezone.zone)
    s = sun(location.observer, date=date, tzinfo=timezone)
    # sunrise = round_time(s['sunrise'])
    # sunset = round_time(s['sunset'])
    sunrise = s['sunrise']
    sunset = s['sunset']
    return sunrise, sunset

if __name__ == "__main__":
    # Example coordinates (e.g., New York City)
    # Location coordinates
    latitude  = 13.0326372542501    # Latitude in degrees (positive for north)
    longitude = 77.5755860908572    # Longitude in degrees (positive for east)
    
    # Date and timezone information
    day, month, year = 31, 12, 2024
    
    date = datetime(year, month, day)
    timezone = 'UTC'
    sunrise, sunset = calculate_sunrise_sunset_astral(latitude, longitude, timezone, date)
    # sunrise_str = sunrise.strftime("%dth %B %Y %I:%M %p")
    # sunset_str = sunset.strftime("%dth %B %Y %I:%M %p")
    sunrise_str = sunrise.strftime("%dth %B %Y %I:%M:%S %p")
    sunset_str = sunset.strftime("%dth %B %Y %I:%M:%S %p")
    print(f"Sunrise: {sunrise_str}")
    print(f"Sunset: {sunset_str}")
```

testcode
```
from solar_calculator import calculate_sunrise_sunset, format_time
from testbed import calculate_sunrise_sunset_astral
from datetime import datetime

if __name__ == "__main__":
    # Location coordinates
    latitude  = 13.0326372542501    # Latitude in degrees (positive for north)
    longitude = 77.5755860908572    # Longitude in degrees (positive for east)
    
    # Date and timezone information
    _, _, year = 31, 12, 2024
    timezone_offset = 0  # UTC offset in hours (e.g., IST is UTC+5.5)
    timezone = 'UTC'
    
    for month in range(3, 5):
        print(f"{'-' * 50}")
        print(f"Month: {month}")
        for day in range(1, 31):
            print(f'Day: {day}')
            # Calculate sunrise and sunset times
            sunrise, sunset = calculate_sunrise_sunset(latitude, longitude, year, month, day, timezone_offset)
            
            if sunrise is not None and sunset is not None:
                # Format and display the results
                print(f"Sunrise Custom: {format_time(sunrise, seconds=True)}")
                print(f"Sunset Custom: {format_time(sunset, seconds=True)}")
            else:
                # Handle cases where sunrise or sunset does not occur
                print("Sunrise and sunset times not available for this location on the given date.")
                
            date = datetime(year, month, day)
            sunrise_astral, sunset_astral = calculate_sunrise_sunset_astral(latitude, longitude, timezone, date)
            sunrise_astral_str = sunrise_astral.strftime("%dth %B %Y %I:%M:%S %p")
            sunset_astral_str = sunset_astral.strftime("%dth %B %Y %I:%M:%S %p")
            print(f"Sunrise Astral: {sunrise_astral_str}")
            print(f"Sunset Astral: {sunset_astral_str}")
            print(f"{'-' * 50}")
```

Results from test
```
PS C:\InfinityX\Automa\Sunset> & C:/Users/themy/AppData/Local/Programs/Python/Python312/python.exe c:/InfinityX/Automa/Sunset/test.py
--------------------------------------------------
Month: 3
Day: 1
Sunrise Custom: 1:05:30 AM
Sunset Custom: 12:58:42 PM
Sunrise Astral: 01th March 2024 01:05:41 AM
Sunset Astral: 01th March 2024 12:58:24 PM
--------------------------------------------------
Day: 2
Sunrise Custom: 1:04:56 AM
Sunset Custom: 12:58:51 PM
Sunrise Astral: 02th March 2024 01:05:07 AM
Sunset Astral: 02th March 2024 12:58:33 PM
--------------------------------------------------
Day: 3
Sunrise Custom: 1:04:22 AM
Sunset Custom: 12:59:00 PM
Sunrise Astral: 03th March 2024 01:04:33 AM
Sunset Astral: 03th March 2024 12:58:41 PM
--------------------------------------------------
Day: 4
Sunrise Custom: 1:03:48 AM
Sunset Custom: 12:59:08 PM
Sunrise Astral: 04th March 2024 01:03:59 AM
Sunset Astral: 04th March 2024 12:58:50 PM
--------------------------------------------------
Day: 5
Sunrise Custom: 1:03:13 AM
Sunset Custom: 12:59:16 PM
Sunrise Astral: 05th March 2024 01:03:24 AM
Sunset Astral: 05th March 2024 12:58:57 PM
--------------------------------------------------
Day: 6
Sunrise Custom: 1:02:37 AM
Sunset Custom: 12:59:24 PM
Sunrise Astral: 06th March 2024 01:02:48 AM
Sunset Astral: 06th March 2024 12:59:05 PM
--------------------------------------------------
Day: 7
Sunrise Custom: 1:02:01 AM
Sunset Custom: 12:59:31 PM
Sunrise Astral: 07th March 2024 01:02:12 AM
Sunset Astral: 07th March 2024 12:59:12 PM
--------------------------------------------------
Day: 8
Sunrise Custom: 1:01:25 AM
Sunset Custom: 12:59:38 PM
Sunrise Astral: 08th March 2024 01:01:36 AM
Sunset Astral: 08th March 2024 12:59:19 PM
--------------------------------------------------
Day: 9
Sunrise Custom: 1:00:48 AM
Sunset Custom: 12:59:44 PM
Sunrise Astral: 09th March 2024 01:00:59 AM
Sunset Astral: 09th March 2024 12:59:26 PM
--------------------------------------------------
Day: 10
Sunrise Custom: 1:00:11 AM
Sunset Custom: 12:59:50 PM
Sunrise Astral: 10th March 2024 01:00:22 AM
Sunset Astral: 10th March 2024 12:59:32 PM
--------------------------------------------------
Day: 11
Sunrise Custom: 12:59:33 AM
Sunset Custom: 12:59:56 PM
Sunrise Astral: 11th March 2024 12:59:44 AM
Sunset Astral: 11th March 2024 12:59:38 PM
--------------------------------------------------
Day: 12
Sunrise Custom: 12:58:55 AM
Sunset Custom: 1:00:02 PM
Sunrise Astral: 12th March 2024 12:59:07 AM
Sunset Astral: 12th March 2024 12:59:43 PM
--------------------------------------------------
Day: 13
Sunrise Custom: 12:58:17 AM
Sunset Custom: 1:00:07 PM
Sunrise Astral: 13th March 2024 12:58:28 AM
Sunset Astral: 13th March 2024 12:59:49 PM
--------------------------------------------------
Day: 14
Sunrise Custom: 12:57:39 AM
Sunset Custom: 1:00:13 PM
Sunrise Astral: 14th March 2024 12:57:50 AM
Sunset Astral: 14th March 2024 12:59:54 PM
--------------------------------------------------
Day: 15
Sunrise Custom: 12:57:00 AM
Sunset Custom: 1:00:18 PM
Sunrise Astral: 15th March 2024 12:57:11 AM
Sunset Astral: 15th March 2024 12:59:59 PM
--------------------------------------------------
Day: 16
Sunrise Custom: 12:56:21 AM
Sunset Custom: 1:00:22 PM
Sunrise Astral: 16th March 2024 12:56:32 AM
Sunset Astral: 16th March 2024 01:00:04 PM
--------------------------------------------------
Day: 17
Sunrise Custom: 12:55:42 AM
Sunset Custom: 1:00:27 PM
Sunrise Astral: 17th March 2024 12:55:53 AM
Sunset Astral: 17th March 2024 01:00:09 PM
--------------------------------------------------
Day: 18
Sunrise Custom: 12:55:02 AM
Sunset Custom: 1:00:31 PM
Sunrise Astral: 18th March 2024 12:55:14 AM
Sunset Astral: 18th March 2024 01:00:13 PM
--------------------------------------------------
Day: 19
Sunrise Custom: 12:54:23 AM
Sunset Custom: 1:00:36 PM
Sunrise Astral: 19th March 2024 12:54:34 AM
Sunset Astral: 19th March 2024 01:00:17 PM
--------------------------------------------------
Day: 20
Sunrise Custom: 12:53:43 AM
Sunset Custom: 1:00:40 PM
Sunrise Astral: 20th March 2024 12:53:55 AM
Sunset Astral: 20th March 2024 01:00:22 PM
--------------------------------------------------
Day: 21
Sunrise Custom: 12:53:03 AM
Sunset Custom: 1:00:44 PM
Sunrise Astral: 21th March 2024 12:53:15 AM
Sunset Astral: 21th March 2024 01:00:26 PM
--------------------------------------------------
Day: 22
Sunrise Custom: 12:52:23 AM
Sunset Custom: 1:00:48 PM
Sunrise Astral: 22th March 2024 12:52:35 AM
Sunset Astral: 22th March 2024 01:00:30 PM
--------------------------------------------------
Day: 23
Sunrise Custom: 12:51:43 AM
Sunset Custom: 1:00:52 PM
Sunrise Astral: 23th March 2024 12:51:55 AM
Sunset Astral: 23th March 2024 01:00:34 PM
--------------------------------------------------
Day: 24
Sunrise Custom: 12:51:03 AM
Sunset Custom: 1:00:55 PM
Sunrise Astral: 24th March 2024 12:51:15 AM
Sunset Astral: 24th March 2024 01:00:37 PM
--------------------------------------------------
Day: 25
Sunrise Custom: 12:50:23 AM
Sunset Custom: 1:00:59 PM
Sunrise Astral: 25th March 2024 12:50:35 AM
Sunset Astral: 25th March 2024 01:00:41 PM
--------------------------------------------------
Day: 26
Sunrise Custom: 12:49:43 AM
Sunset Custom: 1:01:03 PM
Sunrise Astral: 26th March 2024 12:49:55 AM
Sunset Astral: 26th March 2024 01:00:45 PM
--------------------------------------------------
Day: 27
Sunrise Custom: 12:49:03 AM
Sunset Custom: 1:01:06 PM
Sunrise Astral: 27th March 2024 12:49:15 AM
Sunset Astral: 27th March 2024 01:00:49 PM
--------------------------------------------------
Day: 28
Sunrise Custom: 12:48:23 AM
Sunset Custom: 1:01:10 PM
Sunrise Astral: 28th March 2024 12:48:35 AM
Sunset Astral: 28th March 2024 01:00:52 PM
--------------------------------------------------
Day: 29
Sunrise Custom: 12:47:43 AM
Sunset Custom: 1:01:14 PM
Sunrise Astral: 29th March 2024 12:47:55 AM
Sunset Astral: 29th March 2024 01:00:56 PM
--------------------------------------------------
Day: 30
Sunrise Custom: 12:47:03 AM
Sunset Custom: 1:01:17 PM
Sunrise Astral: 30th March 2024 12:47:16 AM
Sunset Astral: 30th March 2024 01:01:00 PM
--------------------------------------------------
--------------------------------------------------
Month: 4
Day: 1
Sunrise Custom: 12:45:44 AM
Sunset Custom: 1:01:25 PM
Sunrise Astral: 01th April 2024 12:45:57 AM
Sunset Astral: 01th April 2024 01:01:08 PM
--------------------------------------------------
Day: 2
Sunrise Custom: 12:45:05 AM
Sunset Custom: 1:01:29 PM
Sunrise Astral: 02th April 2024 12:45:17 AM
Sunset Astral: 02th April 2024 01:01:12 PM
--------------------------------------------------
Day: 3
Sunrise Custom: 12:44:25 AM
Sunset Custom: 1:01:33 PM
Sunrise Astral: 03th April 2024 12:44:38 AM
Sunset Astral: 03th April 2024 01:01:16 PM
--------------------------------------------------
Day: 4
Sunrise Custom: 12:43:46 AM
Sunset Custom: 1:01:37 PM
Sunrise Astral: 04th April 2024 12:43:59 AM
Sunset Astral: 04th April 2024 01:01:20 PM
--------------------------------------------------
Day: 5
Sunrise Custom: 12:43:08 AM
Sunset Custom: 1:01:41 PM
Sunrise Astral: 05th April 2024 12:43:20 AM
Sunset Astral: 05th April 2024 01:01:25 PM
--------------------------------------------------
Day: 6
Sunrise Custom: 12:42:29 AM
Sunset Custom: 1:01:46 PM
Sunrise Astral: 06th April 2024 12:42:42 AM
Sunset Astral: 06th April 2024 01:01:29 PM
--------------------------------------------------
Day: 7
Sunrise Custom: 12:41:51 AM
Sunset Custom: 1:01:50 PM
Sunrise Astral: 07th April 2024 12:42:04 AM
Sunset Astral: 07th April 2024 01:01:34 PM
--------------------------------------------------
Day: 8
Sunrise Custom: 12:41:13 AM
Sunset Custom: 1:01:55 PM
Sunrise Astral: 08th April 2024 12:41:26 AM
Sunset Astral: 08th April 2024 01:01:38 PM
--------------------------------------------------
Day: 9
Sunrise Custom: 12:40:35 AM
Sunset Custom: 1:02:00 PM
Sunrise Astral: 09th April 2024 12:40:48 AM
Sunset Astral: 09th April 2024 01:01:43 PM
--------------------------------------------------
Day: 10
Sunrise Custom: 12:39:58 AM
Sunset Custom: 1:02:05 PM
Sunrise Astral: 10th April 2024 12:40:11 AM
Sunset Astral: 10th April 2024 01:01:49 PM
--------------------------------------------------
Day: 11
Sunrise Custom: 12:39:21 AM
Sunset Custom: 1:02:10 PM
Sunrise Astral: 11th April 2024 12:39:34 AM
Sunset Astral: 11th April 2024 01:01:54 PM
--------------------------------------------------
Day: 12
Sunrise Custom: 12:38:44 AM
Sunset Custom: 1:02:16 PM
Sunrise Astral: 12th April 2024 12:38:58 AM
Sunset Astral: 12th April 2024 01:02:00 PM
--------------------------------------------------
Day: 13
Sunrise Custom: 12:38:08 AM
Sunset Custom: 1:02:21 PM
Sunrise Astral: 13th April 2024 12:38:22 AM
Sunset Astral: 13th April 2024 01:02:05 PM
--------------------------------------------------
Day: 14
Sunrise Custom: 12:37:32 AM
Sunset Custom: 1:02:27 PM
Sunrise Astral: 14th April 2024 12:37:46 AM
Sunset Astral: 14th April 2024 01:02:11 PM
--------------------------------------------------
Day: 15
Sunrise Custom: 12:36:57 AM
Sunset Custom: 1:02:34 PM
Sunrise Astral: 15th April 2024 12:37:11 AM
Sunset Astral: 15th April 2024 01:02:18 PM
--------------------------------------------------
Day: 16
Sunrise Custom: 12:36:22 AM
Sunset Custom: 1:02:40 PM
Sunrise Astral: 16th April 2024 12:36:36 AM
Sunset Astral: 16th April 2024 01:02:24 PM
--------------------------------------------------
Day: 17
Sunrise Custom: 12:35:48 AM
Sunset Custom: 1:02:47 PM
Sunrise Astral: 17th April 2024 12:36:02 AM
Sunset Astral: 17th April 2024 01:02:31 PM
--------------------------------------------------
Day: 18
Sunrise Custom: 12:35:14 AM
Sunset Custom: 1:02:54 PM
Sunrise Astral: 18th April 2024 12:35:28 AM
Sunset Astral: 18th April 2024 01:02:38 PM
--------------------------------------------------
Day: 19
Sunrise Custom: 12:34:41 AM
Sunset Custom: 1:03:01 PM
Sunrise Astral: 19th April 2024 12:34:55 AM
Sunset Astral: 19th April 2024 01:02:46 PM
--------------------------------------------------
Day: 20
Sunrise Custom: 12:34:08 AM
Sunset Custom: 1:03:09 PM
Sunrise Astral: 20th April 2024 12:34:22 AM
Sunset Astral: 20th April 2024 01:02:53 PM
--------------------------------------------------
Day: 21
Sunrise Custom: 12:33:35 AM
Sunset Custom: 1:03:17 PM
Sunrise Astral: 21th April 2024 12:33:50 AM
Sunset Astral: 21th April 2024 01:03:01 PM
--------------------------------------------------
Day: 22
Sunrise Custom: 12:33:04 AM
Sunset Custom: 1:03:25 PM
Sunrise Astral: 22th April 2024 12:33:18 AM
Sunset Astral: 22th April 2024 01:03:10 PM
--------------------------------------------------
Day: 23
Sunrise Custom: 12:32:33 AM
Sunset Custom: 1:03:33 PM
Sunrise Astral: 23th April 2024 12:32:47 AM
Sunset Astral: 23th April 2024 01:03:18 PM
--------------------------------------------------
Day: 24
Sunrise Custom: 12:32:02 AM
Sunset Custom: 1:03:42 PM
Sunrise Astral: 24th April 2024 12:32:17 AM
Sunset Astral: 24th April 2024 01:03:27 PM
--------------------------------------------------
Day: 25
Sunrise Custom: 12:31:32 AM
Sunset Custom: 1:03:51 PM
Sunrise Astral: 25th April 2024 12:31:47 AM
Sunset Astral: 25th April 2024 01:03:36 PM
--------------------------------------------------
Day: 26
Sunrise Custom: 12:31:03 AM
Sunset Custom: 1:04:01 PM
Sunrise Astral: 26th April 2024 12:31:18 AM
Sunset Astral: 26th April 2024 01:03:46 PM
--------------------------------------------------
Day: 27
Sunrise Custom: 12:30:34 AM
Sunset Custom: 1:04:10 PM
Sunrise Astral: 27th April 2024 12:30:49 AM
Sunset Astral: 27th April 2024 01:03:56 PM
--------------------------------------------------
Day: 28
Sunrise Custom: 12:30:06 AM
Sunset Custom: 1:04:20 PM
Sunrise Astral: 28th April 2024 12:30:21 AM
Sunset Astral: 28th April 2024 01:04:06 PM
--------------------------------------------------
Day: 29
Sunrise Custom: 12:29:39 AM
Sunset Custom: 1:04:30 PM
Sunrise Astral: 29th April 2024 12:29:54 AM
Sunset Astral: 29th April 2024 01:04:17 PM
--------------------------------------------------
Day: 30
Sunrise Custom: 12:29:11 AM
Sunset Custom: 1:04:39 PM
Sunrise Astral: 30th April 2024 12:29:28 AM
Sunset Astral: 30th April 2024 01:04:28 PM
--------------------------------------------------
PS C:\InfinityX\Automa\Sunset> 
```

German Astronomy
```
Astronomy Answers
Position of the Sun
[Astronomy Answers] [Dictionary] [AnswerBook] [Universe Family Tree] [Science] [Starry Sky] [Planet Positions] [Calculate] [Colophon]
1. Time ... 2. The Mean Anomaly ... 3. The Equation of Center ... 4. The Perihelion and the Obliquity of the Ecliptic ... 5. The Ecliptical Coordinates ... 6. The Equatorial coordinates ... 7. The Observer ... 8. Solar Transit ... 9. Equation of Time ... 10. Sunrise and Sunset ... 11. The Duration of Sunrise and Sunset ... 12. Earth ... 13. The Netherlands ... 14. Local dependence on location ... 15. Comparison with Planetarium Programs ... 16. Derivation of (Π) and (ε) ... 17. Derivation of ( θ_0 ) ... 18. Planet-based Topographic Coordinates ... 19. Planet-Bound Horizontal Coordinates ... 20. Coordinate Transformations Between All Systems ... 21. North Pole, South Pole, Primary Pole


The position of the Sun in the sky as seen from a planet (such as the Earth) is determined by four things:

the time.
the motion of the planet in its orbit around the Sun, which does not happen at constant speed, because of the eccentricity of the orbit of the planet.
the angle between the axis of rotation of the planet and the plane of the orbit of the planet, which is not equal to 90 degrees. This causes the seasons.
the location of the observer on the planet, which determines how high the Sun can get in the sky.
Below, I provide formulas that take all of these things into account. I neglected small effects here and there, to keep the formulas simple. Remember that you can add multiples of 360 degrees to any angle without changing its direction.

I updated the tables on this page in October 2016 so that they are based on the relevant planet characteristics published by the IAU in 2011, rather than those of 2000. The most important change is that Pluto is no longer regarded as a major planet but now as a dwarf planet, which means that now (since 2006) the other pole of Pluto is regarded as its north pole than before. See Chapter 21.

When in this page an ecliptic or pole or coordinates such as ecliptic longitude or declination are mentioned, then these are appropriate for the planet where the observer is. Such coordinates are based on the orbit and equator of that planet and are not identical to the coordinates of the same names that we use on Earth. For example, the ecliptic of Mars, the ecliptic that is appropriate for Mars, is not the same as the ecliptic of Earth, and an earth-based atlas of the stars is of no use if you have a right ascension and declination calculated for the Sun as seen from Mars.

For examples, we'll calculate the position of the Sun for 1 April 2004 at 12:00 UTC, as seen from 52° north latitude and 5° east longitude (Netherlands) on Earth, and as seen from 14°36' south latitude and 184°36' west longitude (Gusev crater) on Mars.

1. Time
For this kind of astronomical calculations, it is convenient to express the date and time using an unending day numbering scheme. Such a scheme is provided by the Julian Day Number 
. The Julian Day Number Calculation Page explains how you can calculate the Julian Day Number for a date in the Gregorian calendar. For the calculations of the position of the Sun you should express time as Universal Time (UTC), and this includes the Julian Date. For example, JD 2453144.5 corresponds to 0:00 hours UTC on 19 May 2004, which is (for example) equivalent to 1:00 hours Central European Time on 19 May 2004, or 20:00 hours (8 pm) Eastern Standard Time on 18 May 2004.

Only on Earth do the seasons repeat themselves after about one calendar year (in the western ― Gregorian ― calendar), so, only for the Earth, the day number 
 in the calendar year can be used instead of the Julian Day Number: 
 = 1 corresponds to 0:00 UTC on January 1st, 2 to 0:00 UTC on January 2nd, 32 to 0:00 UTC on February 1st, and so on.

2. The Mean Anomaly
Because we see the Sun from the planet, we see the motion of the planet around the Sun reflected in the apparent motion of the Sun along the ecliptic, relative to the stars.

If the orbit of the planet were a perfect circle, then the planet as seen from the Sun would move along its orbit at a fixed speed, and then it would be simple to calculate its position (and also the position of the Sun as seen from the planet). The position that the planet would have relative to its perihelion if the orbit of the planet were a circle is called the mean anomaly, indicated in the formulas as 
.

You can calculate the mean anomaly of the planets (measured in degrees) to reasonable accuracy using the following formula, for a date given as a Julian Day Number (JD) 
:

 

You should take 
 (in degrees) and 
 (in degrees per day) from the following table:

Mercury	174.7948	4.09233445
Venus	50.4161	1.60213034
Earth	357.5291	0.98560028
Mars	19.3730	0.52402068
Jupiter	20.0202	0.08308529
Saturn	317.0207	0.03344414
Uranus	141.0498	0.01172834
Neptune	256.2250	0.00598103
Pluto	14.882	0.00396
For the Earth, you can also use the following formula:

 

where 
 is the time since 00:00 UTC at the beginning of the most recent January 1st, measured in (whole and fractional) days.

The chosen date and time correspond to Julian Day Number 2453097, so 
, and also to 
. For the Earth, using equation 
 we find 
 and with equation 
 we find 
. Because equation 
 is a bit more accurate that equation 
, we'll use equation 
. For Mars, we find 
.

3. The Equation of Center
The orbits of the planets are not perfect circles but rather ellipses, so the speed of the planet in its orbit varies, and therefore the apparent speed of the Sun along the ecliptic also varies throughout the planet's year.

The true anomaly (symbol ν
, nu) is the angular distance of the planet from the perihelion of the planet, as seen from the Sun. For a circular orbit, the mean anomaly and the true anomaly are the same. The difference between the true anomaly and the mean anomaly is called the Equation of Center, written here as 
:

ν
 

To calculate the Equation of Center or the true anomaly from the mean anomaly, you must first solve the Equation of Kepler. This equation cannot be solved in general, but you can find approximations to the solution that are more and more accurate. See the Calculation Page for the Equation of Kepler for more information about this. If the orbit looks more like a circle than like a parabola, then the following approximation is sufficiently accurate:

 

You can find the coefficients 
 through 
 in the following table. They depend on the eccentricity 
 of the planet's orbit. If a particular coefficient is not mentioned, then it is equal to zero (to four decimal places). The orbits of Mercury and Pluto deviate the most from circularity, so they need the most coefficients. The column 
 shows the maximum error that you get if you use the approximation with the coefficients from the table.

Mercury	23.4400	2.9818	0.5255	0.1058	0.0241	0.0055	0.0026
Venus	0.7758	0.0033					0.0000
Earth	1.9148	0.0200	0.0003				0.0000
Mars	10.6912	0.6228	0.0503	0.0046	0.0005		0.0001
Jupiter	5.5549	0.1683	0.0071	0.0003			0.0001
Saturn	6.3585	0.2204	0.0106	0.0006			0.0001
Uranus	5.3042	0.1534	0.0062	0.0003			0.0001
Neptune	1.0302	0.0058					0.0001
Pluto	28.3150	4.3408	0.9214	0.2235	0.0627	0.0174	0.0096
For the Earth, we find

 

and thence ν
. For Mars we find

 

and ν
.

4. The Perihelion and the Obliquity of the Ecliptic
To find the position of the Sun in the sky we need to know what the ecliptic longitude 
 is of the perihelion of the planet, relative to the ecliptic and vernal equinox (the ascending equinox) of the planet. The ecliptic of the planet is the plane of the orbit of the planet, which makes an angle with the orbit (ecliptic) of the Earth (and that angle is called the inclination of the orbit). The vernal equinox of the planet is the point where the Sun passes from south to north through the plane of the equator of the planet. We also need to know the obliquity ε
 of the equator of the planet compared to the orbit of the planet. These two values are listed in the following table for each planet, measured in degrees.


ε
Mercury	230.3265	0.0351
Venus	73.7576	2.6376
Earth	102.9373	23.4393
Mars	71.0041	25.1918
Jupiter	237.1015	3.1189
Saturn	99.4587	26.7285
Uranus	5.4634	82.2298
Neptune	182.2100	27.8477
Pluto	184.5484	119.6075
These values are different than those given on the NASA Fact Sheet (//nssdc.gsfc.nasa.gov/planetary/planetfact.html), because there they are measured relative to the ecliptic of the Earth (the plane of the orbit of the Earth), while here they are measured relative to the ecliptic of the planet (the plane of the orbit of the planet).

If you use the orbital elements of a planet that are measured relative to the ecliptic of the Earth, then the coordinates that you calculate from those elements are relative to the Earth's ecliptic, too. You can then meaningfully compare those coordinates to coordinates of stars and other objects in star atlases and catalogs made for use on Earth, to see where the planet is compared to those stars and other objects, but you cannot use them to calculate where the Sun is relative to the horizon on that planet.

On the other hand, if you use orbital elements that are measured relative to the ecliptic of the planet (as on this page), then you cannot meaningfully compare the resulting coordinates with star atlases and catalogs that were made for use on Earth, but you can use them to calculate where the Sun is relative to the horizon of that planet.

5. The Ecliptical Coordinates
The ecliptical longitude λ
 (lambda) is the position along the ecliptic, relative to the vernal equinox (so relative to the stars). The mean longitude 
 is the ecliptical longitude that the planet would have if the orbit were a perfect circle. That is

 

The ecliptic longitude of the planet as seen from the Sun is equal to

λν
 

If you look at the Sun from the planet, then you're looking in exactly the opposite direction than if you look at the planet from the Sun, so those directions are 180° apart. So, the ecliptic longitude of the Sun, as seen from the planet, is equal to

λ
λν
 

The value of λ
 determines when the (astronomical) seasons begin: when λ
, then spring begins in the northern hemisphere, and autumn in the southern hemisphere. Each next multiple of 90° brings the start of the next season.

The ecliptic latitude β
 (beta) of the Sun, the perpendicular distance of the Sun from the ecliptic, is always so small that we can ignore it here. With this, we now have the ecliptic coordinates of the Sun.

As seen from Earth we find λ
 and as seen from Mars λ
.

6. The Equatorial coordinates
The equatorial coordinate system in the sky is tied to the rotation axis of the planet. The equatorial coordinates are the right ascension α
 (alpha) and the declination δ
 (delta). The declination determines from which parts of the planet the object can be visible, and the right ascension determines (together with other things) when the object is visible.

With these formulas you can calculate the equatorial coordinates from the ecliptic coordinates:

αδ	λεββε
αδ	λβ
δ	βεβελ
 

For the Sun we have β
, so then

α
λ
ελ
δ
λ
ε
 

For future reference we define

α
λ
 

If ε
 is close enough to 0° or 180°, and if we neglect small terms, then we can approximate the relationship between the right ascension α
 and the ecliptic longitude λ
 of the Sun as seen from the planet by

λε	λ
 
ε
 
ε
 
ε
λ
 
ε
 
ε
λ
 
ε
λε
α
λ
λ
λ
λ
λ
 

and the relationship between the declination δ
 and the ecliptic longitude by

λε	ε
 
ε
 
ε
λ
 
ε
 
ε
λ
 
ε
λε
δ
λ
λ
λ
 

with 
, 
, 
, 
, 
, and 
 (measured in degrees) from the following table. The columns 
 and 
 show the greatest errors you make for α
 and δ
 when you use these approximations.

Mercury	−0.0000			0.0000	0.0351			0.0000
Venus	−0.0304			0.0001	2.6367	0.0009		0.0036
Earth	−2.4657	0.0529	−0.0014	0.0003	22.7908	0.5991	0.0492	0.0003
Mars	−2.8608	0.0713	−0.0022	0.0004	24.3880	0.7332	0.0706	0.0011
Jupiter	−0.0425			0.0001	3.1173	0.0015		0.0034
Saturn	−3.2338	0.0909	−0.0031	0.0009	25.7696	0.8640	0.0949	0.0010
Uranus	−42.5874	12.8117	−2.6077	17.6902	56.9083	−0.8433	26.1648	3.34
Neptune	−3.5214	0.1078	−0.0039	0.0163	26.7643	0.9669	0.1166	0.060
Pluto	−19.3248	3.0286	−0.4092	0.5052	49.8309	4.9707	5.5910	0.19
The approximations for Uranus aren't very good, because Uranus lies almost on its side: You'll do better to use the full formulas for Uranus.

As seen from Earth, and using equation 
, we find:

α

and with equation 
 we find

δ

As seen from Mars, we find

α

and

δ

Using equations 
 and 
 we find for the Earth

α
δ
 

and for Mars

α
δ
 

so the approximations yield practically the same results in these cases as the full equations 
 and 
.

7. The Observer
Where a celestial body is in your sky depends on your geographical coordinates (latitude φ
 [phi] north, longitude 
 west), on the position of the body between the stars (its equatorial coordinates α
 and δ
), and on the rotation angle of the planet at your location, relative to the stars. That latter angle is expressed in the sidereal time θ
 (theta). The sidereal time is the right ascension that is on the celestial meridian at that moment. The sidereal time is equal to

θθ
θ
 

with θ
 and θ
 from the next table. The 
 means that on both sides of the 
 you can add or subtract any multiple of 360°. For θ
 the custom is to choose the value between 0° and 360°.


θ
θ
Mercury	132.3282	6.1385025
Venus	104.9067	−1.4813688
Earth	280.1470	360.9856235
Mars	313.3827	350.89198226
Jupiter	145.9722	870.5360000
Saturn	174.3508	810.7939024
Uranus	29.6474	−501.1600928
Neptune	52.4160	536.3128662
Pluto	122.2370	56.3625225
θ
 is the sidereal time (in degrees) at longitude 0° at the instant defined by 
. θ
 is the rate of change of the sidereal time, in degrees per day.

For the Netherlands on Earth we find

θ

and for Gusev on Mars

θ

The position of a celestial body in the sky is specified by its altitude 
 above the horizon, and its azimuth 
. The altitude is 0° at the horizon, +90° in the zenith (straight over your head), and −90° in the nadir (straight down). The azimuth is the direction along the horizon, which we measure from south to west. South has azimuth 0°, west +90°, north +180°, and east +270° (or −90°, that's the same thing). The altitude and azimuth are the horizontal coordinates. To calculate the horizontal coordinates from the equatorial coordinates, you can use the following formulas:

δ
φδδφ
φδφδ
θα
φδφ
 

The 
 is the hour angle, which indicates how long ago (measured in sidereal time) the celestial body passed through the celestial meridian.

For the Netherlands on Earth we find

 

For Gusev on Mars we find

 

At 12:00 UTC on 1 April 2004, the Sun as seen from the Netherlands stands about 5° west of south at 43° above the horizon, and as seen from the Gusev crater on Mars the Sun then stands about 3° south by northwest at 61° above the horizon.

8. Solar Transit
The transit of a celestial body is the moment at which the body passes through the celestial meridian. The transit of the Sun is noon, the middle of the day, at 12 hours solar time. The hour angle 
 of the Sun is then equal to 0. We have

θα
 

Using this and previous equations, you can find the time of transit by making a guess for a value of 
 for which equation 
 holds, calculating θ
 and α
 for it, and then checking whether they satisfy equation 
. If they do not, then you have to adjust 
. All in all, it is a big search operation for the correct value of 
. Moreover, equation 
 does not help you to understand which things are most important to the time of solar transit. You do get such insight if you approximate the solution by neglecting smaller terms. The answer may then not be as accurate as the correct solution, but does indicate clearly what the solution looks like, is usually easier to calculate, and provides an excellent starting guess in the search for the real value of 
.

With equations 
, 
, 
, 
, and 
, and after omitting smaller terms, we find

 
 

where

θ
 
 
 
θ
 
 

If you already have λ
, then you can use it instead of 
 in equation 
: that is a tiny bit more accurate. 
 provides the date and time of a transit of the Sun. 
 shows by how much the time of transit can vary because of the eccentricity 
 of the orbit. 
 indicates how much the time of transit can change because of the obliquity ε
 of the ecliptic. 
 is the average length of the solar day (from one transit to the next). All values are measured in Earth days of 24 hours.

Mercury	45.3497	11.4556		175.9386
Venus	52.1268	−0.2516	0.0099	−116.7505
Earth	0.0009	0.0053	−0.0068	1.0000000
Mars	0.9047	0.0305	−0.0082	1.027491
Jupiter	0.3345	0.0064		0.4135778
Saturn	0.0766	0.0078	−0.0040	0.4440276
Uranus	0.1260	−0.0106	0.0850	−0.7183165
Neptune	0.3841	0.0019	−0.0066	0.6712575
Pluto	4.5635	−0.5024	0.3429	6.387672
To find the date and time of a solar transit near Julian Date 
, you now proceed as follows:


Calculate

 
 
 

and then take for 
 the whole number nearest to 
.

Then calculate

 

This 
 is a reasonable estimate for the date and time of the transit near 
, except that the 
 and 
 corrections are not in it yet.

Calculate 
 and 
 for 
 and then get a better estimate for the date and time of the solar transit from

 

If you want greater precision, then calculate the hour angle 
 of the Sun for 
 and take


 
 

as improved value of 
. You can repeat this until 
 no longer changes.


 is a Julian Date, which counts days and which is equal to a whole number at 12:00 UTC (i.e., noon, UTC). So, a Julian Date like 2453096.9898 that ends in .9898 is .0102 days before the next whole number, i.e., 0.0102 days = 0.0102×24 = 0.245 hours = 0.0102×24×60 = about 15 minutes before 12:00 UTC, i.e., about 11:45 UTC.
For our example, we looked near 
. Which solar transit is closest to that in the Netherlands and in Gusev crater on Mars? For the Netherlands (
) we find

 
 

so 
, so


For 
 we find 
 and then


With that, we find

 

The repetition method provides increased accuracy. For JD 2453096.9766 we find 
, so a better estimate is 
 
. For that moment we find 
 equal to 0 to four digits after the decimal point, which is sufficiently precise. The solar transit at 5° east longitude happens on 1 April 2004 around 11:45 UTC.

For the Gusev crater on Mars (
) we find

 
 

so 
, so


For 
 we find 
 and


With that, we find

 

Now we use the repetition method. For JD 2453096.8945 we find 
, so a better estimate is 
 
. That yields no further change to 
 (to four digits after the decimal point) so we are done. The solar transit in Gusev crater happens on 1 April 2004 around 10:32 UTC.


9. Equation of Time
Above, we calculated when the Sun is highest in the sky, expressed in the Earthly time scale of the Julian Date. With this, you still don't know at what time that is on your clock, and that time scale is definitely not very handy if you are on another planet where the day has a different length from the day on Earth that the Julian Date is based on.

If you are not interested in transforming to times and dates for other planets, then you just want to know at what time the Sun is highest in the sky on your planet according to your clock, which is tied to your solar time and the season on your planet. Solar time is the time determined by the Sun. Three different kinds of solar time are relevant:

The true solar time is a time scale in which the Sun is highest in the sky at exactly 12:00:00 hours each day, but not every day needs be equally long (as measured using something that does not care about the position of the Sun in the sky, such as an atomic clock or planetary phenomena in the sky). You can read true solar time from a well-adjusted sundial.
The mean solar time is a time scale that runs in step with true solar time on average but in which all days are equally long. Averaged over many planet days and planet years, the Sun is highest in the sky at 12:00:00 mean solar time. You can read mean solar time from a well-adjusted clock or watch. This time is often referred to as "Local Time" in these pages.
The official clock time that you use in daily life is the mean solar time of a certain meridian that belongs to your time zone. The difference between official clock time and the mean solar time of a location depends only on the difference between the geographical longitude of that place and the longitude of the meridian that the official clock time is based on.
The answer to the question at what time the Sun is highest in the sky depends a lot on the time scale that you use. In true solar time, the answer is "always at exactly 12 o'clock", but true solar time is otherwise not very convenient at all. You can't use it if the Sun does not shine or if it is hidden behind clouds, and a mechanical clock or watch would have to be a lot more complicated if it had to show true solar time. In mean solar time, the answer is "on average at 12 o'clock", but it may differ from 12 o'clock from day to day. In official clock time, the answer is "on average at a fixed time", but that fixed time depends on your location, though it is usually near 12 o'clock.

So how large is the difference between the mean solar time and the true solar time? This is called the Equation of Time. The Equation of Time is how much you have to add to true solar time ("sundial time") to get mean solar time ("local time").

We found earlier (equation 
) for the moment at which the Sun is highest in the sky (transit) that θα
. For the right ascension α
 of the Sun we found α
. For the sidereal time θ
 measured in degrees we construct a different formula:

θ
 

Here 
 is the mean solar time measured in degrees (0° = midnight, 180° = noon). The explanation for this equation is as follows: The difference between the sidereal time and the solar time reflects the motion of the planet around the Sun so it increases by 360° in a planet year. The sidereal time is tied to the regular rotation of the planet around its axis, which does not notice the varying speed at which the planet orbits around the Sun (associated with eccentricity 
) or the position of the Sun above or below the celestial equatior (associated with the obliquity ε
 of the ecliptic). Therefore, the sidereal time is tied to the mean longitude 
 of the Sun, which increases at a fixed rate by 360° during a planet year, and which is independent of 
 and ε
. During the descending equinox (when λ
) the sidereal time is equal to the mean solar time, so we need the extra 180° in the formula. During a planet day, θ
 increases by 360° (the planet turns once around its axis) plus or minus a small amount because the sidereal time runs a little faster or slower than the mean solar time, but that difference is already caught in 
.

If we now set θ
 equal to α
 then we find

θα
 

or

 

The Equation of Time 
 is equal to the difference between 
 and 180° (because 180° corresponds to 12:00:00 local mean solar time):

 

If you divide the Equation of Time measured in degrees by 15 then you get the Equation of Time measured in hours.

 depends on eccentricity 
 of the planet's orbit and on the mean anomaly 
 of the planet, i.e., on where the planet is compared to its perihelion. 
 depends on the obliquity ε
 of the ecliptic of the planet and on the longitude λ
 of the Sun, i.e., on the season. A first-order approximation is:

λ
 

where we've omitted many smaller terms. The following table shows the amplitudes, i.e., the largest contribution that the 
 and 
 terms give to the Equation of Time for each planet, measured in planet minutes.

Mercury	94.5	0
Venus	3.1	0.1
Earth	7.7	9.9
Mars	42.8	11.4
Jupiter	22.2	0.2
Saturn	25.4	13.0
Uranus	21.2	178.1
Neptune	4.1	14.1
Pluto	114.6	69.3
For the Earth, the contributions of the orbit (
) and the season (
) are about equally large; for Uranus and Neptune the influence of the season is much greater than that of the orbit; and for the other planets the influence of the orbit is much greater than the influence of the seasons.

10. Sunrise and Sunset
For the hour angle that corresponds to 
 we find

δφ
 

This equation has two solutions: if 
 is a solution, then 
 is a solution, too. The solution with 
 is associated with sunset, and the solution with 
 is associated with sunrise. We write 
 for the solution that has 
. By using formulas found earlier and by omitting small terms, we find

λ
φ
λ
φ
φ
λ
φ
φ
φ
 

with

ε
 
ε
 
ε
 
ε
 
ε
 
ε
 

where you must measure ε
 in radians instead of in degrees. You transform from degrees to radians by multiplying the number of degrees by π
.

The values of 
, 
, and 
 are listed for all planets in the following table, measured in degrees. For Uranus and Pluto it is doubtful that the approximations yield reasonable results.

Mercury	0.035		
Venus	2.636	0.001	
Earth	22.137	0.599	0.016
Mars	23.576	0.733	0.024
Jupiter	3.116	0.002	
Saturn	24.800	0.864	0.032
Uranus	28.680	−0.843	8.722
Neptune	26.668	0.967	0.039
Pluto	38.648	4.971	1.864
Sunrise is the moment at which the top of the solar disk touches the horizon in the morning. Sunset is the same, but at night. To calculate the times of sunrise and sunset, you must take into account (besides the effects mentioned earlier) also the following things:

The Sun does not look like a point but like a disk, so when the center of the Sun sets (
), then half of the Sun is still above the horizon. The size of the Sun must be taken into account.
The atmosphere of the planet (if it has one) bends light downward so that a celestial body appears slightly higher in the sky than it would have done without the atmosphere. This refraction of the light is strongest near the horizon, where the light has had to travel furthest through the air. On Earth, at the horizon as seen from sea level, the refraction amounts to some 0.57° on average. This means that a celestial body that with the atmosphere just set would have already been 0.57° below the horizon without the atmosphere.
To compensate for these two effects, you can set 
 equal to the 
 from the following table (measured in degrees), instead of 0°. The refraction has been added in to the value for the Earth, but not for the other planets, because they have either no atmosphere of any consequence, or else they have one that is so dense that you cannot see the Sun at all from the surface. The table also gives the average diameter 
 of the solar disk (in degrees).

Mercury	−0.69	1.38	−0.0120
Venus	−0.37	0.74	−0.0064
Earth	−0.83	0.53	−0.0146
Mars	−0.17	0.35	−0.0031
Jupiter	−0.05	0.10	−0.0009
Saturn	−0.03	0.06	−0.0005
Uranus	−0.01	0.03	−0.0002
Neptune	−0.01	0.02	−0.0002
Pluto	−0.01	0.01	−0.0001
With equation 
 you can then find the hour angle at which the Sun with that α
δ
 rises or sets:

φδ
φδ
 
 

Here again, 
 holds for sunset and 
 for sunrise. The estimated JD of sunrise and sunset are then

 
 
 

You can improve these estimates by repetition similar to how we did it for the transit. For 
 or 
 calculate the hour angle 
 and declination δ
 of the Sun, and then use equation 
 to calculate the hour angle 
 that the Sun at those equatorial coordinates should have to be rising or setting. Then replace the estimates by

 
 
 

until they don't change anymore, to the desired precision.

We found earlier that the transit of the Sun as seen from the Earth occurs at 
. The declination of the Sun is then δ
. According to the above table, we have for the Earth 
. With φ
 we then find from equation 
 that

 

The first estimate for the preceding sunrise is then

 

For that instant of time we calculate the position of the Sun and find 
 and δ
, and hence 
. Then the JD correction is 
 
, so a better estimate is 
.

If we repeat this improvement procedure once more then the correction is equal to 0 (to 4 digits after the decimal point), so we're done.

Sunrise occurs around 05:15 UTC.

The first estimate for the following sunset is

 

For that instant of time we calculate the position of the Sun and find 
 and δ
, and hence 
. Then the JD correction is 
 
, so a better estimate is 
.

If we repeat this improvement procedure once more then the correction is equal to 0 (to 4 digits after the decimal point), so we're done.

Sunset occurs on 1 April 2004 around 18:15 UTC.

We found earlier that the transit of the Sun as seen from Gusev crater on Mars occurs at 
. The declination of the Sun is then δ
. According to the above table, we have for Mars 
. With φ
 we then find from equation 
 that

 

The first estimate for the preceding sunrise is then

 

For that instant of time we calculate the position of the Sun and find 
 and δ
, and hence 
. Then the JD correction is 
 
, so that first estimate was already accurate enough.

Sunrise in Gusev crater happens on 1 April 2004 around 04:27 UTC.

The first estimate for the following sunset is

 

For that instant of time we calculate the position of the Sun and find 
 and δ
, and hence 
. Then the JD correction is 
 
, so a better estimate is 
. Another repeat does not improve this result anymore.

Sunset in Gusev crater happens on 1 April 2004 around 16:37 UTC.

If you watch sunrise or sunset from great altitude, then there are two more corrections to be made:

The horizon appears to be lower, by an amount in degrees that, on Earth, is approximately equal to the square root of the altitude in kilometers. For an altitude of 10 km (6 mi), the horizon is about 3.2° lower.
There is more refraction of sunlight, because the sunlight almost touches the ground but then goes back up through the atmosphere to get to the observer. The total refraction is at most twice as great as it is at sea level, and I expect that the refraction from jet cruising altitude (10 km) above the Earth is almost at that maximum, so it is about 0.5° more than at sea level.
The correction 
 to equation 
 that is necessary because of the large Sun and the atmosphere is approximately equal to

φ
δ
 
 

For the Netherlands and Belgium (with φ
 about 50°) this is about 5 to 6 minutes.

11. The Duration of Sunrise and Sunset
To calculate the duration of sunrise or sunset you must calculate the moment that the top of the solar disk touches the horizon, as described above (with 
), and also the moment that the bottom of the solar disk touches the horizon. To calculate that last one, you should use 
 instead of 
.

12. Earth
For places on Earth, we can combine various approximations into the following (with results measured in hours UTC):

 
 
 
 
φ
φ
φ
 
 
φ
φ
φ
 
 
 
φ
φ
φ
 
 
φ
φ
φ
 
 

13. The Netherlands
For a spot at φ
 and 
 (in the middle of the Netherlands) we find, approximately, measured in Central European Time (CET, equal to UTC plus one hour):

 

These formulas, though approximations, are reasonably accurate. A comparison, for all days of the year 2000, between the results of the full formulas and the results of the last approximate formulas shows differences of at most 3.5 minutes for the times of sunrise and sunset, and less than 1 minute for the time of transit.

14. Local dependence on location
From the preceding formulas we can deduce how the times of sunrise, transit, and sunset depend on the location (near any chosen spot). The times of sunrise, transit, and sunset (measured in a fixed time zone) get earlier by 4 minutes for each degree that you go to the east. This corresponds to 
φ
 
 seconds per kilometer. The time of transit (i.e., high noon) doed not depend on the geographical latitude. In addition, sunrise gets earlier and sunset later by a number of minutes per degree that you go north that is approximately equal to 
φ
 
 and that is equivalent to about 
φ
 
 seconds per kilometer. The shift of the times of sunrise and sunset with latitude is greatest at the beginning of summer and winter. At the beginning of spring and autumn it is 0, so then the times of sunrise and sunset do not depend on the latitude (just like the time of transit).

15. Comparison with Planetarium Programs
The planetarium programs Redshift 3, Redshift 5, Stellarium 0.15.0, the NASA web page JPL HORIZONS, and the current procedure (AA) yield the following values for the azimuth and altitude of the Sun in the sky at JD 2451545.0 (1-1-2000, 12:00 UTC) as seen from the location with latitude and longitude equal to zero on each of the planets:

Redshift 5	Redshift 3	AA	Stellarium 0.15.0	JPL HORIZONS
Mercury	89 59 36	−4 28 58	90 00	−4 21	270 01 32	−4 29 32	89 59 42.8	−4 22 03.8	90.0256	−4.4830
Venus	263 39 25	−70 00 13	263 39	−70 00	83 39 27	−69 59 26	92 14 33.2	−14 59 18.5	263.6545	−70.0006
Earth	178 03 33	+66 57 05	177 37	+66 57	357 18 25	+66 56 23	178 02 03.4	+66 57 37.5	178.0722	+66.9528
Mars	233 17 03	+44 46 23	233 21	+44 41	53 03 40	+45 04 15	233 20 09.0	+44 43 21.0	233.2109	+44.8716
Jupiter	275 36 38	−56 47 58	273 19	+22 27	93 19 50	+23 06 54	283 53 42.5	+77 07 12.8	273.3132	+22.3831
Saturn	115 08 60	+33 21 18	115 08	+33 18	294 58 03	+32 59 53	115 08 41.7	+33 21 41.8	115.1490	+33.3541
Uranus	223 07 12	+45 26 37	223 09	+45 25	44 13 02	+45 46 42	213 51 11.4	−51 47 54.9	223.1205	+45.4433
Neptune	217 28 20	−54 09 28	218 03	−54 32	37 49 37	−54 44 36	217 55 18.7	−53 54 11.7	217.4714	−54.1581
Pluto	125 32 00	−43 59 37	125 32	−44 00	125 35 22	−42 06 27	119 07 54.6	+30 48 30.8	305.5559	−42.1332
And here are similar results for JD 2453097.0 (1-4-2004, 12:00 UTC):

Redshift 5	Redshift 3	AA	Stellarium 0.15.0	JPL HORIZONS
Mercury	89 52 29	−87 19 22	89 01	−87 06	269 19 40	−87 19 34	89 56 12.4	−87 12 34.2	89.3290	−87.3182
Venus	266 46 35	+35 02 25	266 47	+35 03	86 46 41	+35 02 37	264 45 51.5	−59 43 08.3	266.7781	+35.0387
Earth	11 08 57	+85 07 29	11 10	+85 07	194 28 06	+85 05 14	11 13 57.3	+85 07 14.8	11.1353	+85.1259
Mars	77 36 54	−63 14 20	77 46	−62 54	257 30 37	−63 34 42	77 38 47.5	−63 07 28.1	77.5625	−63.3588
Jupiter	92 10 30	+46 10 15	91 36	+20 10	271 35 20	+19 04 53	96 19 20.0	+76 11 48.7	91.5977	+19.6703
Saturn	231 05 23	+47 32 41	231 24	+47 11	50 42 48	+47 56 56	231 16 56.8	+47 20 24.3	231.0880	+47.5457
Uranus	143 35 09	−72 11 32	144 19	−72 22	321 08 50	−72 43 10	109 18 45.9	+41 28 47.3	143.5871	−72.1924
Neptune	173 45 46	−61 31 02	172 54	−61 57	352 58 27	−61 58 26	174 30 03.9	−61 33 40.9	173.7614	−61.5171
Pluto	135 36 03	−41 02 35	135 37	−41 04	135 41 23	−39 00 02	127 19 52.2	+27 17 46.0	315.6817	−39.0309
If you take into account that I measure azimuth from the south but the other sources measure azimuth from the north (corresponding to a difference of 180°), then the results fit to within about a degree, except that the position from Jupiter as calculated by Redshift 5 is very different from that calculated by the other sources, and that the azimuth of the Sun from Pluto calculated by JPL HORIZONS differs by 180° from that of the other sources. I therefore suspect that Redshift 5 is in error for Jupiter. My azimuth for Pluto differs about 180° from that of JPL because I count azimuth from the south and JPL counts it from the north. The other sources are older than, or presumably haven't been updated for, the change in which of Pluto's poles is regarded as the north pole.

16. Derivation of 
 and ε
Table 1 mentions 
 en ε
. How do you calculate those?

The orientation of the orbit of a planet in space is expressed by three independent angles relative to a chosen fundamental plane and a chosen primary direction in that fundamental plane. These angles are part of the orbital elements. We use the inclination 
, the longitude 
 of the ascending node, and the argument ω
 of the perihelion as the three angles. The inclination is the angle between the orbit and the fundamental plane. The ascending node is one of the two intersections between the orbit and the fundamental plane. The longitude of the ascending node is the angle between the primary direction and the ascending node, as seen from the relevant focus of the orbit. Likewise, the argument of the perihelion is the angle between the ascending node and the perihelion of the orbit.

The calculation of the position of the Sun in the sky of a planet is easiest if the fundamental plane is chosen to be equal to the plane of the orbit of the planet, and the primary direction is chosen to be the intersection of the fundamental plane and the extension of the planet's equator into the sky. In the remainder of this text, we assume that that choice has been made. In the case of the Earth, the fundamental plane is the plane of the ecliptic, and the primary direction is the vernal equinox, and the coordinates are with respect to "the equator and equinox of the date".

Then the inclination of the orbit is zero, and the nodes aren't defined, and the ecliptic latitude of the Sun is zero. In that case the only angle needed to specify the orientation of the orbit is the angle between the primary direction and the perifocus of the orbit, measured in the plane of the orbit. That is the planet-based ecliptic longitude of the perifocus. We use symbol 
 for that angle.

Once the planet-based ecliptic coordinates of the Sun are known, then we convert them into planet-based equatorial coordinates. For this we need to know the angle ε
 between the fundamental plane (chosen to be the plane of the orbit of the planet) and the equator of the planet. For the Earth, that angle is called the obliquity of the ecliptic.

So, to calculate the Sun's equatorial coordinates relative to the planet's equator, we need to know 
 and ε
 of that planet.

How do we calculate ε
? ε
 is the angle between the plane of the planet's orbit and the plane of the planet's equator. This is also the angle between the primary pole of the planet (perpendicular to the planet's equator) and the primary pole of the planet's orbit (perpendicular to the planet's orbit). The planet and its orbit each have two poles. One of those two has been chosen as the primary pole. Usually it is the north pole. See Chapter 21.

We calculate the direction of the planet's primary pole from its Earth-based equatorial coordinates, which are available on the internet. And we calculate the direction of the planet orbit's primary pole from its Earth-based orbital elements.

How do we calculate 
? 
 is the angle between the planet's primary direction and the perifocus of the planet's orbit. We can calculate the direction of the perifocus of the planet in Earth-based coordinates from the widely available Earth-based orbital elements of the planet, so we need to calculate the Earth-based coordinates of the planet's primary direction. That primary direction is one of the intersections of the plane of the planet's orbit (the fundamental plane) and the plane of the extension of the planet's equator, so it is perpendicular to the primary pole of the planet and to the primary pole of the orbit of the planet.

Here are the mathematical details. We assume that the following quantities are given:

α
right ascension of north pole
δ
declination of north pole
longitude of ascending node
inclination of the orbit
ω
argument of perifocus
rotation angle at epoch
rotation speed at epoch
ε
obliquity of the ecliptic of Earth
All of these orbital elements except for 
 are assumed specified with respect to the ecliptic (= orbit) or equator of the Earth, not those of the planet.

We use the following notation:



 is an orthogonal matrix of 3 by 3 elements that represents the conversion of cartesian coordinates from coordinate system 
 to coordinate system 
.


 for any 
 is a column vector with 3 elements that represents the cartesian coordinates of a direction in space, with the coordinates specified in coordinate system 
. The three elements of the vector are indicated as 
. The length of all vectors named below is equal to 1.


 are column vectors that represent three mutually perpendicular coordinate axes of coordinate system 
, with the coordinates specified in coordinate system 
.


ε
 is the angle over which to rotate the equator of coordinate system 
 to get coordinate system 
. The axis around which to rotate depends on the context.

The following directions are of interest. They are defined in terms of lines and planes in space, which should be imagined extended out to infinity. Parallel lines in space then correspond to two opposite points on the celestial sphere, of which one is called the primary one. Parallel planes in space then correspond to a great circle on the celestial sphere. All lines perpendicular to those planes correspond to two opposite points on the celestial sphere, which are called the poles of the plane, of which one is called the primary pole of the planes. Intersecting planes correspond to two great circles on the celestial sphere, which intersect in two points, of which one is called the primary intersection,



The primary pole of the object. We define the primary pole to be the one from above which the object seems to rotate around its axis in the counterclockwise direction. For most objects that is what the IAU calls the north pole, but for Venus and Uranus it is what the IAU calls the south pole. See Sec. 21.


The primary pole of the orbit of the object. The two poles of the orbit are the opposite directions perpendicular to the plane of the orbit of the object. The below formulas and values of 
 and 
 define which of the two poles is the primary one.


The primary direction of the orbit of the planet. This is one of the two intersections of the plane of the equator of the planet with the plane of the orbit of the planet. The below formulas define which of the two intersections is the primary one. For the Earth, the vernal equinox is the primary direction.


The primary location of the planet. This is the location with planetographic longitude and latitude equal to zero, which is the intersection of the prime meridian of the planet (longitude zero) and the equator of the planet (latitude zero).


The primary intersection of the equator of the planet with the equator of the Earth. The below formulas define which of the two intersections is the primary one. For the Earth, the primary intersection is not defined.

The coordinate systems that are associated with the Earth have a lowercase letter for symbol, and the coordinate systems that are associated with the other planet have the corresponding uppercase letter for symbol. The coordinate systems are:



 is the Earth-based equatorial coordinate system. It uses the north pole of the Earth as its z axis, and the Earth-based primary direction (the vernal equinox) as its x axis.


 is likewise the planet-based equatorial coordinate system, with the primary pole of the planet as its z axis, and the planet-based primary direction as its x axis.


 is the Earth-based ecliptic coordinate system. It uses the primary pole of the Earth orbit as its z axis, and the Earth-based primary direction (the vernal equinox) as its x axis.


Likewise, 
 is the planet-based ecliptic coordinate system, with the primary pole of the planet orbit as its z axis, and the planet-based primary direction as its x axis.


 is the planet-based hybrid equatorial coordinate system, which uses the primary pole of the planet as its z axis, and uses the primary intersection of the equator of the Earth and the equator of the planet as its x axis. 
 is needed because the IAU uses it to define the rotation angle 
 of the planet.


 is the planet-based topographic coordinate system, which uses the primary pole of the planet as its z axis, and uses the primary location on the planet as its x axis. Coordinate system 
 rotates around its z axis at a rate of once per sidereal planet day, relative to coordinate systems 
 and 
.

We derive formulas for converting between all of those coordinate systems and the Earth-based ecliptic coordinate system, which means we can convert between all of them. All of the described coordinate systems are right-handed.

The relevant quantities to begin with for all major planets and Pluto are listed in the following table. All angles are given in degrees. 
 is given in degrees per Earth day (of 86400 seconds). α
δ
 come from Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009 and Erratum to: Reports of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2006 & 2009, found at //astrogeology.usgs.gov/groups/iau-wgccre. ω
 come from [Meeus] and are based on the VSOP87 model of Betagnon en Francou, except for Pluto. All values are given for the J2000.0 epoch (JD 2451545).

α
δ
ω
Mercury	281.0097	61.4143	48.330893	7.004986	29.125226	329.5469	6.1385025
Venus	272.76	67.16	76.679920	3.394662	54.883787	160.20	−1.4813688
Earth	0	90	174.873174	0	288.064174	190.147	360.9856235
Mars	317.68143	52.88650	49.558093	1.849726	286.502141	176.630	350.89198226
Jupiter	268.056595	64.495303	100.464441	1.303270	273.866868	284.95	870.5360000
Saturn	40.589	83.537	113.665524	2.488878	339.391263	38.90	810.7939024
Uranus	257.311	−15.175	74.005947	0.773196	98.999212	203.81	−501.1600928
Neptune	299.36	43.46	131.784057	1.769952	276.339634	253.18	536.3128492
Pluto	132.993	−6.163	110.307	17.140	113.768	302.695	56.3625225
For the rotation angle 
 of the planet at the time given by Julian Day Number 
 we have

 

As an example, we'll calculate εθ
 for Mars. I show all coordinates rounded to 7 digits after the decimal marker, and all angles to 4 digits after the decimal marker. For the obliquity of the ecliptic ε
 I assume the value 23.4392911°.


Fig. 1: Coordinate systems C, c, Q, q
Fig. 1: Coordinate systems C, c, Q, qFigure 1 shows the Earth-based equatorial coordinate system 
 in fuchsia, the Earth-based ecliptic coordinate system 
 in red, the planet-based ecliptic coordinate system 
 in yellow, and the planet-based equatorial coordinate system 
 in blue.

The black circle represents the celestial sphere of infinite size, to which the Sun and stars and planets appear to be attached.


All displayed ellipses are actually projections of circles on the celestial sphere, seen at an angle. Each circle represents the intersection of the celestial sphere and a flat plane of interest. The part of each circle on the near side of the celestial sphere is displayed as a solid curve, and the part on the far side of the celestial sphere is displayed as a dotted curve.


All displayed straight lines from the center of the celestial sphere are projections of lines from the center of the sphere to a point on the surface of the sphere. Each of those lines represents the infinite extension of an axis or pole of interest.


Actually, the circles and straight lines from the center represent entire families of parallel planes and lines in space. Because the celestial sphere is infinitely large, parallel lines in space seem to meet in a single point on the celestial sphere, just like parallel railroad tracks seem to meet at the horizon. Likewise, parallel planes in space seem to meet in a single great circle on the celestial sphere. In the current discussion we're not interested in differences in the origins (the zero points in space) of the coordinate systems, but only in their orientations in space. For the purposes of this discussion, we move the coordinate systems so that their origins coincide, without rotating them.


The fuchsia curve represents the plane of the equator of the Earth. 
 represents the primary pole of the Earth, which is the north pole. 
 is the primary direction of the orbit of the Earth, one of the two intersections of the equator of the Earth and the orbit of the Earth. In the case of the Earth, it is called the vernal equinox. 
 is perpendicular to both 
 and 
 and completes the Earth-based equatorial coordinate system.


The red curve represents the plane of the orbit of the Earth, i.e., the ecliptic. 
 represents the primary pole of the orbit of the Earth, i.e., the north pole of the ecliptic. 
 completes the Earth-based ecliptic coordinate system. ε
 is the angle over which to rotate the equatorial coordinate system 
 around its x axis to get the ecliptic coordinate system 
. In the case of the Earth, it is called the obliquity of the ecliptic.


The yellow curve represents the plane of the orbit of the planet. 
 represents the primary pole of the orbit of the planet. 
 is the primary direction of the orbit of the planet, one of the two intersections of the equator of the planet and the orbit of the planet. 
 completes the planet-based ecliptic coordinate system.


The blue curve represents the plane of the equator of the planet. 
 represents the primary pole of the planet. 
 is the primary direction, and 
 completes the planet-based equatorial coordinate system.


ε
ε
 is the angle over which to rotate the Earth-based equatorial coordinate system 
 around its x axis to get the Earth-based ecliptic coordinate system 
. It is called the obliquity of the ecliptic.


εε
 is the angle over which to rotate the planet-based equatorial coordinate system 
 around its x axis to get the planet-based ecliptic coordinate system 
. One might call it the planet-based obliquity of the ecliptic.


Point 
 is the ascending node of the orbit of the planet (relative to the Earth's orbit). Angle 
 is the angle from the Earth-based primary direction (the vernal equinox) to the ascending node.


ε
 is the inclination of the orbit of the planet (relative to the Earth's orbit). It is the angle over which to rotate the plane of the Earth's orbit around the ascending node to get the plane of the planet's orbit.


Point 
 is the perifocus of the orbit of the planet: the point at which it is closest to the Sun.


ω
 is the argument of the perifocus of the planet. It is the angle over which to rotate the ascending node around the z axis of the orbit of the planet to get 
.


 is the planet-based ecliptic longitude of the perifocus of the planet. It is the angle over which to rotate the planet-based primary direction around the z axis of the orbit of the planet to get 
.




 are the coordinates of the primary pole of the planet, expressed in Earth-based equatorial coordinates.

α
δ
α
δ
δ
 
 


For Mars we find
 



Converting from Earth-based equatorial coordinates to Earth-based ecliptic coordinates goes via a rotation around the x axis by angle ε
:

ε
ε
ε
ε
ε
 
 

and converting in the opposite direction, from Earth-based ecliptic coordinates to Earth-based equatorial coordinates goes via the opposite rotation:


ε
ε
ε
ε
ε
 
 

 are the coordinates of the primary pole of the planet, expressed in Earth-based ecliptic coordinates.


α
δ
ε
δ
ε
α
δ
ε
δ
ε
α
δ
 
 


For Mars we find
 

which corresponds to Earth-based ecliptic longitude 352.9076° and latitude +63.2821°.



 are the coordinates of the primary pole of the planet orbit, expressed in Earth-based ecliptic coordinates.

 
 
 

where 
 represents rotation around the z axis, similar to 
 described before:


 
 


For Mars we find
 

which corresponds to Earth-based ecliptic longitude 319.5581° and latitude +88.1503°.



ε
 is the obliquity of the planet's orbit, which is the angle between the primary pole of the planet and the primary pole of the orbit of the planet. It has a value between 0 and 180°. Its cosine is equal to the inner product of the two poles:

ε	
ε
δ
ε
α
δ
α
δ
ε
δ
ε
α
δ
 
 
 


but going via the cosine leads to loss of precision for small values of ε
 (see the page about Distances in the Sky). For this, a better alternative is to go via


 
ε
 
 


If after that you are interested in ε
 or ε
, then they can be calculated from


ε	
 
ε	
 
 


For Mars we find
ε		
 
 
ε		
 
 
 
 
 
 
 
 
 
 
ε		
ε		
 
 
ε		
 
 
 


 are the coordinates of the primary direction of the planet, expressed in Earth-based ecliptic coordinates. The primary direction is perpendicular to the primary pole of the planet and also to the primary pole of the orbit of the planet. We calculate this using the vector cross product of the two poles. The cross product is perpendicular to both factors. The order of the two factors determines which of the two perpendicular directions is chosen.

 
ε
 
 


The division by the length of the cross product or ε
 is necessary because the length of this vector cross product isn't necessarily equal to 1.



For Mars we find
 
 
 
 
 
 

which corresponds to Earth-based ecliptic longitude 265.0421° and latitude −1.0740°.


Now that we have 
 and 
, we can complete the coordinate basis 
 for converting from Earth-based ecliptic coordinates to planet-based ecliptic coordinates. We combine the three vectors as rows of the matrix.

 
 


where the 
 symbol represents the transposition operator which turns rows of a matrix into columns and columns into rows.


The coordinate basis 
 for converting from planet-based ecliptic coordinates to Earth-based ecliptic coordinates is the inverse of 
, which inverse is equal to


 


For Mars we find
 
 
 


which corresponds to Earth-based ecliptic longitude 355.0703° and latitude −1.5058°. Then


 
 
 



 are the coordinates of the perifocus of the orbit of the planet, expressed in Earth-based ecliptic coordinates.

ω
 
ωω
ωω
ω
 
 


For Mars we find
 
 
 
 
 

which corresponds to Earth-based ecliptic longitude 336.0684° and latitude −1.7735°.



 are the coordinates of the perifocus of the orbit of the planet, expressed in planet-based ecliptic coordinates. We find it by converting 
 from Earth-based ecliptic coordinates to planet-based ecliptic coordinates.

 

The perifocus is part of the orbit of the planet, so it lies in the planet-based plane of the ecliptic, so its planet-based ecliptic z coordinate and latitude are equal to 0.



For Mars we find
 
 
 

which corresponds to Earth-based ecliptic longitude 71.0041° and latitude 0°.



 is the planet-based ecliptic longitude of the perifocus.

 


For Mars we find


We can now finish the planet-based equatorial coordinate basis. We already had 
, and the primary direction is shared between the planet-based equatorial coordinate system and the planet-based ecliptic coordinate system, so 
 and we already had 
. We complete the coordinate basis 
.

 
 


For Mars we find
 
 
 
 
 
 

17. Derivation of θ

Fig. 2: Coordinate systems C, E, Q, q, T
Fig. 2: Coordinate systems C, E, Q, q, TFigure 2 shows the Earth-based equatorial coordinate system 
 in fuchsia, the planet-based equatorial coordinate system 
 in blue, the planet-based ecliptic coordinate system 
 in yellow (all three as in Figure 1), the planet-based hybrid equatorial coordinate system 
 in olive, and the planet-based topographic coordinate system 
 in maroon.
Coordinate systems 
, 
, and 
 were described in detail earlier.

 represents the primary intersection of the equator of the Earth and the equator of the planet. This is the basis for the planet rotation angle 
 used by the IAU.

 represents the primary location on the planet (with topographic latitude and longitude equal to zero). It rotates around the rotation axis of the planet (
) once per sidereal planet day. The planet's prime meridian is indicated by the dashed (front) and dotted (rear) ellipse passing through 
 and 
.

The sidereal time angle θ
 says which planet-based right ascension passes through the zenith of an observer at the planet's prime location 
 at the epoch (
). This is equal to how far the planet's prime location 
 has rotated beyond the planet-based primary direction 
.

The angle υθ
 says how far the planet's primary intersection 
 is beyond the planet's primary direction 
.

How do you calculate θ
?


 are the coordinates of the primary intersection of the planet's equator with the Earth's equator, expressed in Earth-based equatorial coordinates, as defined by the IAU.

α
α
 
α
α
 
 


For Mars we find
 

which corresponds to Earth-based ecliptic longitude 47.6814° and latitude 0°.


 are the coordinates of the primary intersection expressed in Earth-based ecliptic coordinates.

ε
α
ε
α
ε
α
 
 


For Mars we find
 

which corresponds to Earth-based ecliptic longitude 45.21813 and latitude −17.1049°.



We now have the x axis 
 of the planet-based hybrid equatorial coordinate system, in terms of Earth-based ecliptic coordinates. The z axis is the planet's primary pole

 


Now that we have two of the three axes, we can complete the coordinate basis 
 for converting from Earth-based ecliptic coordinates to planet-based hybrid equatorial coordinates.


 
 


Normalization of the cross product isn't needed here because the two factors are perpendicular and of length 1.



For Mars we find
 
 
 
 
 
 



 are the coordinates of the primary intersection of the equators of Earth and the planet, expressed in planet-based equatorial coordinates.

 


The z coordinate of that vector is equal to 0 by definition, because the primary intersection is on the celestial equator of the planet.



For Mars we find
 
 
 

which corresponds to Earth-based ecliptic longitude 136.7527° and latitude 0°.


υ
 is the planet-based right ascension of the primary intersection 
. This is the angle between the planet's primary direction and the primary intersection.

υ
 


For Mars we find
υ


Calculate θ
 from 
 and υ
:

θ
υ
 


For Mars we find
θ

18. Planet-based Topographic Coordinates
We derive the transformation of coordinates between planet-based topographic coordinate system 
 and Earth-based ecliptic coordinate system 
. Coordinate system 
 is the only one of the coordinate systems discussed on this page that rotates. The orientation of coordinate system 
 depends on the rotation angle 
.

 

The conversion from 
 to 
 goes as follows.


 are the coordinates of the primary pole of the planet, expressed in Earth-based ecliptic coordinates.
 


 are the coordinates of the primary location on the planet, expressed in Earth-based ecliptic coordinates. We find 
 from 
 by converting that to the planet-based equatorial coordinate system 
, then rotating over angle 
 around the z axis of 
 (the rotation axis of the planet), and then converting back to the Earth-based ecliptic coordinate system.
 
 



We find 
 from 
 and 
, which are at right angles.
 


Now that we have all three coordinate axes we can construct the coordinate basis 
 for converting Earth-based ecliptic coordinates to planet-based topographic coordinates, and 
 for the opposite conversion.

 
 

 and 
 vary with time, because they depend on 
 which varies with time.

For Mars for 
 we find

 
 
 
 
 
 
 
 
 
 
 
 

19. Planet-Bound Horizontal Coordinates
From the planet-bound topographic coordinate system 
 we can deduce the planet-bound horizontal coordinate system 
 for location 
 at latitude 
 and longitude 
 on the planet, assuming the planet to be spherical. If the primary pole is the north pole, then 
 indicates north latitude and 
 indicates east longitude.

The z axis of 
 points to the zenith of an observer at that location. The y axis of 
 at that location points along the surface toward the primary pole. The x axis is perpendicular to the y and z axes and makes the coordinate system right-handed. If the north pole is the primary pole, then the x axis points to the east.

Coordinate system 
 is equal to coordinate system 
 of the primary location on the planet, except that the x axis of 
 is the z axis of 
, the z axis of 
 is the y axis of 
, and the y axis of 
 is the x axis of 
.

For location 
 on the planet, we find 
 from 
 as follows:


Rotate by angle 
 around the z axis so that the x axis points from the center of the planet to the meridian of the location.

Rotate by angle 
 around the y axis such that the x axis points from the center of the planet to the location.

Substitute the axes as follows: 
.
Combined this yields:

 
 
 
 

Let's check. The direction from the center of the planet to location 
 at latitude 
 and longitude 
 is, expressed in topographic coordinates:

 
 

and that is indeed equal to the z axis (3rd column) of 
, so equal to the direction of the zenith of point 
, expressed in topographic coordinates.

20. Coordinate Transformations Between All Systems

Fig. 3: Coordinate Transformations
Fig. 3: Coordinate TransformationsWe now have a network of coordinate transformations with which we can convert coordinates from each of the coordinate systems 
 to all other systems from that list. Figure 3 shows the coordinate systems. A line between two coordinate systems in that figure means we have described a direct transformation between the two systems, earlier on this page. By combining two direct transformations we can connect all other coordinate systems with each other. For example:



21. North Pole, South Pole, Primary Pole
The poles of a celestial object are the two intersections of the rotation axis of the object and the surface of the object. We call one of them the north pole and the other one the south pole, but which is which? Not everyone agrees about that. There are basically two ways to choose which pole is the north pole and which is the south pole:


The north pole is the pole that points to the north side of a suitable fundamental plane.

The north pole is the pole seen from above which the object appears to rotate counterclockwise around its axis.
The IAU uses the first of these definitions for the planets of the Solar System. The fundamental plane is then the Invariable Plane of the Solar System ― of which the poles are the average rotation axis for all orbital motion and rotation in the Solar System. The north side of that plane is the side to which the north pole of the Earth points. According to the IAU definition, the north poles of all planets of the Solar System point to the north side of the Invariable Plane of the Solar System.

For all other celestial objects, including dwarf planets and asteroids and comets in our Solar System, and for all celestial objects outside of the Solar System (including major planets), the IAU uses the second definition.

The great advantage of the second definition is that you need to look only at the celestial object itself, and do not need to search for some fundamental plane with which to compare the pole. The orientation of the fundamental plane might nog even be known very accurately.

For most planets in the Solar System both definitions yield the same result ― but not for Venus and Uranus. If you look down from above the IAU-defined north pole of Venus or Uranus, then those planets rotate clockwise around their axis.

Both definitions don't yield the same results for Pluto, either. Until 2006, Pluto was regarded as a small one of the major planets, to which the first definition applies, so the pole of Pluto that points to the north of the Invariable Plane of the Solar System was regarded as its north pole. Since 2006, Pluto is regarded as a dwarf planet, to which the second definition applies, so now the pole from above which Pluto appears to rotate in the counterclockwise direction is regarded as its north pole ― which is the opposite pole from what was regarded as its north pole until 2006.

Because it isn't always clear which pole is the north pole, I use the term "primary pole" instead in the above text, so that you cannot be led astray if your idea of which pole is the north pole were different from mine.


[AA]
languages: [en] [nl]

//aa.quae.nl/en/reken/zonpositie.html;
Last updated: 2021-07-19
```

Now, why is there a discrepency bw my code results and astrals? Please make mine more accurate (as accurate as possible). Keep in mind that acos, atan, etc.. all math functions are proper with 100% accuracy, so they are not the problem ( i tested them against math module, it gets 100% accuracy).
Return full code.