And the Evening and the Morning Were the First Day

This is a topic that is near and dear to my heart, because it goes so far beyond the bounds of nerdery that it boggles my own mind to imagine that I spent my time on this.

When does the sun rise? When does it set?

This changes depending on where you are on the Earth, or if you're on a mountain. So, what time? Don't cheat.

Start with these variables (assuming spherical planet, normal orbital, massively dumbed down):
\[c_E, r_E, n, h, y, a, G, e, \Pi, m_S, m_E, \varepsilon, \phi, \theta\]
Where $c_E$ is the circumference of your world, $r_E$ is the radius, $n$ is the number of days since the beginning of the year, $h$ is the observer's height above sea level in feet, $y$ is the number of days in the year, $a$ is the distance between the world and the star, $G$ is the gravitational constant, $e$ is the orbital eccentricity, $\Pi$ is the argument of perihelion, $m_S$ is the mass of the sun, $m_E$ is the mass of the world, $\varepsilon$ is the world's axial tilt, $\phi$ is the observer's latitude, and $\theta$ is the observer's longitude.

Got all that? Now for the hard parts.

Find the standard graviational parameter $\mu = G \cdot (m_S + m_E)$ of the system. The mean anomaly $M$ is then found by:
\[M = \sqrt{\mu \over a^3} {n \over y} {24 \cdot 60 \cdot 60 \cdot \textrm{ years} \over \textrm{seconds}}\]
To find the true anomaly $\nu$, we must use first-kind Bessel functions $J_\alpha$ and $\beta^m$:
\[\beta^m = \left({e \over 2}\right)^m \left(1 + m \sum_{t=1}^\infty {(2t + m - 1)! \over t! (t + m)!} \left({e \over 2}\right)^{2t}\right)\]
\[\nu = M + 2 \sum_{s=1}^\infty\left({1 \over s}\left(J_s(s \cdot e) + \sum_{p=1}^\infty\left(\beta(p) (J_{s-p}(s \cdot e) + J_{s+p}(s \cdot e))\right)\right)\sin(s \cdot M)\right)\]
Make sense? Great! That's the hard part. I don't recommend doing it by hand.
With those solved, we can find the equation of the center $C$.
\[C = \nu - M\]
And then the ecliptic longitude of the sun $\lambda$.
\[\lambda = (M + C + \Pi + 180^\circ) \textrm{ mod } 360^\circ\]
Then the declination of the sun $\delta$.
\[\delta = \sin^{-1}\left(\sin(\lambda)\cdot\sin(\varepsilon)\right)\]
The height correction $h_{corr}$ is given by:
\[h_{corr} = -1.15^\circ \sqrt{h\over\textrm{ft}} {1 \over 60}\]
Then the hour angle $\omega$ is given by:
\[\omega = \cos^{-1}\left({\sin(-0.83^\circ + h_{corr}) - \sin\phi\sin\delta\over \cos\phi\cos\delta}\right)\]
Almost there. Now for the equation of time, which sounds pretty high fantasy to me.
\[\Delta t_{ey} = -2 e \sin M + \tan^2{\varepsilon\over 2} \sin(2 M + 2 \lambda)\]
I also add a correction for the time zone, but this will depend on whether your world has time zones. Otherwise, the calcuation will yield sunrise/set relative to your equivalent of Greenwich Observatory.
\[TZ = -{\theta \over 15^\circ}\]
The final bit of information is the calculation of the Julian transit (or solar noon). This is when the sun is directly overhead.
\[J_t = 0.5 + {\theta \over 360^\circ} + \Delta t_{ey} + TZ\]
Basically done. The hour angle $\omega$ is the difference between the sunset/rise and $J_t$!. $J_{set}$ and $J_{rise}$ will be the hour of set and rise. Minute and seconds can be derived trivially.
\[J_{set} = \left(J_t + {\omega\over360^\circ}\right) 24\textrm{ hr}\]
\[J_{rise} = \left(J_t - {\omega\over360^\circ}\right) 24\textrm{ hr}\]

The last step?

Draw your analemma (a plot of $\delta - \phi$ vs $\Delta t_{ey}$).