Seasonal changes in hours of daylight I will post my own answer to this question unless someone else posts the same answer first, but I am curious to know what other points of view might lead to different ways of answering it.
Temporarily (for the duration of this question) assume for simplicity that the earth is perfectly spherical, that the earth's orbit is a perfect circle, and that the sun appears in the sky as a point rather than as a disk.
The earth has an "axial tilt" of about $23.56^\circ$.  (I'd like to see a graphic showing the ecliptic as tilted and the earth's axis as vertical.  That's an equally valid viewpoint and it might make certain things clearer for psychological reasons.)
How do we find the number of hours of daylight at a specified latitude at a specified date?
Appendix:
Two observations:


*

*The incorrect idea that the number of hours would be $12+f(\text{latitude}) \times \sin\Big( \text{calendar date} \in [0,2\pi)\Big)$ seems plausible if you don't think about it: circular motion; periodic oscillation, etc. . . . . . .  But it's wrong.

*A latitude halfway between the equator and the tropic of cancer has two dates when the sun's angular elevation above the horizon is $90^\circ$: one before and one after the June solstice.  I thought maybe that meant it would have two dates of longest hours of daylight.  Not so: like everywhere else south of the arctic circle and north of the equator, it has the longest hours of daylight at the June solstice.
 A: A point at latitude $\lambda$ goes around a circle of radius $\cos\lambda$ once every $24$ hours, in a plane at distance $\sin\lambda$ from the plane of the equator.
Suppose the earth is at a point when the angle between its axis and the ecliptic is $\tau\in(-23.56^\circ,+23.56^\circ)$.  At latitude $\lambda$ the sun rises $\gamma\cdot\dfrac{24}{2\pi}$ hours after solar midnight and sets $(2\pi-\gamma)\cdot\dfrac{24}{2\pi}$ hours after solar midnight.  Bisect the earth with a plane orthogonal to the line from the sun to the earth.  The sun rises at a point at distance $\sin\gamma\cos\lambda$ from that plane.  But if you look at a triangle in the plane through the sun and the north and south poles, you see that that distance is $\tan\tau\sin\lambda$.  Therefore
$$
\sin\gamma\cos\lambda = \tan\tau\sin\lambda,
$$
or
$$
\sin\gamma = \tan\tau\tan\lambda. \tag 1
$$
(The equality $(1)$ is valid between the arctic and antarctic circles; at more extreme latitudes it would involve a sine exceeding $1$ in absolute value.) Next observe that $\tan\tau=\tan23.56^\circ\sin\theta$, where $\theta\in[0,2\pi)$ is the calendar date, with $\theta=0$ at the equinox.  Hence
$$
\sin\gamma=\tan\lambda\tan23.56^\circ\sin\theta.
$$
Therefore
$$
\text{length of daylight} = (\text{24 hours})\cdot\left(\frac 1 2 + \frac{\arcsin(\tan\lambda\tan23.56^\circ\sin\theta)}{\pi}\right).
$$
This approaches


*

*the shape of a sine wave as $\lambda\to0$;

*a sawtooth as $\lambda\to(90^\circ-23.56^\circ)=\text{latitude of the arctic circle}$.


At $45^\circ$ north (which is where I am sitting right now), the maximum discrepancy between this and a sine wave is less than two minutes.  If you look at the hours for Reykjavik, you see the near-sawtooth: At the end of May, the rate of change is 5 & 1/2 minutes per day, not much less than 6 & 1/2 at the March equinox!
North of the arctic circle, you're taking the arcsine of a number greater than $1$.
A: I'm going to give the same answer as Michael, but phrased somewhat differently: I'll work in a coordinate system whose $xy$-plane is the plane of the orbit of the earth, and whose $z$-axis is the perpendicular to this plane. The summer solstice occurs when the earth's center has largest possible $x$-coordinate, and the polar axis of the earth is therefore tilted so that at the summer solstice, it's in the xz-plane, with the north pole nearer to the sun, and the south pole farther from the sun. 
Letting $\theta = \dfrac{2\pi(\text{date} - \text{summerSolstice})}{365}$, we have the position of the center of the earth as 
$$
C(\theta) = (R \cos \theta, 0, R \sin \theta)
$$
where $R$ is the radius of the earth's orbit, which will not be needed further.
$$
\newcommand{\vv}{\mathbf v}
\newcommand{\uu}{\mathbf u}
\newcommand{\ww}{\mathbf u}
\newcommand{\ey}{\mathbf {e_2}}
$$
The vector along the earth's axis, from S to N, is 
$$
\vv = (-\sin \alpha, 0, \cos \alpha),
$$
where $\alpha = 23.56 \text{ degrees}$ is the inclination of the axis. (You can make $\alpha$ a function of time to include nutation, etc., but I'm not going to do that).  The vector perpendicular to this, in the $xz$ plane, pointing mostly away from the sun at the summer solstice, is
$$
\uu = (\cos \alpha, 0, \sin \alpha).
$$
A third vector, perpendicular to both of these, is the $y$-axis, $\ey$. A typical point $P$ on the earth can be described by its displacement from the earth's center, i.e., by
$$
\ww = P - C(\theta)
$$
If the latitude of $P$, measured between $\pi/2$ at the south pole to $\pi/2$ at the north pole, is $\lambda$, then the vector $\ww$ looks like 
$$
\ww = r\sin (\lambda) [ \sin(t) \ey + \cos(t) \uu] + r\cos (\lambda) \vv
$$
where $t$ is the time of day, measured from $0$ to $2 \pi$, and $r \ll R$ is the radius of the earth. 
Sunrise/sunset occur when the vector $\ww$ is orthogonal to the ray from the origin to the point $P$, i.e., when $\ww$ is orthgonal to $\ww  + C(\theta)$, i.e.
$$
\ww \cdot \ww  + \ww \cdot C(\theta) = 0.
$$
At this point, there are two quite reasonable approximations to be made. 


*

*Pretend that the time $t$ and the date $\theta$ are independent. We fix $\theta$ at the angle for mid-day, on the grounds that it's changing very slowly. 

*The dot product of $\ww$ with itself is small compared to the dot product of $C(\theta)$ with $\ww$, since the magnitude of $C(\theta)$ is the radius of the earth's orbit, while the magnitude of $\ww$ is the radius of the earth, so pretend that the $\uu \cdot \uu$ term is missing in the equation. 
The first of these (and perhaps the second) is used in Michael Hardy's solution. As I recall from my study of celestial navigation, the impact is fairly small, although it does vary from day to day once we no longer assume a circular orbit, so it's not to be ignored if you want ultimate accuracy. The second introduces an error in angle no larger than $r/R$, which is something like $1/2000$ radian, although this is a coarse overestimate; that's about 1.7 minutes of arc, which would turn into about 6 minutes of time. 
At any rate, one can write $t$ in terms of $\theta$: multiply $\theta$ by 365; divide by $2 \pi$, and then take the result mod 1, and multiply by $2 \pi$ to get $t$. 
Once you've done that, the equation above involves only $\theta$, and you can solve for when it's zero (although you probably have to resort to numerical means; exact trig isn't likely to work out nicely). 
The only real difference of this from the previous solution is that time and date get related, and the computation is done with dot products rather than classical geometry. 
