skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.4k stars 211 forks source link

How to calculate the brightness of the satellite at a specified time? #660

Closed lingduxingxi closed 2 years ago

lingduxingxi commented 2 years ago

Predicting future satellite brightness: But I found it impossible,The code is as follows:

#!/usr/bin/python
# -*- coding: UTF-8 -*-

import math,ephem,datetime,time
from datetime import datetime
def apparent_magnitude_test():
    distanceToSatellite = iss.range / 1000 #This is in KM
    #phaseAngleDegrees = phase_angle #Angle from sun->satellite->observer
    pa = phaseAngleDegrees * 0.0174533#Convert the phase angle to radians
    intrinsicMagnitude = -1.8#-1.8 is std. mag for iss
    term_1 = intrinsicMagnitude
    term_2 = 5.0 * math.log(distanceToSatellite / 1000.0,10)
    arg = math.sin(pa) + (math.pi - pa) * math.cos(pa)
    term_3 = -2.5 * math.log(arg,10)
    apparentMagnitude = term_1 + term_2 + term_3
    print(apparentMagnitude)

def apparent_magnitude():
    distanceToSatellite = 2267#This is in KM
    # phaseAngleDegrees = 113.1#Angle from sun->satellite->observer
    pa = phaseAngleDegrees * 0.0174533#Convert the phase angle to radians
    intrinsicMagnitude = 3.2#-1.8 is std. mag for iss
    term_1 = intrinsicMagnitude
    term_2 = 5.0 * math.log(distanceToSatellite / 1000.0,10)
    arg = math.sin(pa) + (math.pi - pa) * math.cos(pa)
    term_3 = -2.5 * math.log(arg,10)
    apparentMagnitude = term_1 + term_2 + term_3
    # Was this correct? Go to the shell window to check:
    print(apparentMagnitude)

def phase_angle():
    # SSA Triangle.  We have side a and b and angle C.  Need to solve to find side c 
    au=149597870.7 #km
    a = sun.earth_distance * au - ephem.earth_radius/1000 #distance sun from observer (Km)
    b = iss.range / 1000 # distance to ISS from observer (Km)
    sun.compute(cadi)
    angle_c = ephem.separation( (iss.az, iss.alt), ( sun.az, sun.alt) ) 
    c = math.sqrt( math.pow(a,2) + math.pow(b,2) - 2*a*b*math.cos( angle_c) )
    # now we find the "missing" angles (of which angle A is the one we need)
    angle_a = math.acos((math.pow(b,2) + math.pow( c,2) - math.pow(a,2)) / (2 * b * c)) 
    angle_b = math.pi - angle_a - angle_c #note: this is basically ZERO - not a big surprise really - and I don't need this anyway.
    # phase_angle = angle_a # This is the angle we need.  BINGO!!
    phaseAngleDegrees=angle_a
    print(phaseAngleDegrees)
    return phaseAngleDegrees

# date = datetime.now()
from datetime import datetime
date = datetime(2021, 11, 18, 11, 29, 57)

print(date)
sun=ephem.Sun()
sun.compute(date)
# print(sun.compute(date))

satellite_name="TIANHE"             
satellite1="1 48274U 21035A   21300.07730318  .00016357  00000-0  19379-3 0  9995"
satellite2="2 48274  41.4698 198.7958 0006485 267.0658 229.4040 15.60844829 28326"
iss=ephem.readtle(satellite_name,satellite1,satellite2)

# print(iss)
cadi=ephem.Observer()
cadi.lon,cadi.lat,cadi.elevation='102.9505','24.5801',1722
iss.compute(cadi)
print(iss)
# phaseAngleDegrees=0
phaseAngleDegrees = phase_angle()
apparent_magnitude()

What went wrong?

brandon-rhodes commented 2 years ago

Thank you for the edits! The formatting of your program is much improved — I can read it very easily now — and it now sounds more friendly. I think that contributors (me or someone else) will soon have time to try running your code, now that we can see the format (which is so important with Python, because it pays attention to indent!).

The comments in the code will help us understand what you are doing. I will reply when I have the chance!

(Oh — it looks like you replied twice, by accident, so I will remove one of the two replies.)

(And now that your question is in better shape, I’ll delete our comments about the earlier version, so that contributors simply focus on your question as they come and read!)

bmxp commented 2 years ago

In far does your code relate to skyfield? I just see ephem used there.

brandon-rhodes commented 2 years ago

I have had the chance to run the script now — it prints:

2021-11-18 11:29:57
<ephem.EarthSatellite 'TIANHE' at 0x7f80f05d0760>
2.88979102172
3.7357496733

What magnitude were you expecting? Is 3.7 more dim than expected, or more bright?

brandon-rhodes commented 2 years ago

(Also, in case it helps — but maybe you have already seen this:)

https://space.stackexchange.com/questions/5142/how-to-calculate-the-brightness-of-a-passing-satellite

lingduxingxi commented 2 years ago

(Also, in case it helps — but maybe you have already seen this:)

https://space.stackexchange.com/questions/5142/how-to-calculate-the-brightness-of-a-passing-satellite

thank you

brandon-rhodes commented 2 years ago

I am going to close this, as several questions were never answered, and without them the issue won't be able to make further progress. But feel free to make further comments once the issue is closed.