Response to KH

''' Ken .. see also https://pastebin.com/Rqb1YBsJ .. late addition.
Modified by 1billritchie@gmail.com from Ken Harris' post https://pastebin.com/5k8bwnHX .
Written for Ken Harris' eyes, but nothing private herein.
Skyfield V1.41 used .... my pip says that 1.41 is latest! (But 1.42 available on some installations)

This is my first python attempt, other than Hello World's, etc to grasp the absolute basics.
Not really mine ... mainly rearrangements of Ken's work.
As a python beginner, have changed variables, etc to more verbose ones solely to help my understanding of the code.

Purpose is to quantify a difference in arc distance outputs between this and two comparisons,
    (a) https://friendsofthevigilance.org.uk/Astron/Astron.html
    (b) http://fer3.com/arc  /  DATA  /  Analyze Lunars Online
    Both written independently, the latter by Frank Reed, probably the world's leading expert on lunars.
    Both usually agree within 0.1 arc minutes.

OBSERVATIONS.
    1.
    Skyview arc distances differ from comparisons.
    (Trend to increase difference with increasing arc values, with Markab as exception to this).
    Below table illustrates this.
    NB The value of Moon radius has been changed trivially to agree with Astron  (line 154)

      Near Limbs Arc Comparison. 2022 March 5th at 03:00:00 from 38.0N 122.0W 0m altitude
      ===================================================================================
                    	                    Skyfield 	        	           Fer3
    Moon - Body	         Skyfield          reformatted        Astron      Lunars Analyzer
    Regulus	          136deg 13' 01.0"	    136° 13.0'		136° 07.9'		136° 07.8'
    Pollux	          100deg 00' 19.1"	    100° 00.3'	    099° 56.4'		099° 56.4'
    Betelgeuse	       74deg 48' 41.0"	    074° 48.7'		074° 45.7'		074° 45.6'
    Aldebaran	       56deg 05' 19.0"	    056° 05.3'		056° 02.4'		056° 02.3'
    Markab	           29deg 49' 53.1"	    029° 49.9'	    029° 43.9'		029° 43.8'

    Markab - Regulus  149deg 27' 35.8"	    149° 27.6'	    149° 10.9'		n/a
    Markab - Regulus  149deg 27' 35.8"	    149° 27.6'	    149° 27.6'		n/a
    (geocentric)
    ======================================================================================

    2.
    Markab - Regulus, geocentric. Agrees with comparisons.
    Markab - Regulus, topocentric. Skyfield result same as geocentric! Is refraction the problem?
    Chosen Markab low alt & azimuths opposite, so refraction subtractive from arc distance.
    Markab (altitude 2° 13') appears 16.1' higher due refraction.
    Regulus (altitude 25° 43') appears 2.1' higher.
    Combined refraction is 18.2'. This reduces Skyfield's 149° 27.6' to 149° 9.4', much closer to desired 149° 10.9'

    3.
    Changing the altitude argument from 0m (or 6m) to (say) 1000m has (correctly) no effect with star pair,
    but correct small (about 0.5') with Moon - Star pair.  Again suggests refraction not working WITH ARC calculation.

    4.
    Refraction IS WORKING with altaz().
    EG: print(f'Altitude : {Apparent1.altaz("standard")[0]}') gives correct refracted values.

    5.
    "Apparent1.separation_from(Apparent2)" gives arc from observer position without refraction, seemingly working
    by vector math of observer and two body positions. Whilst it must be possible to use altaz() to first
    create hypothetical body positions that compensate for refraction, I think my traditional code,
    starting at line 122, is much easier. Results agree with comparisons.

    6.
    Skyfield's fantastic accuracy claims are probably correct for the effects considered. However, I
    have yet to find out what effects are covered. There is more to accuracy than calculation accuracy!
    Does aberration also include the observer's (say 1400km/h) velocity on a rotating Earth? Do the WGS84 methods
    compensate for parallax in AZIMUTH?  How can the accuracy be specified if the user can change body radii?
    Does it use barycentric proper motion for binary stars, particularly Rigil Kentaurus and Sirius? Does it
    include gravitational deflections of the vertical which, on Earth, are far greater than the included relativistic effects?
    When I find answers to such questions, I will use it to improve Astron's methods.

line 69
'''
#### Load essentials ####
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

