search

skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.44k stars 214 forks source link

Calculating apparent magnitude for comets and minor planets #416

Closed Bernmeister closed 4 years ago

Bernmeister commented 4 years ago

(have edited this ticket to enhance the test program, update results and give a better explanation)

I have created a test program (see the end) to display the magnitude for a comet / minor planet both in PyEphem and Skyfield. Magnitude can be either absolute or apparent. I would like the apparent magnitude as, I believe, that is what one uses to determine whether an object is within viewing distance.

Running the test program yields:

PyEphem: 3.7.6.0
Skyfield: 1.33

PyEphem 88P/Howell 
    Magnitude (from PyEphem): 14.89 
    Earth Body Distance AU: 1.76101815700531 
    Apparent Magnitude (calculated): 14.89285058374846 

Skyfield 88P/Howell 
    Absolute Magnitude (H from data file): 11.0 
    Earth Body Distance AU: 1.761017531999252 
    Apparent Magnitude (calculated): 11.987619922504665 

PyEphem 1 Ceres 
    Magnitude (from PyEphem): 8.97 
    Absolute Magnitude (H from data file): 3.3399999141693115 
    Earth Body Distance AU: 2.8166115283966064 
    Apparent Magnitude (calculated): 8.967859508070683 

Skyfield (1) Ceres 
    Absolute Magnitude (H from data file): 3.34 
    Earth Body Distance AU: 2.8166120923794224 
    Apparent Magnitude (calculated): 8.967860459749872

PyEphem uses different data formats for comets versus minor planets. Looking at each of the data formats, it appears the gk model is used for comets, whereas the HG model is used for minor planets. Both models are described here. I have have taken the formulae specified and put together two functions, one for each model (see test code below). Using Skyfield to load comet / minor planet data and compute their position, I then use my two functions to compute apparent magnitude. Using PyEphem, each object loaded, computed and the magnitude determined, but I also calculate my own absolute magnitude using the data.

Looking at the results:

In short, there is a discrepancy in the apparent magnitude for the comet between PyEphem and Skyfield, but the numbers match exactly (as expected) for the minor planet apparent magnitude between PyEphem and Skyfield.

Would like some feedback please if the calculation and logic is sound and I am more than happy for the code to be reused if anyone finds it useful.

The test program, adjust for the current date and ensure you obtain the latest data for the comet and minor planet. See the comments in the code.

#!/usr/bin/env python3

# Refer to https://github.com/skyfielders/python-skyfield/issues/416

# Download and save 
#     https://www.minorplanetcenter.net/iau/Ephemerides/Comets/Soft00Cmt.txt
#     https://www.minorplanetcenter.net/iau/Ephemerides/Comets/Soft03Cmt.txt
#     https://www.minorplanetcenter.net/iau/Ephemerides/Bright/2018/Soft00Bright.txt
#     https://www.minorplanetcenter.net/iau/Ephemerides/Bright/2018/Soft03Bright.txt

# External verification:
# https://in-the-sky.org/data/object.php?id=0088P
# https://theskylive.com/88p-info
# https://heavens-above.com/comet.aspx?cid=88P
# https://www.calsky.com/cs.cgi/Comets/3?

import datetime, ephem, io, math, skyfield.api, skyfield.constants, skyfield.data.mpc

def pyephemCometMinorPlanet( utcNow, latitude, longitude, name, data, isComet ):
    observer = ephem.Observer()
    observer.date = ephem.Date( utcNow )
    observer.lat = str( latitude )
    observer.lon = str( longitude )

    body = ephem.readdb( getData( name, data ) )
    body.compute( observer )

    sun = ephem.Sun()
    sun.compute( observer )

    if isComet:
        apparentMagnitude = getApparentMagnitude_gk( body._g, body._k, body.earth_distance, body.sun_distance )
        print( "PyEphem", name,
               "\n\tMagnitude (from PyEphem):", body.mag,
               "\n\tEarth Body Distance AU:", body.earth_distance,
               "\n\tApparent Magnitude (calculated):", apparentMagnitude,
               "\n" )

    else:
        apparentMagnitude = getApparentMagnitude_HG( body._H, body._G, body.earth_distance, body.sun_distance, sun.earth_distance )
        print( "PyEphem", name,
               "\n\tMagnitude (from PyEphem):", body.mag,
               "\n\tAbsolute Magnitude (H from data file):", body._H,
               "\n\tEarth Body Distance AU:", body.earth_distance,
               "\n\tApparent Magnitude (calculated):", apparentMagnitude,
               "\n" )

