Further to KH

'''
This is BR_Mods2.py
Sequel to https://pastebin.com/YT28KrG0 (re BR_Mods.py) that should be read first.

My beginners learning curve continues. I have found the Skyfield method "from_altaz". It answers my comment
"Whilst it must be possible to use altaz() to first create hypothetical body positions that compensate for refraction..."

I have added code (lines 86-90) to create such hypothetical positions.
Skyfield's "separation_from" method then calculates the arc distance perfectly.
The traditional code is on lines 74-84

Output compares KH's original, BR's conventional modifications and now the revised Skyfield code.
The latter two agree.


line 16
'''


import skyfield.api, skyfield.data.hipparcos
import numpy, datetime, math, argparse, sys
SkyfieldLoad = skyfield.api.Loader('~/skyfield-data')
Ephemeris = SkyfieldLoad('de421.bsp') # 16.7 MB download
with SkyfieldLoad.open(skyfield.data.hipparcos.URL) as f0:
    hippa = skyfield.data.hipparcos.load_dataframe(f0)
from skyfield.api import wgs84, Angle
SkyfieldTimescale = SkyfieldLoad.timescale() # timescale object with built-in UT1-tables


def Calculate(Arguments) :
    EventUTC, GeoPos, Body1Name, Body2Name = Arguments  # Unpack Arguments

    # Assign and show time
    t = SkyfieldTimescale.ut1(EventUTC.year, EventUTC.month, EventUTC.day,
            EventUTC.hour, EventUTC.minute, EventUTC.second + EventUTC.microsecond / 1e6)
    print(); print(); print(f't : {t.ut1_strftime()}')
    print('================================')

    # Observer Position
    ObsBodycentricPos = NameToBodyData('earth')     # Observer at Earth centre
    if GeoPos:
        ObsPos = ObsBodycentricPos + GeoPos         # Observer at specified 3D position (height relative to earth's WGS84 surface)
    else:
        ObsPos = ObsBodycentricPos
    ObsBaryPos = ObsPos.at(t)                       # Barycentric Observer BCRS position

    # First body
    Body1Data = NameToBodyData(Body1Name)
    Apparent1 = ObsBaryPos.observe(Body1Data).apparent()
    print(f'Body1Name : {Body1Name}')
    print(f'RA : {Apparent1.radec("date")[0]}')
    print(f'Dec : {Apparent1.radec("date")[1]}')
    if GeoPos:
        print(f'Azimuth : {Apparent1.altaz("standard")[1]}')
        print(f'Altitude : {Apparent1.altaz("standard")[0]}')
    print('================================')

    # Second body
    Body2Data = NameToBodyData(Body2Name)
    Apparent2 = ObsBaryPos.observe(Body2Data).apparent()
    print(f'Body2Name : {Body2Name}')
    print(f'RA : {Apparent2.radec("date")[0]}')
    print(f'Dec : {Apparent2.radec("date")[1]}')
    print(f'SD : {SemiDiameter(Apparent2).arcminutes():.4f}')
    if GeoPos:
        print(f'Azimuth : {Apparent2.altaz("standard")[1]}')
        print(f'Altitude : {Apparent2.altaz("standard")[0]}')
    print('================================')

    # Get refracted values (dist unused)
    Ha1, Az1, dist1 = Apparent1.altaz("standard")
    Ha2, Az2, dist2 = Apparent2.altaz("standard")

    # Calculate arc between centres - traditional method
    ArcBetweenCentres = math.degrees(
        (math.acos(
            math.sin(Ha1.radians) *
            math.sin(Ha2.radians) +
            math.cos(Ha1.radians) *
            math.cos(Ha2.radians) *
            math.cos(Az1.radians - Az2.radians)
            )
        )
    )

    # Calculate arc between centres - revised Skyfield method
    # Hypothetical positions to compensate for refraction
    HypoRefrPos1 = ObsBaryPos.from_altaz(alt_degrees=Ha1.degrees, az_degrees=Az1.degrees)
    HypoRefrPos2 = ObsBaryPos.from_altaz(alt_degrees=Ha2.degrees, az_degrees=Az2.degrees)
    Separation = HypoRefrPos1.separation_from(HypoRefrPos2).degrees

    # Show Arc between centre results
    print(f'Original code:    Arc between centres: {Apparent1.separation_from(Apparent2)}')
    print('Trad method:      Arc between centres:', skyfield.units.Angle(degrees=ArcBetweenCentres))
    print('Revised Skyfield: Arc between centres:', skyfield.units.Angle(degrees=Separation)) ; print()

    # Show arc between near limbs
    # Original code
    NearArc1 = Apparent1.separation_from(Apparent2).degrees - SemiDiameter(Apparent1).degrees - SemiDiameter(Apparent2).degrees
    # Trad arc
    NearArc2 = ArcBetweenCentres - SemiDiameter(Apparent1).degrees - SemiDiameter(Apparent2).degrees
    # Revised arc
    Near3Arc = Separation - SemiDiameter(Apparent1).degrees - SemiDiameter(Apparent2).degrees
    print('Original code:    Arc between near limbs: ', skyfield.units.Angle(degrees=NearArc1))
    print('Trad method:      Arc between near limbs: ', skyfield.units.Angle(degrees=NearArc2))
    print('Revised Skyfield: Arc between near limbs: ', skyfield.units.Angle(degrees=Near3Arc)) ; print()


