It’s common knowledge that the length of a day (time from sunrise to sunset) varies over the year and with one’s location on the earth. Given only the following, can you derive an approximate equation for the length of a day (in hours) of any specific location on the earth, on any given day of the year?
1) The inclination of the earth’s rotational axis is 23.45 degrees
2) The length of a day is exactly 24 hours
3) The length of a year is exactly 365 days
4) Location on earth is given by the latitude
You are allowed to ignore secondary (but real!) effects such as the earth’s non-circular orbit, the sun being a disc, refraction of sunlight by the atmosphere, etc. To allow easier comparison of different solutions, let’s also assume that North Latitude is positive, and that the Winter Solstice in the Northern Hemisphere is “day 0” (Hint: therefore also day 365!) of the year.
The first thing you need is the distance of the sun from the celestial equator in degrees; this is called the declination, so we'll use d to represent it. So we'll use D to represent the day past the December solstice.
The declination, d, depends on the 23.45° tilt of the ecliptic (the sun's path on the celestial sphere) and the sun's celestial longitude, l, which is (D*360/365 - 90°) as it is measured from the March equinox.
Here we start with the spherical trig. (Later we'll simplify for a differnt formula to skip over this spherical trig portion, but even then we'll need the spherical trig in a later step.) We know the 23.45° angle where the ecliptic meets the equator, and the 90° angle where the sought side of size d meets the equator, and l, which is the hypotenuse of this spherical right triangle. The spherical law of sines tells us that sin(l)/sin(90°) = sin(d)/sin(23.45°), so
d = arcsin(sin(l)*sin(23.45°))
This arc, d, is the complement of the distance of the sun from the north celestial pole, which is one side of the triangle we'll need to use. That triangle's vertices are the north pole, the sun and the zenith of the location at the time of sunrise or sunset. The relevant information for this triangle is:
Zenith to sun at sunrise/set = 90°
North pole to sun = 90° - d (this is more than 90° when d is negative)
North pole to observer's zenith = 90° - L, where L is the observer's latitude.
Here we get to the unavoidable spherical trig, via the spherical law of cosines, as we need to know the angle SPZ (Sun, Pole, Zenith), called the hour angle, so we'll call its measure h, rather than P:
cos(90°) = cos(90°-L)*cos(90°-d) + sin(90°-L)*sin(90°-d)*cos(h)
where the 90° on the LHS is the desired distance from the zenith.
Rearranging, substituting 0 for cos(90) and switching sin and cos to avoid complementing:
cos(h) = -sin(L)*sin(d) / (cos(L)*cos(d)) = -tan(L)*tan(d)
If h is measured in degrees and you want time measured in hours,
t = h/15
but that's only half a day, from local noon to sunset (or sunrise), so we need twice that: 2*h/15.
Working backward:
Day length =
2*arccos(-tan(L)*tan(d)) / 15
or
2*arccos(-tan(L)*tan(arcsin(sin(l)*sin(23.45°)))) / 15
or
2*arccos(-tan(L)*tan(arcsin(sin(D*360/365 - 90)*sin(23.45°)))) / 15
Simplification:
As mentioned before, the first application of spherical trig could be substituted with a simplified formula for the distance of the sun north of the celestial equator: 23.45*sin((D-90)*(360/365))--. Then you'd have
2*arccos(-tan(L)*tan(23.45*sin((D-90)*(360/365)))) / 15
clearvars,clc
L=40;
for D=0:10:360 % in days
formula1=2*acosd(-tand(L)*tand(asind(sind(D*360/365 - 90)*sind(23.45)))) / 15;
formula2=2*acosd(-tand(L)*tand(23.45*sind((D-90)*(360/365)))) / 15;
fprintf('%4d %6.2f %6.2f\n',D,formula1,formula2)
end
compares the two formulae at latitude 40° North:
day original less exact
formula formula
0 9.15 9.15
10 9.21 9.22
20 9.36 9.37
30 9.60 9.61
40 9.90 9.92
50 10.26 10.28
60 10.66 10.68
70 11.08 11.11
80 11.51 11.55
90 11.95 12.00
100 12.38 12.45
110 12.82 12.89
120 13.24 13.32
130 13.64 13.72
140 14.01 14.08
150 14.33 14.39
160 14.59 14.63
170 14.76 14.78
180 14.84 14.85
190 14.82 14.81
200 14.69 14.68
210 14.47 14.46
220 14.18 14.16
230 13.83 13.81
240 13.44 13.42
250 13.03 13.00
260 12.60 12.56
270 12.16 12.11
280 11.73 11.66
290 11.29 11.22
300 10.86 10.78
310 10.46 10.38
320 10.08 10.00
330 9.74 9.68
340 9.47 9.42
350 9.27 9.25
360 9.17 9.16
The main difference from published times comes from atmospheric refraction and the radius of the sun's disk, which can add about five minutes. The eccentricity of the earth's orbit only changes the dates slightly at which the given dates the times are achieved.
The atmospheric refraction and the sun's apparent radius can be incorporated by using (.25° for the sun's apparent radius and .55° for atmospheric refraction):
cos(90.8°) = cos(90°-L)*cos(90°-d) + sin(90°-L)*sin(90°-d)*cos(h)
to determine the hour angle h.
h = arccos((cos(90.8°) - cos(90°-L)*cos(90°-d)) / (sin(90°-L)*sin(90°-d)) )
For comparison the formula incorporating the correction for refraction is appended on the right:
clearvars,clc
L=40;
for D=0:10:360 % in days
formula1=2*acosd(-tand(L)*tand(asind(sind(D*360/365 - 90)*sind(23.45)))) / 15;
formula2=2*acosd(-tand(L)*tand(23.45*sind((D-90)*(360/365)))) / 15;
l=(D*360/365 - 90);
d = asind(sind(l)*sind(23.45));
h= acosd((cosd(90.8) - cosd(90-L)*cosd(90-d)) / (sind(90-L)*sind(90-d)) );
t=2*h/15;
fprintf('%4d %6.2f %6.2f %6.2f\n',D,formula1,formula2,t)
end
day original less exact original with
formula formula correction for
refraction
0 9.15 9.15 9.32
10 9.21 9.22 9.37
20 9.36 9.37 9.52
30 9.60 9.61 9.75
40 9.90 9.92 10.06
50 10.26 10.28 10.41
60 10.66 10.68 10.80
70 11.08 11.11 11.22
80 11.51 11.55 11.65
90 11.95 12.00 12.08
100 12.38 12.45 12.52
110 12.82 12.89 12.96
120 13.24 13.32 13.38
130 13.64 13.72 13.79
140 14.01 14.08 14.16
150 14.33 14.39 14.49
160 14.59 14.63 14.75
170 14.76 14.78 14.93
180 14.84 14.85 15.01
190 14.82 14.81 14.98
200 14.69 14.68 14.85
210 14.47 14.46 14.63
220 14.18 14.16 14.33
230 13.83 13.81 13.98
240 13.44 13.42 13.59
250 13.03 13.00 13.17
260 12.60 12.56 12.74
270 12.16 12.11 12.30
280 11.73 11.66 11.87
290 11.29 11.22 11.43
300 10.86 10.78 11.01
310 10.46 10.38 10.60
320 10.08 10.00 10.23
330 9.74 9.68 9.90
340 9.47 9.42 9.62
350 9.27 9.25 9.43
360 9.17 9.16 9.33
|
|
Posted by Charlie
on 2023-07-04 17:31:41 |
Solution (spoiler)
