micropython-samples icon indicating copy to clipboard operation
micropython-samples copied to clipboard

I wrote new function for sun_moon: alt_az_ra_dec()

Open sgbmzm opened this issue 11 months ago • 1 comments

Further to some of what I asked for https://github.com/peterhinch/micropython-samples/issues/42

At the beginning of the sun_moon file, need to add atan2, degrees, asin to the math import.

# A function I built based on the sin_alt function, together with AI and examples from the Astral library
# The function returns the altitude of the Sun or Moon in degrees, azimuth from north in degrees, right ascension in hours (decimal - without minutes and seconds), and declination in degrees
# The function also calculates the local sidereal time and the hour angle of the Sun or Moon, but these are not returned

def alt_az_ra_dec(self, hour, sun=True):
    """
    Returns the altitude (Alt) and azimuth (Az) of the Sun in degrees.
    """
    func = minisun if sun else minimoon
    mjd = (self.mjd - 51544.5) + hour / 24.0
    t = mjd / 36525.0
    x, y, z = func(t)  # Cartesian coordinates of the Sun or Moon

    tl = self.lstt(t, hour) + self.long  # Local sidereal time
    sin_alt = self.sglat * z + self.cglat * (x * cos(radians(tl)) + y * sin(radians(tl)))
    alt = degrees(asin(sin_alt))  # Altitude of the Sun in degrees
    
    # Declination calculation in degrees
    rho = sqrt(x * x + y * y)  # Projection of the vector on the XY plane
    dec = degrees(atan2(z, rho))  # Declination calculation in degrees
    
    # Right ascension (RA) calculation
    ra = ((48.0 / (2 * pi)) * atan(y / (x + rho))) % 24  # Right ascension in hours, in a decimal fracture

    # Azimuth (Az) calculation
    hourangle = radians(tl) - radians(ra * 15)  # Local sidereal time minus the right ascension of the star gives its hour angle (ra * 15 converts to degrees)
    hourangle_hours = (degrees(hourangle) % 360) / 15.0  # Hour angle in hours, in a decimal fracture
    
    sh = sin(hourangle)
    ch = cos(hourangle)
    sd = sin(radians(dec))
    cd = cos(radians(dec))
    sl = self.sglat  # ==sin(radians(lat))
    cl = self.cglat  # ==cos(radians(lat))

    x = -ch * cd * sl + sd * cl
    y = -sh * cd
    az = degrees(atan2(y, x)) % 360  # Azimuth in degrees
    
    return alt, az, ra, dec

use:

# Calculating the current hour in a decimal fracture
#current_hour = (hour + (minute / 60) + (second / 3600)) - utc_offset_hours
#s_alt, s_az, s_ra, s_dec = riset.alt_az_ra_dec(current_hour, sun=True)
#m_alt, m_az, m_ra, m_dec = riset.alt_az_ra_dec(current_hour, sun=False)

sgbmzm avatar Feb 24 '25 19:02 sgbmzm

In reply to your various posts. From astronomy/README.md

I am not an astronomer. If there are errors in the fundamental algorithms I am unlikely to be able to offer an opinion, still less a fix.

From the main README.md:

[samples] are intended as pointers for programmers rather than being complete solutions. Egregious bugs will be fixed but I may not accept feature requests.

Re astronomy, my source was ""Astronomy on the Personal Computer" by Montenbruck and Pfleger". If you can get your hands on a copy it may have the answer to some of your queries.

peterhinch avatar Feb 25 '25 13:02 peterhinch