def skyfieldCometMinorPlanet( utcNow, latitude, longitude, name, data, isComet ):
    timeScale = skyfield.api.load.timescale( builtin = True )
    topos = skyfield.api.Topos( latitude_degrees = latitude, longitude_degrees = longitude )
    ephemeris = skyfield.api.load( "de421.bsp" )
    sun = ephemeris[ "sun" ]
    earth = ephemeris[ "earth" ]

    with io.BytesIO( getData( name, data ).encode() ) as f:
        if isComet:
            dataframe = skyfield.data.mpc.load_comets_dataframe( f ).set_index( "designation", drop = False )
            body = sun + skyfield.data.mpc.comet_orbit( dataframe.loc[ name ], timeScale, skyfield.constants.GM_SUN_Pitjeva_2005_km3_s2 )

        else:
            dataframe = skyfield.data.mpc.load_mpcorb_dataframe( f ).set_index( "designation", drop = False )
            body = sun + skyfield.data.mpc.mpcorb_orbit( dataframe.loc[ name ], timeScale, skyfield.constants.GM_SUN_Pitjeva_2005_km3_s2 )

    t = timeScale.utc( utcNow.year, utcNow.month, utcNow.day, utcNow.hour, utcNow.minute, utcNow.second )
    alt, az, earthSunDistance = ( earth + topos ).at( t ).observe( sun ).apparent().altaz()
    ra, dec, sunBodyDistance = ( sun ).at( t ).observe( body ).radec()
    alt, az, earthBodyDistance = ( earth + topos ).at( t ).observe( body ).apparent().altaz()

    apparentMagnitude = getApparentMagnitude_HG( dataframe.loc[ name ][ "magnitude_H" ], dataframe.loc[ name ][ "magnitude_G" ], earthBodyDistance.au, sunBodyDistance.au, earthSunDistance.au )

    print( "Skyfield", name,
           "\n\tAbsolute Magnitude (H from data file):", dataframe.loc[ name ][ "magnitude_H" ],
           "\n\tEarth Body Distance AU:", earthBodyDistance.au,
           "\n\tApparent Magnitude (calculated):", apparentMagnitude,
           "\n" )

# https://stackoverflow.com/a/30197797/2156453
def getData( orbitalElementName, orbitalElementData ):
    return next( ( s for s in orbitalElementData if orbitalElementName in s ), None )

# https://www.clearskyinstitute.com/xephem/help/xephem.html#mozTocId564354
def getApparentMagnitude_gk( g_absoluteMagnitude, k_luminosityIndex, bodyEarthDistanceAU, bodySunDistanceAU ):
    return g_absoluteMagnitude + \
           5 * math.log10( bodyEarthDistanceAU ) + \
           2.5 * k_luminosityIndex * math.log10( bodySunDistanceAU )

# Calculate apparent magnitude (returns None on error).
#
# https://www.clearskyinstitute.com/xephem/help/xephem.html#mozTocId564354
# https://www.britastro.org/asteroids/dymock4.pdf
def getApparentMagnitude_HG( H_absoluteMagnitude, G_slope, bodyEarthDistanceAU, bodySunDistanceAU, earthSunDistanceAU ):
    beta = math.acos( \
                        ( bodySunDistanceAU * bodySunDistanceAU + \
                          bodyEarthDistanceAU * bodyEarthDistanceAU - \
                          earthSunDistanceAU * earthSunDistanceAU ) / \
                        ( 2 * bodySunDistanceAU * bodyEarthDistanceAU ) \
                    )

    psi_t = math.exp( math.log( math.tan( beta / 2.0 ) ) * 0.63 )
    Psi_1 = math.exp( -3.33 * psi_t )
    psi_t = math.exp( math.log( math.tan( beta / 2.0 ) ) * 1.22 )
    Psi_2 = math.exp( -1.87 * psi_t )

    # Have found a combination of G_slope, Psi_1 and Psi_2 can lead to a negative value in the log calculation.
    try:
        apparentMagnitude = H_absoluteMagnitude + \
                            5.0 * math.log10( bodySunDistanceAU * bodyEarthDistanceAU ) - \
                            2.5 * math.log10( ( 1 - G_slope ) * Psi_1 + G_slope * Psi_2 )

    except:
        apparentMagnitude = None

    return apparentMagnitude

utcNow = datetime.datetime.strptime( "2020-11-26", "%Y-%m-%d" )
latitude = -33
longitude = 151

print( "PyEphem:", ephem.__version__ )
print( "Skyfield:", skyfield.__version__ )
print()

with open( "Soft03Cmt.txt" ) as f:
    orbitalElementData = f.readlines()
    pyephemCometMinorPlanet( utcNow, latitude, longitude, "88P/Howell", orbitalElementData, True )

with open( "Soft00Cmt.txt" ) as f:
    orbitalElementData = f.readlines()
    skyfieldCometMinorPlanet( utcNow, latitude, longitude, "88P/Howell", orbitalElementData, True )

with open( "Soft03Bright.txt" ) as f:
    orbitalElementData = f.readlines()
    pyephemCometMinorPlanet( utcNow, latitude, longitude, "1 Ceres", orbitalElementData, False )

with open( "Soft00Bright.txt" ) as f:
    orbitalElementData = f.readlines()
    skyfieldCometMinorPlanet( utcNow, latitude, longitude, "(1) Ceres", orbitalElementData, False )
brandon-rhodes commented 4 years ago

Thanks for providing this sample code! I will probably not have time to look today, but likely will have time to try it out tomorrow and see about the difference in results between this code and the other library.

Bernmeister commented 4 years ago

I contact Elwood Downey of XEphem fame and described the issue of differences between values for apparent magnitude of comets. His response:

But two comments come to mind. One is that comets can change rapidly so the different brightness models may be updated often and separately. The other is there are in fact other definitions of H and G than the one I chose to implement at the time, so it is possible an alternate algorithm has come into current vogue.

So maybe my implementation is flawed from the outset...might have to go back to looking for the latest model(s)!

brandon-rhodes commented 4 years ago

@Bernmeister — Looking at this with fresh eyes after a few months, here’s my guess:

  1. The parameters 11.0 6.0 in Soft00Cmt.txt are g and k, not H and G. You need to pass them to your getApparentMagnitude_gk() routine.
  2. That would also agree with the predicted gray curve at http://astro.vanbuitenen.nl/comet/88 which is working from the same MPC data and predicts a magnitude of around 14.5 or 15.0 for November.
  3. As also shown in that graph, the MPC parameters happen to be an out of date prediction. The informative blurb down at the bottom of https://in-the-sky.org/data/object.php?id=0088P flatly recommends avoiding MPC magnitude curves if better data is available, I note.
  4. Note that either g/k or HG would fit the descriptions at https://minorplanetcenter.net/iau/info/CometOrbitFormat.html describing the parameters as Absolute magnitude and Slope parameter, so I think we have to use other data to decide which model is used by each file? You might email the MPC and ask that those descriptions be updated to name the parameters to clarify the model they are supposed to power.
Bernmeister commented 4 years ago

@brandon-rhodes Very confused now...I looked at _COMET_COLUMNS in mpc.py and specified are not g and k, but H and G:

    ('magnitude_H', (91, 95)),
    ('magnitude_G', (96, 100)),

Also I checked the MPC website and all I could find in regards to format for comets was here which is also present in mpc.py. This suggests the format is H, G, not g, k.

Yes, might drop the MPC a note to dig a little deeper...

Bernmeister commented 4 years ago

Did a debug of the test code above. In the function pyephemCometMinorPlanet() when processing Soft03Cmt.txt for 88P/Howell, I added the line ra = body.a_ra at the end of the function (applied knowledge: now that I know compute() is lazy!).

If I inspect the comet/body I see that _g = 11.0, _k = 6.0 and for _G and _H:

str: Traceback (most recent call last): File "/home/bernard/.eclipse/360744286_linux_gtk_x86_64/plugins/org.python.pydev.core_7.5.0.202001101138/pysrc/_pydevd_bundle/pydevd_resolver.py", line 214, in _get_py_dictionary attr = getattr(var, name) RuntimeError: this object has g/k magnitude coefficients

Still digging and more confused!

brandon-rhodes commented 4 years ago

Very confused now...I looked at _COMET_COLUMNS in mpc.py…

Those were simply my own guesses, which I am now suspecting were, alas, incorrect.

Also I checked the MPC website and all I could find in regards to format for comets was here which is also present in mpc.py. This suggests the format is H, G, not g, k.

What in that text suggests HG? I don't myself see any mention of HG.

If I inspect the comet/body I see that _g = 11.0, _k = 6.0 and for _G and _H … Still digging and more confused!

Could you elaborate on the source of confusion? Thanks to Elwood’s excellent data file design, the data file explicitly declares that the magnitude parameters are "g/k" parameters. Look for the little letter "g" right before the 11 and 6:

88P/Howell,e,4.3836,56.6820,235.9072,3.105692,0.1800804,0.56432755,11.0554,11/27.0/2020,2000,g 11.0,6.0

The data format is documented here, including the fact that the g specifies a g/k magnitude model:

http://www.clearskyinstitute.com/xephem/help/xephem.html#mozTocId564354

