Download or view cambridgetempFourier.frink in plain text format
// This file gives the average temperature in Cambridge, MA for any moment of
// the year. It's based on some Fourier analysis I did of temperatures over
// a 30-year span. Thanks to Sarang Gupta and David Bergman for providing
// raw temperature data. It includes all Fourier terms with magnitude greater
// than 0.2 degree F. This data is very likely overfitted to the noise in the
// original data (which you'll see if you plot it using
// plotCambridgeTemp.frink.) You would think that it would probably be better
// and faster to remove more terms, but when you remove too many high and
// low terms, the Fourier series creates unreliable high and low temperatures.
cambridgeTemp[date] :=
{
startOfYear = beginningOfYear[date, "US/Eastern"]
// t is the time elapsed since the beginning of the year, as a unit with
// dimensions of time. We can thus get continuous variations over the
// course of a day.
t = date - startOfYear
/* The coefficients below can be converted into the period of each
oscillation. For example, the term with the coefficient 1/(8784 hours)
corresponds to a frequency of 1/(8784 hours) or once/year. (There are
8784 hours in 366 days.) The period can be obtained by inverting.
The coefficient 1/(8784 hours) corresponds to a period of (8784 hours)/1
or 366 days.
The next-largest coefficient with coefficient 366/(8784 hours)
corresponds to a frequency of 366/(8784 hours) or a period of
(8784 hours)/366 or 1 day, (which is what you would hope you would see!)
This gives a daily swing of +/- 9.25 degrees F around the daily average.
This means that the average daily temperature varies about 18.5 degrees
F from coldest to hottest.
So, we see that:
* the average temperature over a whole year is 57.94 F,
* the annual temperature average varies around this by plus and
minus 21.0998 degrees F (see the component with frequency of
1/(8784 hours)),
* the daily temerature varies around the daily average by plus and minus
9.25 degrees F (see the component with frequency of 366/(8784 hours),
which corresponds to a period of (8784 hours)/366 which is equal to
24 hours.)
* Interestingly, the third-largest component has a period of
(8784 hours) / 732, which is 12 hours or 1/2 of a day.
(Previous versions of this file got the calculation wrong and stated
that the period was 12 days. 12 hours is much less surprising.)
* We can calculate the maximum daily temperature rate increase by
taking the derivative of the daily term (see the component with
frequency of 366/(8784 hours), which corresponds to a period T of
(8784 hours)/366 which is equal to 24 hours.
It has an amplitude "a" of 9.25 degrees F)
To make it simpler, we treat the term as a sin function, ignoring the
phase:
a sin[2 pi f t]
or
a sin[2 pi t/T]
Taking the derivative of this, we get
2 pi a T^-1 cos[2 pi t / T]
and taking the derivative at t = 0 s, the cos term becomes 1, giving:
2 pi a / T
or
2 pi a f
Which, for the daily term (T=24 hours), gives
2.422 degF/hour
Thus, at the time of day that the temperature is increasing the fastest,
we should expect it to be increasing by 2.422 degF per hour. This is
an interesting result. The Earth would probably have smaller daily
temperature swings if it rotated faster, as it did in the past. This
result gives you a way of potentially quantifying that.
* What is encoded by the phase? For example, in the annual term,
21.099806 cos[2 pi 1/(8784 hours) t - -2.906732]
The "-2.906732" encodes the phase. The cosine function normally
has its maximum at 0, but subtracting a number moves the peak to
the *right*. Adding a number moves the peak to the left. Since the
case here has two negative signs, it's like adding, and moves the
peak to the left. The phase will be between -pi and pi.
(Again, the term with the coefficient 1/(8784 hours)
corresponds to a frequency of 1/(8784 hours) or once/year. (There are
8784 hours in 366 days.) The period can be obtained by inverting.
The coefficient 1/(8784 hours) corresponds to a period T of
(8784 hours)/1
or 366 days.)
The amount the cosine curve is shifted *right* is given by:
-T phase / (2 pi)
Which, for the above case, gives
-year * -2.906732 / (2 pi)
which is -168.92 days.
To make this always positive, you want to take the result mod T:
((-T phase) / (2 pi)) mod T
which gives 196.31 days after the beginning of the year (corresponding
to July 16 08:26 AM)
*/
return F[57.942890 +
21.099806 cos[2 pi 1/(8784 hours) t - -2.906732] +
1.356891 cos[2 pi 2/(8784 hours) t - 2.233304] +
1.336084 cos[2 pi 3/(8784 hours) t - 2.193206] +
0.223166 cos[2 pi 4/(8784 hours) t - -1.090373] +
0.521739 cos[2 pi 5/(8784 hours) t - 1.560452] +
0.492775 cos[2 pi 6/(8784 hours) t - 1.439721] +
0.263230 cos[2 pi 8/(8784 hours) t - 0.417370] +
0.556417 cos[2 pi 10/(8784 hours) t - 2.246992] +
0.527294 cos[2 pi 11/(8784 hours) t - 1.969223] +
0.254602 cos[2 pi 13/(8784 hours) t - -2.251578] +
0.217195 cos[2 pi 16/(8784 hours) t - -0.332239] +
0.316010 cos[2 pi 17/(8784 hours) t - 0.773907] +
0.444743 cos[2 pi 18/(8784 hours) t - 2.923433]+
0.265131 cos[2 pi 19/(8784 hours) t - 3.096506] +
0.386864 cos[2 pi 22/(8784 hours) t - 2.064502] +
0.213935 cos[2 pi 23/(8784 hours) t - 3.138032] +
0.209864 cos[2 pi 25/(8784 hours) t - -2.698274] +
0.210136 cos[2 pi 26/(8784 hours) t - -0.895948] +
0.386686 cos[2 pi 28/(8784 hours) t - -0.020524] +
0.263972 cos[2 pi 29/(8784 hours) t - -1.972257] +
0.273588 cos[2 pi 30/(8784 hours) t - -1.537358] +
0.382904 cos[2 pi 31/(8784 hours) t - -0.290935] +
0.459307 cos[2 pi 32/(8784 hours) t - 1.758215] +
0.248827 cos[2 pi 33/(8784 hours) t - 2.639307] +
0.321330 cos[2 pi 34/(8784 hours) t - 0.734143] +
0.261203 cos[2 pi 36/(8784 hours) t - -0.526883] +
0.382089 cos[2 pi 40/(8784 hours) t - 1.726138] +
0.277952 cos[2 pi 42/(8784 hours) t - -1.709813] +
0.351550 cos[2 pi 47/(8784 hours) t - -2.988918] +
0.228715 cos[2 pi 50/(8784 hours) t - 0.349665] +
0.290932 cos[2 pi 55/(8784 hours) t - -1.062905] +
0.250276 cos[2 pi 56/(8784 hours) t - 0.899870] +
0.215393 cos[2 pi 57/(8784 hours) t - 1.793273] +
0.345179 cos[2 pi 62/(8784 hours) t - -0.519825] +
0.208761 cos[2 pi 63/(8784 hours) t - -1.109292] +
0.246095 cos[2 pi 73/(8784 hours) t - -0.547546] +
0.297484 cos[2 pi 76/(8784 hours) t - -0.249699] +
0.223039 cos[2 pi 79/(8784 hours) t - -2.698298] +
0.297275 cos[2 pi 87/(8784 hours) t - -1.849753] +
0.281730 cos[2 pi 100/(8784 hours) t - 1.104291] +
0.299084 cos[2 pi 113/(8784 hours) t - -2.426574] +
0.242557 cos[2 pi 363/(8784 hours) t - -2.234784] +
0.283761 cos[2 pi 364/(8784 hours) t - 1.835604] +
0.429462 cos[2 pi 365/(8784 hours) t - 2.729831] +
9.249708 cos[2 pi 366/(8784 hours) t - -0.555028] +
1.055960 cos[2 pi 367/(8784 hours) t - 1.931133] +
0.392247 cos[2 pi 368/(8784 hours) t - -2.936348] +
0.243464 cos[2 pi 369/(8784 hours) t - 1.750387] +
0.426907 cos[2 pi 731/(8784 hours) t - -1.320502] +
1.805608 cos[2 pi 732/(8784 hours) t - -2.232321] +
0.299388 cos[2 pi 733/(8784 hours) t - -2.079456] +
0.369901 cos[2 pi 1098/(8784 hours) t - -0.200385] +
0.228087 cos[2 pi 1099/(8784 hours) t - 2.390151] +
0.280292 cos[2 pi 1464/(8784 hours) t - -1.966137]]
}
Download or view cambridgetempFourier.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 20217 days, 15 hours, 17 minutes ago.