Eventually I ran across a post discussing the merits of different simulation environments, and learned about GMAT, or the "General Mission Analysis Tool". This is a free, open-source package, developed by NASA, that has been used to design real missions. Since August of last year I have been learning the ins and outs of this environment, and using it to search for Snoopy. Unfortunately the acronym for this package is shared with another important "GMAT", the test that people take to get into graduate school. If you go searching the internet for this package, you may find you have better luck spelling out its name in full.
You can download GMAT here. There is a youtube channel with tutorials on how to use it here. Documentation for GMAT exists in various places, but I find this to be the easiest to use. As the name implies, it is a general tool, that can simulate objects in Earth orbit, Mars missions, etc.
One of the first things you will want to do with GMAT for any lunar simulation work is to upgrade the lunar gravity model. The gravity field of the moon is notoriously "lumpy", due to mass concentrations, or "mascons". The gravity model accounts for this lumpiness by decomposing the field using spherical harmonics, and there are a few things you need to know about it. The model that is installed by default is a decent model, named LP165P, but you should upgrade to a GRAIL model. GRAIL models derive from data gathered by a pair of lunar satellites, named Ebb and Flow, that tracked each other around the moon for several years. The resulting data has been processed into a series of ever-more-detailed gravity models. I have been using this one. In order to use it you download the file, save it as "jggrx_0420a_sha.tab", and copy it into the lunar gravity file directory in your GMAT installation. On my installation it is at ...\AppData\Local\GMAT\R2018a\data\gravity\luna\
|This image was derived from GRAIL data and shows the local "lumpy" variations in the moon's gravity field. The "lumps" are likely from heavy objects that hit the moon and then stayed close to it's surface, creating mass concentrations or "mascons"
The model above has "420" in the name because it can run with the degree and order set as high as 420. This is how many spherical harmonics can be used to simulate the field, and more is better in terms of the fidelity. However, more is also slower, in terms of simulation time. I once ran a simulation with a high-fidelity model, with the degree and order set to 1000. This simulation ran for about a week on my computer, simulating a month in orbit. This is too slow for most of what I do, and it does not seem to be required. By running the same simulation with higher degree/order, I have noticed that there are diminishing returns, i.e. very small differences in the results, above degree/order setting of 150/150.