I think the only problem here is that Soft00Cmt.txt required me to guess at the meaning of its magnitude parameters, and without thinking to check them against the same data file in PyEphem format (which I didn't know they offered!), I guessed wrong about which model they were using.

Bernmeister commented 4 years ago

@brandon-rhodes Think I've nailed it...

I currently use the following data files from the MPC within PyEphem (one for comets and four for various minor planets):

https://minorplanetcenter.net/iau/Ephemerides/Comets/Soft03Cmt.txt
https://minorplanetcenter.net/iau/Ephemerides/Bright/2018/Soft03Bright.txt
https://minorplanetcenter.net/iau/Ephemerides/CritList/Soft03CritList.txt
https://minorplanetcenter.net/iau/Ephemerides/Distant/Soft03Distant.txt
https://minorplanetcenter.net/iau/Ephemerides/Unusual/Soft03Unusual.txt

Ditto for Skyfield:

https://minorplanetcenter.net/iau/Ephemerides/Comets/Soft00Cmt.txt
https://minorplanetcenter.net/iau/Ephemerides/Bright/2018/Soft00Bright.txt
https://minorplanetcenter.net/iau/Ephemerides/CritList/Soft00CritList.txt
https://minorplanetcenter.net/iau/Ephemerides/Distant/Soft00Distant.txt
https://minorplanetcenter.net/iau/Ephemerides/Unusual/Soft00Unusual.txt

For PyEphem, in section 7.1.2.2 there is Field 12 (for when Field 2 is 'e'). This field either starts with 'g' or 'H' representing the model to be used.

If we take C/1995 O1 (Hale-Bopp) from Soft03Cmt.txt, note the 'g' and values of -2.0 and 4.0 (g and k). Comparing the same body in Soft00Cmt.txt, there are also the same values -2.0 and 4.0, hence comet data in the MPC format is g, k.

In mpc.py you define _COMET_COLUMNS with the fields magnitude_H and magnitude_G but really should be magnitude_g and magnitude_k.

As for the minor planet data files, things check out, so no issues there.

I have modified the test code to now account for misunderstanding that comet data for Skyfield is in H, G format, whereas it should have been g, k. Code is below...

#!/usr/bin/env python3

# Refer to https://github.com/skyfielders/python-skyfield/issues/416

# Download and save 
#     https://www.minorplanetcenter.net/iau/Ephemerides/Comets/Soft00Cmt.txt
#     https://www.minorplanetcenter.net/iau/Ephemerides/Comets/Soft03Cmt.txt
#     https://www.minorplanetcenter.net/iau/Ephemerides/Bright/2018/Soft00Bright.txt
#     https://www.minorplanetcenter.net/iau/Ephemerides/Bright/2018/Soft03Bright.txt

# External verification:
# https://in-the-sky.org/data/object.php?id=0088P
# https://theskylive.com/88p-info
# https://heavens-above.com/comet.aspx?cid=88P
# https://www.calsky.com/cs.cgi/Comets/3?

import datetime, ephem, io, math, skyfield.api, skyfield.constants, skyfield.data.mpc

def pyephemCometMinorPlanet( utcNow, latitude, longitude, name, data, isComet ):
    observer = ephem.Observer()
    observer.date = ephem.Date( utcNow )
    observer.lat = str( latitude )
    observer.lon = str( longitude )

    body = ephem.readdb( getData( name, data ) )
    body.compute( observer )

    sun = ephem.Sun()
    sun.compute( observer )

    if isComet:
        apparentMagnitude = getApparentMagnitude_gk( body._g, body._k, body.earth_distance, body.sun_distance )
        print( "PyEphem", name,
               "\n\tMagnitude (from PyEphem):", body.mag,
               "\n\tAbsolute Magnitude (g from data file):", body._g,
               "\n\tEarth Body Distance AU:", body.earth_distance,
               "\n\tApparent Magnitude (calculated):", apparentMagnitude,
               "\n" )

    else:
        apparentMagnitude = getApparentMagnitude_HG( body._H, body._G, body.earth_distance, body.sun_distance, sun.earth_distance )
        print( "PyEphem", name,
               "\n\tMagnitude (from PyEphem):", body.mag,
               "\n\tAbsolute Magnitude (H from data file):", body._H,
               "\n\tEarth Body Distance AU:", body.earth_distance,
               "\n\tApparent Magnitude (calculated):", apparentMagnitude,
               "\n" )

def skyfieldCometMinorPlanet( utcNow, latitude, longitude, name, data, isComet ):
    topos = skyfield.api.Topos( latitude_degrees = latitude, longitude_degrees = longitude )
    ephemeris = skyfield.api.load( "de421.bsp" )
    sun = ephemeris[ "sun" ]
    earth = ephemeris[ "earth" ]

    timeScale = skyfield.api.load.timescale( builtin = True )
    t = timeScale.utc( utcNow.year, utcNow.month, utcNow.day, utcNow.hour, utcNow.minute, utcNow.second )
    alt, az, earthSunDistance = ( earth + topos ).at( t ).observe( sun ).apparent().altaz()

    with io.BytesIO( getData( name, data ).encode() ) as f:
        if isComet:
            dataframe = skyfield.data.mpc.load_comets_dataframe( f ).set_index( "designation", drop = False )
            body = sun + skyfield.data.mpc.comet_orbit( dataframe.loc[ name ], timeScale, skyfield.constants.GM_SUN_Pitjeva_2005_km3_s2 )

        else:
            dataframe = skyfield.data.mpc.load_mpcorb_dataframe( f ).set_index( "designation", drop = False )
            body = sun + skyfield.data.mpc.mpcorb_orbit( dataframe.loc[ name ], timeScale, skyfield.constants.GM_SUN_Pitjeva_2005_km3_s2 )

    ra, dec, sunBodyDistance = ( sun ).at( t ).observe( body ).radec()
    alt, az, earthBodyDistance = ( earth + topos ).at( t ).observe( body ).apparent().altaz()

    if isComet:
        # The fields "magnitude_H" and "magnitude_G" really should be called "magnitude_g" and "magnitude_k".
        apparentMagnitude = getApparentMagnitude_gk( dataframe.loc[ name ][ "magnitude_H" ], dataframe.loc[ name ][ "magnitude_G" ], earthBodyDistance.au, sunBodyDistance.au )
        print( "Skyfield", name,
               "\n\tAbsolute Magnitude (g from data file):", dataframe.loc[ name ][ "magnitude_H" ],
               "\n\tEarth Body Distance AU:", earthBodyDistance.au,
               "\n\tApparent Magnitude (calculated):", apparentMagnitude,
               "\n" )

    else:
        apparentMagnitude = getApparentMagnitude_HG( dataframe.loc[ name ][ "magnitude_H" ], dataframe.loc[ name ][ "magnitude_G" ], earthBodyDistance.au, sunBodyDistance.au, earthSunDistance.au )

        print( "Skyfield", name,
               "\n\tAbsolute Magnitude (H from data file):", dataframe.loc[ name ][ "magnitude_H" ],
               "\n\tEarth Body Distance AU:", earthBodyDistance.au,
               "\n\tApparent Magnitude (calculated):", apparentMagnitude,
               "\n" )

# https://stackoverflow.com/a/30197797/2156453
def getData( orbitalElementName, orbitalElementData ):
    return next( ( s for s in orbitalElementData if orbitalElementName in s ), None )

# https://www.clearskyinstitute.com/xephem/help/xephem.html#mozTocId564354
def getApparentMagnitude_gk( g_absoluteMagnitude, k_luminosityIndex, bodyEarthDistanceAU, bodySunDistanceAU ):
    return g_absoluteMagnitude + \
           5 * math.log10( bodyEarthDistanceAU ) + \
           2.5 * k_luminosityIndex * math.log10( bodySunDistanceAU )

# Calculate apparent magnitude (returns None on error).
#
# https://www.clearskyinstitute.com/xephem/help/xephem.html#mozTocId564354
# https://www.britastro.org/asteroids/dymock4.pdf
def getApparentMagnitude_HG( H_absoluteMagnitude, G_slope, bodyEarthDistanceAU, bodySunDistanceAU, earthSunDistanceAU ):
    beta = math.acos( \
                        ( bodySunDistanceAU * bodySunDistanceAU + \
                          bodyEarthDistanceAU * bodyEarthDistanceAU - \
                          earthSunDistanceAU * earthSunDistanceAU ) / \
                        ( 2 * bodySunDistanceAU * bodyEarthDistanceAU ) \
                    )

    psi_t = math.exp( math.log( math.tan( beta / 2.0 ) ) * 0.63 )
    Psi_1 = math.exp( -3.33 * psi_t )
    psi_t = math.exp( math.log( math.tan( beta / 2.0 ) ) * 1.22 )
    Psi_2 = math.exp( -1.87 * psi_t )

    # Have found a combination of G_slope, Psi_1 and Psi_2 can lead to a negative value in the log calculation.
    try:
        apparentMagnitude = H_absoluteMagnitude + \
                            5.0 * math.log10( bodySunDistanceAU * bodyEarthDistanceAU ) - \
                            2.5 * math.log10( ( 1 - G_slope ) * Psi_1 + G_slope * Psi_2 )

    except:
        apparentMagnitude = None

    return apparentMagnitude

utcNow = datetime.datetime.strptime( "2020-11-29", "%Y-%m-%d" )
latitude = -33
longitude = 151

print( "PyEphem:", ephem.__version__ )
print( "Skyfield:", skyfield.__version__ )
print()

with open( "Soft03Cmt.txt" ) as f:
    orbitalElementData = f.readlines()
    pyephemCometMinorPlanet( utcNow, latitude, longitude, "88P/Howell", orbitalElementData, True )

with open( "Soft00Cmt.txt" ) as f:
    orbitalElementData = f.readlines()
    skyfieldCometMinorPlanet( utcNow, latitude, longitude, "88P/Howell", orbitalElementData, True )

with open( "Soft03Bright.txt" ) as f:
    orbitalElementData = f.readlines()
    pyephemCometMinorPlanet( utcNow, latitude, longitude, "1 Ceres", orbitalElementData, False )

with open( "Soft00Bright.txt" ) as f:
    orbitalElementData = f.readlines()
    skyfieldCometMinorPlanet( utcNow, latitude, longitude, "(1) Ceres", orbitalElementData, False )
brandon-rhodes commented 4 years ago

I am glad we are now in agreement — and that your magnitude predictions come out correctly when the data file is interpreted as containing g/k parameters. I have just landed a commit (see the link right above this comment) to rename those two columns, so those should change with the next release of Skyfield (which will probably come within the next two weeks). I am trying to think of a way to also support code using the old names but I'm not sure there are many good options, and hopefully the columns are not widely used since they are not yet documented.

But if I think of a good way to support legacy code, I'll follow up with a further commit.

Bernmeister commented 2 years ago
  1. Note that either g/k or HG would fit the descriptions at https://minorplanetcenter.net/iau/info/CometOrbitFormat.html describing the parameters as Absolute magnitude and Slope parameter, so I think we have to use other data to decide which model is used by each file? You might email the MPC and ask that those descriptions be updated to name the parameters to clarify the model they are supposed to power.

As a result of lots of testing and revisiting, partly initiated by my discussion, I contacted the MPC to clarify the magnitude parameters for comets and minor planets:

I'd like to confirm please the naming of the magnitude parameters for the formats for both Minor Planets and Comets. My question arises as I have been using XEphem/PyEphem which switches between the g,k magnitude model or the H,G magnitude model, depending on the format of the XEphem specific data files, https://www.minorplanetcenter.net/iau/Ephemerides/Soft03.html. I now plan on moving to the MPC formatted data files and wanted to clarify the parameter names.

For Minor Planets, https://www.minorplanetcenter.net/iau/info/MPOrbitFormat.html, there are defined fields for absolute magnitude (H) at columns 9 - 13 and slope parameter (G) at columns 15 - 19. As a matter of redundancy, I assume these columns are indeed intended for the H,G magnitude model as named?

For Comets, https://www.minorplanetcenter.net/iau/info/CometOrbitFormat.html, specifically the second section entitled "Ephemerides and Orbital Elements Format", there are defined fields for absolute magnitude at columns 92 - 95 and slope parameter at columns 97 - 100. As there is no mention of 'H' or 'G' on the column use, are these columns specifically only for the H,G magnitude model?

MPC's response:

For minor planets the parameters do correspond to H and G. For comets they correspond to m0 and (2.5n) in the formula: m = m0 + 2.5nlog(r) + 5log(delta) where m is the apparent magnitude, m0 the absolute magnitude, r the distance from the Sun, and delta the distance from the Earth.

For minor planets, the H and G parameters match with _MPCORB_COLUMNS in mpc.py.

For comets, there is a match for g in _COMET_COLUMNS but it is ambiguous to me if (2.5*n) above corresponds to k or 2.5 * k.

I have gone down a rabbit hole (several actually) in calculating apparent magnitude for comets / minor planets, using both PyEphem and Skyfield. When comparing various bodies between the backends and other external sources, I get mismatches. I suspect that the magnitude parameters supplied by the MPC are somewhat dated and therefore insufficient to calculate a current apparent magnitude (but still perfectly fine for calculating rise/set/az/alt/ra/dec). Several websites use recent observations to determine the apparent magnitude and I have found, at first glance at least, these websites match with each other (could be they are all using the same data source, I suppose).

In short, @brandon-rhodes , the fields you have defined for comets might need to be looked at (is it k or 2.5 * k), but for now I'd leave things as they are. Given that Skyfield doesn't calculate apparent magnitude for comets / minor plants (yet), the onus is on the end user to beware...

jurezakrajsek commented 2 years ago

@Bernmeister I use this function to calculate the magnitudes of comets for all calculation of magnitude at cobs.si and it seems to calculate magnitudes accurately so the k is not (2.5k) but just k

def obj_magnitude(self):
    r = self.comet_obj.distance()
    ra, dec, delta = self.observed_comet.radec()
    mag_g = self.row.magnitude_g
    mag_k = self.row.magnitude_k

    if self.num > 1:
        return np.asarray(
            [mag_g + 5 * math.log10(di) + 2.5 * mag_k * math.log10(ri) for ri, di in zip(r.au, delta.au)])
    else:
        return mag_g + 5 * math.log10(delta.au) + 2.5 * mag_k * math.log10(r.au)

The g an k parameters in the MPC files are usually way off and the predictions out of those parameters are not the best. Since I have lot of comet observations in the cobs database I do a fit of g and k from the observations and use those parameters for the magnitude estimation.

You can also download the MPC file (in 1-line or ephem format) from cobs.si, but i would need to check if the g and k parameters in the output file are beeing substituted for the ones calculated by cobs.

Clear skies