#### Calculation function, called from main(). Was within ParseArguments().
def Calculate(Arguments) :
    # Unpack Arguments
    EventUTC, GeoPos, Body1Name, Body2Name = 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
    Observer = NameToBodyData('earth')      # Observer at Earth centre
    if GeoPos: Observer += GeoPos           # Observer at specified 3D position (height relative to earth's WGS84 surface)
    Barycentric = Observer.at(t)            # Barycentric Observer BCRS position

    # First body
    Body1Data = NameToBodyData(Body1Name)           # print(f'# Body1Data : {Body1Data}')
    Astrometric1 = Barycentric.observe(Body1Data)   # Corrected for light time
    Apparent1 = Astrometric1.apparent()             # relativistic effects and aberration -> Apparent GCRS position
    print(f'Body1Name : {Body1Name}')
    print(f'RA : {Apparent1.radec("date")[0]}')
    print(f'Dec : {Apparent1.radec("date")[1]}')
    print(f'SD : {SemiDiameter(Apparent1).arcminutes():.4f}')
    if GeoPos:
        print(f'Azimuth : {Apparent1.altaz("standard")[1]}')
        print(f'Altitude : {Apparent1.altaz("standard")[0]}')
    print('================================')

    # Second body
    Body2Data = NameToBodyData(Body2Name) # print(f'# Body2Data : {Body2Data}')
    Astrometric2 = Barycentric.observe(Body2Data) # Astrometric ICRS / ICRF position and velocity
    Apparent2 = Astrometric2.apparent() # deflection & aberration : Apparent GCRS position and velocity
    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]}')

    # ADDED CODE STARTS HERE
    # Calculate arc between centres - traditional method
    Ha1 = Apparent1.altaz("standard")[0].radians
    Az1 = Apparent1.altaz("standard")[1].radians
    Ha2 = Apparent2.altaz("standard")[0].radians
    Az2 = Apparent2.altaz("standard")[1].radians
    ArcBetweenCentres = math.degrees(
        (math.acos(
            math.sin(Ha1) *
            math.sin(Ha2) +
            math.cos(Ha1) *
            math.cos(Ha2) *
            math.cos(Az1 - Az2)
            )
        )
    )
    print(f'Skyfield:    Arc between centres: {Apparent1.separation_from(Apparent2)}')
    print('Trad method: Arc between centres:', skyfield.units.Angle(degrees=ArcBetweenCentres)) ; print()
    # Skyfield arc between near limbs
    NearArc1 = Apparent1.separation_from(Apparent2).degrees - SemiDiameter(Apparent1).degrees - SemiDiameter(Apparent2).degrees
    # Trad arc
    NearArc2 = ArcBetweenCentres - SemiDiameter(Apparent1).degrees - SemiDiameter(Apparent2).degrees
    print('Skyfield:    Arc between near limbs: ', skyfield.units.Angle(degrees=NearArc1))
    print('Trad method: Arc between near limbs: ', skyfield.units.Angle(degrees=NearArc2)) ; print()
# ADDED CODE ENDS

#### Other Functions ####
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
        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()))

'''
# unused code
def PrintBodyInfo(Astrometric, Apparent): # print lots of debug info
    #print(f'## Astrometric : {Astrometric}')
    #print(f'## Apparent : {Apparent}')
    #print(f'## Apparent.center : {Apparent.center}')
    #print(f'## Apparent.target : {Apparent.target}')
    #print(f'## Ephemeris.names()[Apparent.center] : {Ephemeris.names()[Apparent.center]}') # doesn't work w/ Ephemeris + latlon
    #print(f'## Ephemeris.names()[Apparent.target] : {Ephemeris.names()[Apparent.target]}') # doesn't work w/ stars
    #print(f'# Astrometric.radec("date") : {Astrometric.radec("date")}')
'''

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

'''
Scratch pad for copy/paste command lines
python BR_Mods.py -t 2022-03-05T03:00:00 -n1 Aldebaran -n2 Moon -g 38.0,-122,0
# python BR_Mods.py
# python BR_Mods.py -t 2022-03-05T03:00:00 -n1 Markab -n2 Regulus -g 38.0,-122,1000
'''