ISS brightness

[sgbmzm]

I followed the beautiful things that are presented on the following page: https://rhodesmill.org/skyfield/earth-satellites.html One thing I was missing: the ISS brightness estimate (as an example) This is available on the website https://www.heavens-above.com/ and I think it is important information. Although it may not be possible to calculate it precisely, but even a general estimate would be better than nothing. So I tried looking here https://stackoverflow.com/questions/19739831/is-there-any-way-to-calculate-the-visual-magnitude-of-a-satellite-iss And here: https://astronomy.stackexchange.com/questions/28744/calculating-the-apparent-magnitude-of-a-satellite/28765#28765 But I actually couldn't figure out what to do. Can you add such a feature to Skyfield, or an example of how to do it? Thanks!

[brandon-rhodes]

To clarify the question: you would like the formula shown in that second link, translated from C# to Python, because it's probably a similar calculation to the one used (but apparently not documented?) by Heavens Above?

[mworion]

Hi, I did some work on that, please look here:

https://github.com/mworion/MountWizzard4/blob/master/mw4/gui/mainWmixin/tabSat_Search.py

It‘s far from prefect, but I hope this might help.

[sgbmzm]

What I tried to do is: I used chat.openai.com to translate into Python what is presented there in the C# language. The result I got is:

import math

distance_to_satellite = 485  # This is in KM
phase_angle_degrees = 113.1  # Angle from sun->satellite->observer
pa = phase_angle_degrees * 0.0174533  # Convert the phase angle to radians
intrinsic_magnitude = -1.8  # -1.8 is std. mag for iss

term_1 = intrinsic_magnitude
term_2 = 5.0 * math.log10(distance_to_satellite / 1000.0)

arg = math.sin(pa) + (math.pi - pa) * math.cos(pa)
term_3 = -2.5 * math.log10(arg)

apparent_magnitude = term_1 + term_2 + term_3

But I don't know what about phase_angle_degrees = 113.1 # Angle from sun->satellite->observer How to set Angle from sun->satellite->observer in Skyfield.

[brandon-rhodes]

Try the phase_angle() method, and see whether it works:

https://rhodesmill.org/skyfield/api-position.html#skyfield.positionlib.ICRF.phase_angle

[mworion]

In my code example I used it too.Am 02.02.2024 um 15:33 schrieb Brandon Rhodes @.***>: Try the phase_angle() method, and see whether it works: https://rhodesmill.org/skyfield/api-position.html#skyfield.positionlib.ICRF.phase_angle

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you commented.Message ID: @.***>

[brandon-rhodes]

In my code example I used it too.

Thank you for confirming that phase_angle() seems to work for this application! If @sgbmzm also gets good results, then I'll be happy to add a generic magnitude routine to Skyfield using it. We will see what they report back.

[mworion]

I have to correct myself, this was an older try, actually I'm calculating the visible phase like:

def calcSatSunPhase(sat, loc, ephemeris, tEv):
        """
        https://stackoverflow.com/questions/19759501
            /calculating-the-phase-angle-between-the-sun-iss-and-an-observer-on-the
            -earth
        https://space.stackexchange.com/questions/26286
            /how-to-calculate-cone-angle-between-two-satellites-given-their-look-angles
        :param sat:
        :param loc:
        :param ephemeris:
        :param tEv:
        :return:
        """
        earth = ephemeris['earth']

        vecObserverSat = (sat - loc).at(tEv)
        vecSunSat = (earth + sat).at(tEv)
        phase = vecObserverSat.separation_from(vecSunSat)
        return phase

Which might be the same as the phase_angle now. I implemented it some time about and skyfield improved since then... Michel

[sgbmzm]

Unfortunately I could not use phase_angle() I get an error:

ValueError: you can only observe() a body whose vector center is the Solar System Barycenter, but this vector has the center 399 EARTH

from skyfield.api import load, wgs84
ts = load.timescale()
stations_url = 'http://celestrak.org/NORAD/elements/stations.txt'
satellites = load.tle_file(stations_url)

by_name = {sat.name: sat for sat in satellites}
iss = by_name['ISS (ZARYA)']

eph = load('de441S.bsp')
bluffton = wgs84.latlon(+40.8939, -83.8917)
t0 = ts.utc(2024, 1, 23)
t1 = ts.utc(2024, 1, 24)

