Calculating length of summer solstice day using tilt and latitude I want to be able to predict how long the summer solstice lasts at any part of world given the maximum tilt angle $φ = 0.4101524 \; \mathrm{rad}$ and the latitude $θ \; \mathrm{rad}$. I specifically want to do this using a geometric approach. The problem is, I never had experience with spherical trig before, and the values I got are somewhat off from their real values. So I want to know if that's because I didn't do the math correctly, and/or because of physical reasons (e.g. the Earth isn't perfectly spherical, the tilt changes slightly etc.).
Here is how I did it:

Let the Earth be represented as a perfect sphere of radius $R$ moving at a constant angular velocity. Let the Earth's rotational axis be in the z direction, and the axis of the circle forming the dark/light boundary be tilted by angle $φ$ relative to the Earth's axis. During the summer solstice, our location of interest (say, in the northern hemisphere) will trace out a circular path of radius $r &lt; R$ perpendicular to the z-axis. These three circles, $C_{\text{rot}}$, $C_{\text{bound}}$ and $C_{\text{path}}$ create a spherical triangle. Using the arc length $s$ of this triangle along $C_{\text{path}}$, we can calculate the length of the day in hours using: $$24 \times \frac{πr + 2s}{2πr}$$

The point where $C_{\text{rot}}$ and $C_{\text{path}}$ meet is located at $$P(\text{rot/path}) = (R\cosθ, 0, R\sinθ)$$

To find the point where $C_{\text{bound}}$ and $C_{\text{path}}$ meet, we first observe that $$y = \sqrt{R^2 - x^2}\sinφ$$ and $$z = \sqrt{R^2 - x^2}\cosφ$$
Which we can solve for $x$ by setting $z$ equal to the $z$ from the previous point
$$x = R\sqrt{1 - \frac{\sin^2θ}{\cos^2φ}}$$
Plugging this back into $y$ and $z$, we finally get
$$P(\text{bound/path}) = \left (R\sqrt{1 - \frac{\sin^2θ}{\cos^2φ}},R\sinθ\tanφ,R\sinθ \right )$$
Now these two points are on the same plane in the z-axis, which is the plane of $P_{\text{path}}$. The distance $d$ between these points is 
$$d = \sqrt{\left (R\sqrt{1 - \frac{\sin^2θ}{\cos^2φ}} - R\cosθ \right)^2 + R^2\sinθ^2\tan^2φ}$$
From this we can get $s$ by using the arc length formula:
$$s = 2\arcsin \left (\frac{d}{2r} \right )r$$
I used Mathematica to compare the calculated and expected values in 2017:

New York: $θ = 0.7105724077$; length (calculated) = 14 hours 56 minutes; length (expected) = 15 hours, 5 minutes

London: $θ = 0.8989737191$; length (calculated) = 14 hours 25 minutes; length (expected) = 16 hours, 38 minutes

As you can see, there is a noticeable difference (10+ minutes).
If you are reducing everything to plain geometry, you can compute the fraction of the "sunny" arc of the circle of latitude for a given location. (Hint: draw a diagram showing the tropic of Cancer, the circle of latitude, the terminator and the axis, all as straight lines.) The height of this segment of a circle of latitude (at $\lambda$) is

$$h = \frac{\cos(\lambda) + \sin(\lambda)\times \tan(\epsilon)}{\cos(\lambda)}$$

The angle of this segment is

$$a = 2\times \arccos(1 - h)$$

I get 


New York: 14.9297092807291
London: 16.4195095586004


which is fairly close to the "expected" numbers you have provided. The difference is very likely due to Mathematica using a geoid rather than a sphere and including atmospheric refraction.