I believe what you'll want to do is create an SPK kernel for the ephemeris of the point on the surface of the moon and then whatever ID you give that point should be the center of the topocentric frame. That is my understanding from this info:
"This method of definition is particularly well suited for defining topocentric frames on the surface of the Earth. For example, suppose you have an SPK (ephemeris) file that specifies the location of some surface point on the Earth, and that the SPK ID code of this point is 399100. Moreover suppose you have the geodetic co-latitude (COLAT) and longitude (LONG) measured in degrees for this point. (Note that the co-latitude is the complement of latitude: latitude + co-latitude = 90 degrees.)" from: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/frames.html
To create the SPK kernel, you can get the inertial position vector by rotating the vector defined by the lat/long coordinates into an inertial frame, and then doing finite differences to calculate the velocity. This is how I did it in Python (using the SPICE python wrapper called spiceypy):
'''
Create SPK kernel from Lunar lat-lon coordinates
'''
from sys import path
path.append( '/home/alfonso/AWP/python_tools' )
import plotting_tools as pt
import planetary_data as pd
import spice_tools as st
import numerical_tools as nt
import spiceypy as spice
import numpy as np
fn = '/mnt/c/Users/alfon/AWP/misc/'
r2d = np.pi / 180.0
if __name__ == '__main__':
spice.furnsh( '../../spice/lsk/naif0012.tls' )
spice.furnsh( '../../spice/pck/pck00010.tpc' )
lat = 50.0 * r2d
lon = 100.0 * r2d
r = 1737.1
et0 = spice.str2et( '2021-05-17' )
etf = spice.str2et( '2022-01-01' )
delta_t = 1000
ets = np.arange( et0, etf, delta_t )
states = np.zeros( ( len( ets ), 6 ) )
pos_iau_moon = spice.latrec( r, lon, lat )
for n in range( len( ets ) ):
states[ n, :3 ] = np.dot(
spice.pxform( 'IAU_MOON', 'J2000', ets[ n ] ),
pos_iau_moon )
for n in range( len( ets ) - 1 ):
states[ n, 3: ] = ( states[ n + 1, :3 ] - states[ n, :3 ] ) / delta_t
pos_np1 = np.dot(
spice.pxform( 'IAU_MOON', 'J2000', ets[ -1 ] + delta_t ),
pos_iau_moon )
states[ -1, 3: ] = ( pos_np1 - states[ -1, :3 ] ) / delta_t
handle = spice.spkopn( 'lunar_point.bsp', 'SPK_file', 0 )
point_id = 301100
center = 301
frame = 'J2000'
degree = 5
spice.spkw09( handle, point_id, center, frame,
ets[ 0 ], ets[ -1 ], '0', 5, len( ets ),
states.tolist(), ets.tolist() )
spice.spkcls( handle )
spice.furnsh( 'lunar_point.bsp' )
states_j2000 = st.calc_ephemeris( point_id, ets, 'J2000', 301 )
states_iau_moon = st.calc_ephemeris( point_id, ets, 'IAU_MOON', 301 )
latlons_j2000 = nt.inert2latlon( states_j2000[ :, :3 ], 'J2000', 'IAU_MOON', ets )
latlons_iau_moon = nt.bf2latlon( states_iau_moon[ :, :3 ] )
pt.plot_groundtracks(
[ latlons_j2000, latlons_iau_moon ],
{
'labels' : [ 'J2000', 'IAU_MOON' ],
'plot_coastlines': False,
'filename' : fn + 'lunar_point_groundtracks.png',
'dpi' : 300
} )
pt.plot_orbits( [ states_j2000[ :, :3 ], states_iau_moon[ :, :3 ] ],
{
'labels' : [ 'J2000', 'IAU_MOON' ],
'cb_radius': 1737.1,
'cb_cmap' : 'gist_yarg',
'azimuth' : 41.0,
'elevation': 30.0,
'axes_mag' : 1.0,
'filename' : fn + 'lunar_point_3d.png',
'dpi' : 300
} )
In this case I used IAU_MOON but you could just replace that with MOON_ME and the ephemeris times that you are interested in.
The imports at the top (plotting_tools, planetary_data, spice_tools, numerical_tools) are all personal libraries, but you don't need them for the main part of the program (writing the BSP kernel), only after to get some visual confirmation that the BSP was written correctly. But for reference, here is the calc_ephemeris function in spice_tools (just a convenience wrapper to spkezr and spkgeo):
def calc_ephemeris( target, ets, frame, observer ):
'''
Wrapper for spkezr and spkgeo
'''
if type( target ) == str:
return np.array( spice.spkezr( target, ets, frame, 'NONE', observer )[ 0 ] )
else:
n_states = len( ets )
states = np.zeros( ( n_states, 6 ) )
for n in range( n_states ):
states[ n ] = spice.spkgeo( target, ets[ n ], frame, observer )[ 0 ]
return states
I have part of the plotting_tools file posted here under src/python_tools: https://github.com/alfonsogonzalez/AWP
I show how to write a bunch of orbital mechanics software with Python in YouTube videos (if you just search "orbital mechanics with python" on YouTube you will find them, not sure if I'm allowed to post links).
And you're right, you could write this BSP kernel with the position w.r.t a body-fixed frame (IAU_MOON, MOON_ME), and I am just realizing that that would definitely be easier. At first, I did it in the inertial frame because calculating derivatives in the inertial frame are straight forward and you don't have to take into account the angular velocity of the frame since its inertial. But actually in this case the velocity of the point on the surface is just 0, so you're right I think you should do that instead. I just redid it that way and got what we would expect:
