How do I calculate sunrise/sunset times for non-Earth planets?

[snoopy]

EDIT: I was right after all!

I'm probably being hasty, posting before I test my results, but I think I can solve for the hour angle now, which I believe I can use to get time using the method from my previous post.

Here we go:

1. Get the Mean Anomaly. You can do this by simply picking a value between 0 and 360, but it is the time since the planet last passed perihelion moving at its mean motion. In math terms:

M = sqrt((G/(m1+m2))/a^3)*t
where
M is the mean anomaly
G is the gravitational constant in m^3 kg^-1 s^-2 (which you could define in terms of your custom time scale)
m1 and m2 are the masses of your orbiting bodies in kg
a is the semi-major axis of your planets orbit in meters
t is the time in seconds (don't forget to convert to custom time if you did that for the gravitational constant as well)

2. Use the mean anomaly to calculate the equation of center, which is the difference between true anomaly (actual angle along the ecliptic) and mean anomaly (which treats the ecliptic as a circle). This is done with the (accurate to 4 decimal places) approximation

C = (2e-0.25e^3)*sinM + 1.25e^2*sin2M + (13/12)e^3*sin3M
where
C is the equation of center
e is the eccentricity of the orbit (which you will need to decide on yourself)
M is the mean anomaly (I think it needs to be in degrees, even though it would appear to output an answer in radians. This seems wrong, but it gave me the correct answer when I was following an explanation. This will definitely be something I double check)

3. Find the ecliptic longitude using

λ = M + C + Π + 180
where
λ is the ecliptical longitude
Π is the longitude of perihelion (pericenter, periapsis) in degrees (you obtain this by adding the Argument of Perihelion and the Longitude of the Ascending Node, both of which you will have from defining your orbit. If you do it in Space Engine, remember that the values in that are relative to our sun. You will need to get the difference in those values between your sun and your planet)
180 because it is the sun viewed from the planet, not the planet viewed from the sun

4. Find the declination of the sun using

sinδ = sinλ + sinA
where
δ is the declination
A is the axial tilt in degrees

5. Calculate the hour angle using

cosω = (sin(-0.83)-sin(Phi)*sinδ) / (cos(Phi)*cosδ)
where
ω is the hour angle
-0.83 degrees because of atmospheric refraction (assumes Earth-like atmosphere)
Phi is the latitude of the observer

And there it is, generalized for other planets. Correctly, I hope.

I found this link very helpful: http://aa.quae.nl/en/reken/zonpositie.html