Where is the sun relative to my lunar rover? I'm a robotics researcher writing orienteering code for a solar-powered lunar rover.
I have two similar problems:


Given the position of a robot on the moon in lunar latitude/longitude, where is the sun? That is, which direction should I orient my solar panel to maximize incidence?
For testing purposes, I need to solve the same problem on Earth. Given Earth lat/lon, where is the sun?


To be honest, I'm not even sure what coordinate system to use to express the answer to these questions. I was considering a coordinate frame centered on the rover with East, North, and Up axes, but I'm open to a better suggestion!

So far, I've found some information about VSOP87 and ELP2000, and I've started reading "Astronomical Algorithms", but I'd appreciate help from the experts.

Edit

Here's what I've got so far. I can get the coordinates of the Earth-Moon barycenter using VSOP2013. If I can get from the barycenter to the center of the Earth, I can probably use EPS2000/82 to find the position of the Moon in heliocentric coordinates.

Sun->VSOP->Earth-Moon Barycenter->Earth->EPS2000->Moon

At very least, I think this would get me the moon coordinates in the sun frame. I'm hoping it would also give me the orientation of the Moon in the sun frame?
If you like Python then see Skyfield; Observing from a Moon location which I believe will do exactly what you need! Also see for example this answer to Ephemeris for lunar body-centered body-fixed coordinates? written by Skyfield's author.

For example, you posted your question at 2020-04-20 15:21:36 UTC. Let's say your rover has just driven up to Surveyor III at 3° 0′ 58.03″ S, 23° 25′ 4.48″ W  (-3.01612°, -23.41791°) (read more about Surveyor III in How did Apollo-12 manage to land next to Surveyor-3? First “Space-Tourists”?)

At that moment the Sun was below the horizon. The output of the script below is

RA:  01h 54m 39.28s
Dec:  +11deg 45' 38.6"
alt:  -38deg 30' 34.6" degrees (below the horizon)
az:  265deg 50' 21.5" degrees around the horizon from north



  I can get the coordinates of the Earth-Moon barycenter using VSOP2013


factoids: VSOP is what PyEphem uses under the hood. It is currently in maintenance mode maintained by the same person who wrote Skyfield. 



#!/usr/bin/env python
# -*- coding: utf-8 -*-
# adapted from https://rhodesmill.org/skyfield/planetary.html

from skyfield.api import PlanetaryConstants, Loader

load = Loader('~/Documents/fishing/SkyData')  # avoids multiple copies of large files
ts = load.timescale(builtin=True)
t = ts.utc(2020, 4, 20, 15, 21, 36)

eph = load('de421.bsp')
earth, moon, sun = [eph[thing] for thing in ('earth', 'moon', 'sun')]

pc = PlanetaryConstants()
pc.read_text(load('moon_080317.tf'))
pc.read_text(load('pck00008.tpc'))
pc.read_binary(load('moon_pa_de421_1900-2050.bpc'))

frame = pc.build_frame_named('MOON_ME_DE421')
Surveyor_III = moon + pc.build_latlon_degrees(frame, latitude_degrees = -3.01612,
                                              longitude_degrees = -23.41791)

apparent = Surveyor_III.at(t).observe(sun).apparent()
ra, dec, distance = apparent.radec()
print('RA: ', ra)
print('Dec: ', dec)

alt, az, distance = apparent.altaz()
if alt.degrees &gt;= 0:
    print('alt: ', alt, 'degrees (above the horizon)')
else:
    print('alt: ', alt, 'degrees (below the horizon)')
print('az: ', az, 'degrees around the horizon from north')