t, events = iss.find_events(bluffton, t0, t1, altitude_degrees=30.0)
iss_phase_angle = eph['earth'].at(t).observe(iss).apparent().phase_angle(eph['sun'])
iss_phase_angle

[brandon-rhodes]

Interesting! I guess phase_angle() has only ever been used for planets before, and not Earth satellites.

The way to make that message go away is to .observe(eph['earth'] + iss) so that Skyfield knows the whole vector between the Solar System Barycenter and the ISS. Which, to my surprise, then lead to another error:

ValueError: zero-size array to reduction operation maximum which has no identity

So I ran print(t) and saw:

<Time tt=[]>

It turns out that .observe() didn't like computing for no specific times. So I've just committed a fix, so that .observe() and .apparent() both successfully compute empty position arrays when given an empty time. (Just for completeness; not because it really helps your situation, of course.)

We can tweak your script to involve an actual position by pivoting to just using t0 for now:

iss_phase_angle = eph['earth'].at(t0).observe(eph['earth'] + iss).apparent().phase_angle(eph['sun'])
print(iss_phase_angle)

And I get:

104deg 23' 09.7"

So maybe you can proceed with that. But I really don't like making people .observe() a satellite, because it incurs so much needless expense, so while you work with that angle to see if it gives you a good magnitude, I'll think through how .phase_angle() might be made to work even if the position is centered at the Earth.

[sgbmzm]

Thanks. I tried. the brightness does not match what is written on https://www.heavens-above.com/ So maybe the code in the C# language really doesn't match the formula of https://www.heavens-above.com/

[brandon-rhodes]

Does Heavens Above publish its math anywhere, or have a list of citations from which you could maybe find which source they draw their equations?

[sgbmzm]

Unfortunately I don't know. But at least, thanks to you, I could test what was brought in the C# language, that was also an important thing.

[brandon-rhodes]

Out of curiosity, what is the H.A. magnitude, and which did your routine compute? Could you share your Python code snippet? I'd be happy to read over it and double-check that everything is being called correctly.

[sgbmzm]

Here is a simple example that follows the examples given in https://rhodesmill.org/skyfield/earth-satellites.html plus the calculation translated from the C# language Set the location as you wish, and compare against heavens-above

from skyfield.api import load, wgs84
import math

ts = load.timescale()
stations_url = 'http://celestrak.org/NORAD/elements/stations.txt'
satellites = load.tle_file(stations_url)

by_name = {sat.name: sat for sat in satellites}
iss = by_name['ISS (ZARYA)']

eph = load('de441.bsp')
bluffton = wgs84.latlon(+40.8939, -83.8917)
t0 = ts.utc(2024, 2, 4)
t1 = ts.utc(2024, 2, 8)

t, events = iss.find_events(bluffton, t0, t1, altitude_degrees=30.0)
sunlit = iss.at(t).is_sunlit(eph)
event_names = 'rise above 30°', 'culminate', 'set below 30°'

for ti, event, sunlit_flag in zip(t, events, sunlit):
    name = event_names[event]
    state = ('in shadow', 'in sunlight')[sunlit_flag]
    iss_topocentric = (iss - bluffton).at(ti)
    alt, az, distance = iss_topocentric.altaz()
    iss_phase_angle = eph['earth'].at(ti).observe(eph['earth'] + iss).apparent().phase_angle(eph['sun'])
    distance_to_satellite = distance.km  # This is in KM
    iss_phase_angle = eph['earth'].at(ti).observe(eph['earth'] + iss).apparent().phase_angle(eph['sun'])
    pa = iss_phase_angle.radians  # Convert the phase angle to radians
    intrinsic_magnitude = -1.8  # -1.8 is std. mag for iss
    term_1 = intrinsic_magnitude
    term_2 = 5.0 * math.log10(distance_to_satellite / 1000.0)
    arg = math.sin(pa) + (math.pi - pa) * math.cos(pa)
    term_3 = -2.5 * math.log10(arg)
    apparent_magnitude = term_1 + term_2 + term_3
    magnitude = apparent_magnitude

    print(ti.utc_strftime('%Y %b %d %H:%M:%S'), name, state,"; magnitude:", round(magnitude,1))