#### Other Functions follow ####
def SemiDiameter(BodyData): # Apparent
    if BodyData.target == 10: # sun
        radius = 696340.0 # km
    elif BodyData.target == 301: # moon
        # radius = 1737.1 # km
        radius = 1737.92
    else:
        radius = 0.0
    # print(f'## BodyData.distance().km : {BodyData.distance().km}')
    return skyfield.units.Angle(radians=numpy.arcsin(radius / BodyData.distance().km))

def NameToBodyData(Name) :
    Name = Name.lower()
    if Name in ['sun', 'earth', 'moon']:
        BodyData = Ephemeris[Name] # bodycentric
    elif Name in ['venus', 'mars', 'jupiter', 'saturn']:
        BodyData = Ephemeris[Name + ' barycenter'] # barycentric (why?)
    else:
        if Name == 'hamal': HipNumber = 9884
        elif Name == 'aldebaran': HipNumber = 21421
        elif Name == 'betelgeuse': HipNumber = 27989
        elif Name == 'sirius': HipNumber = 32349
        elif Name == 'pollux': HipNumber = 37826
        elif Name == 'regulus': HipNumber = 49669
        elif Name == 'spica': HipNumber = 65474
        elif Name == 'antares': HipNumber = 80763
        elif Name == 'altair': HipNumber = 97649
        elif Name == 'fomalhaut': HipNumber = 113368
        elif Name == 'markab': HipNumber = 113963
        elif Name == 'zubenelgenubi': HipNumber = 72622
        BodyData = skyfield.api.Star.from_dataframe(hippa.loc[HipNumber])
        # UnboundLocalError if Name not in above
    return BodyData

def ParseArguments(Arguments):
    Arguments = vars(Arguments) # usings vars() to create dict
    # --time
    if Arguments['time'] is None:
        EventUTC = datetime.datetime.utcnow()
    else:
        EventUTC = datetime.datetime.strptime(Arguments['time'], '%Y-%m-%dT%H:%M:%S')
    # --geo
    if Arguments['geo']:
        GeoArgs = Arguments['geo'].split(',')
        # use "parse" ? : r0 = '{lat:g},{:s}{lng:g}{}' # ignore Altitude
        # ITRS ECEF # XXX height used for "dip" calc (in altaz) ?
        GeoPos = skyfield.api.wgs84.latlon(float(GeoArgs[0]), float(GeoArgs[1]), float(GeoArgs[2]))
    else:
        GeoPos = None
    # -n1, n2
    Body1Name = Arguments['Body1Name']
    Body2Name = Arguments['Body2Name']
    # Do the hard work elsewhere
    return EventUTC, GeoPos, Body1Name, Body2Name

def main():
    parser = argparse.ArgumentParser(description="almanac")
    parser.add_argument("-n1", "--Body1Name", default = 'Aldebaran', type = str, help="heavenly body")  # was -n
    parser.add_argument("-n2", "--Body2Name", default = 'Moon', type = str, help="heavenly body")       # was -N
    parser.add_argument("-t", "--time", default = '2022-03-05T03:00:00', help="time : %Y-%m-%dT%H:%M:%S")
    parser.add_argument("-g", "--geo", default = '38.0,-122.0,0', help="assumed position : lat,lng,meters") # SF area
    parser.add_argument("-T", "--temp", default = "standard", help="temperature_C (for refraction)")
    Calculate(ParseArguments(parser.parse_args()))

if __name__ == '__main__':
    sys.exit(main())

'''
Scratch pad for copy/paste command lines
python BR_Mods2.py -t 2022-03-05T03:00:00 -n1 Aldebaran -n2 Moon -g 38.0,-122.0,0
python BR_Mods2.py -t 2022-02-22T11:05:48 -n1 Zubenelgenubi -n2 Moon -g 32.0,-109,1400
python BR_Mods2.py -t 2022-02-22T12:27:22 -n1 Zubenelgenubi -n2 Moon -g 32.0,-109,1400
'''