Shallow Thoughts : : Dec

Plotting the Analemma

My SJAA planet-observing column for January is about the Analemma and the Equation of Time.

The analemma is that funny figure-eight you see on world globes in the middle of the Pacific Ocean. Its shape is the shape traced out by the sun in the sky, if you mark its position at precisely the same time of day over the course of an entire year.

The analemma has two components: the vertical component represents the sun's declination, how far north or south it is in our sky. The horizontal component represents the equation of time.

The equation of time describes how the sun moves relatively faster or slower at different times of year. It, too, has two components: it's the sum of two sine waves, one representing how the earth speeds up and slows down as it moves in its elliptical orbit, the other a function the tilt (or "obliquity") of the earth's axis compared to its orbital plane, the ecliptic.

The Wikipedia page for Equation of time includes a link to a lovely piece of R code by Thomas Steiner showing how the two components relate. It's labeled in German, but since the source is included, I was able to add English labels and use it for my article.

But if you look at photos of real analemmas in the sky, they're always tilted. Shouldn't they be vertical? Why are they tilted, and how does the tilt vary with location? To find out, I wanted a program to calculate the analemma.

Calculating analemmas in PyEphem

The very useful astronomy Python package PyEphem makes it easy to calculate the position of any astronomical object for a specific location. Install it with: easy_install pyephem for Python 2, or easy_install ephem for Python 3.

import ephem
observer = ephem.city('San Francisco')
sun = ephem.Sun()
sun.compute(observer)
print sun.alt, sun.az

The alt and az are the altitude and azimuth of the sun right now. They're printed as strings: 25:23:16.6 203:49:35.6 but they're actually type 'ephem.Angle', so float(sun.alt) will give you a number in radians that you can use for calculations.

Of course, you can specify any location, not just major cities. PyEphem doesn't know San Jose, so here's the approximate location of Houge Park where the San Jose Astronomical Association meets:

observer = ephem.Observer()
observer.name = "San Jose"
observer.lon = '-121:56.8'
observer.lat = '37:15.55'

You can also specify elevation, barometric pressure and other parameters.

So here's a simple analemma, calculating the sun's position at noon on the 15th of each month of 2011:

    for m in range(1, 13) :
        observer.date('2011/%d/15 12:00' % (m))
        sun.compute(observer)

I used a simple PyGTK window to plot sun.az and sun.alt, so once it was initialized, I drew the points like this:

    # Y scale is 45 degrees (PI/2), horizon to halfway to zenith:
    y = int(self.height - float(self.sun.alt) * self.height / math.pi)
    # So make X scale 45 degrees too, centered around due south.
    # Want az = PI to come out at x = width/2.
    x = int(float(self.sun.az) * self.width / math.pi / 2)
    # print self.sun.az, float(self.sun.az), float(self.sun.alt), x, y
    self.drawing_area.window.draw_arc(self.xgc, True, x, y, 4, 4, 0, 23040)

So now you just need to calculate the sun's position at the same time of day but different dates spread throughout the year.

And my 12-noon analemma came out almost vertical! Maybe the tilt I saw in analemma photos was just a function of taking the photo early in the morning or late in the afternoon? To find out, I calculated the analemma for 7:30am and 4:30pm, and sure enough, those were tilted.

But wait -- notice my noon analemma was almost vertical -- but it wasn't exactly vertical. Why was it skewed at all?

Time is always a problem

As always with astronomy programs, time zones turned out to be the hardest part of the project. I tried to add other locations to my program and immediately ran into a problem.

The ephem.Date class always uses UTC, and has no concept of converting to the observer's timezone. You can convert to the timezone of the person running the program with localtime, but that's not useful when you're trying to plot an analemma at local noon.

At first, I was only calculating analemmas for my own location. So I set time to '20:00', that being the UTC for my local noon. And I got the image at right. It's an analemma, all right, and it's almost vertical. Almost ... but not quite. What was up?

Well, I was calculating for 12 noon clock time -- but clock time isn't the same as mean solar time unless you're right in the middle of your time zone.

You can calculate what your real localtime is (regardless of what politicians say your time zone should be) by using your longitude rather than your official time zone:

    date = '2011/%d/12 12:00' % (m)
    adjtime = ephem.date(ephem.date(date) \
                    - float(self.observer.lon) * 12 / math.pi * ephem.hour)
    observer.date = adjtime

Maybe that needs a little explaining. I take the initial time string, like '2011/12/15 12:00', and convert it to an ephem.date. The number of hours I want to adjust is my longitude (in radians) times 12 divided by pi -- that's because if you go pi (180) degrees to the other side of the earth, you'll be 12 hours off. Finally, I have to multiply that by ephem.hour because ... um, because that's the way to add hours in PyEphem and they don't really document the internals of ephem.Date.

Set the observer date to this adjusted time before calculating your analemma, and you get the much more vertical figure you see here. This also explains why the morning and evening analemmas weren't symmetrical in the previous run.

This code is location independent, so now I can run my analemma program on a city name, or specify longitude and latitude.

PyEphem turned out to be a great tool for exploring analemmas. But to really understand analemma shapes, I had more exploring to do. I'll write about that, and post my complete analemma program, in the next article.

Tags: analemma, astronomy, science, programming, python, writing
[ 20:54 Dec 29, 2011    More science/astro | permalink to this entry | ]

Calculating the Solstice and shortest day

Today is the winter solstice -- the official beginning of winter.

The solstice is determined by the Earth's tilt on its axis, not anything to do with the shape of its orbit: the solstice is the point when the poles come closest to pointing toward or away from the sun. To us, standing on Earth, that means the winter solstice is the day when the sun's highest point in the sky is lowest.

You can calculate the exact time of the equinox using the handy Python package PyEphem. Install it with: easy_install pyephem for Python 2, or easy_install ephem for Python 3. Then ask it for the date of the next or previous equinox. You have to give it a starting date, so I'll pick a date in late summer that's nowhere near the solstice:

>>> ephem.next_solstice('2011/8/1')
2011/12/22 05:29:52

That agrees with my RASC Observer's Handbook: Dec 22, 5:30 UTC. (Whew!)

PyEphem gives all times in UTC, so, since I'm in California, I subtract 8 hours to find out that the solstice was actually last night at 9:30. If I'm lazy, I can get PyEphem to do the subtraction for me:

ephem.date(ephem.next_solstice('2011/8/1') - 8./24)
2011/12/21 21:29:52

I used 8./24 because PyEphem's dates are in decimal days, so in order to subtract 8 hours I have to convert that into a fraction of a 24-hour day. The decimal point after the 8 is to get Python to do the division in floating point, otherwise it'll do an integer division and subtract int(8/24) = 0.

The shortest day

The winter solstice also pretty much marks the shortest day of the year. But was the shortest day yesterday, or today? To check that, set up an "observer" at a specific place on Earth, since sunrise and sunset times vary depending on where you are. PyEphem doesn't know about San Jose, so I'll use San Francisco:

>>> import ephem
>>> observer = ephem.city("San Francisco")
>>> sun = ephem.Sun()
>>> for i in range(20,25) :
...   d = '2011/12/%i 20:00' % i
...   print d, (observer.next_setting(sun, d) - observer.previous_rising(sun, d)) * 24
2011/12/20 20:00 9.56007901422
2011/12/21 20:00 9.55920379754
2011/12/22 20:00 9.55932991847
2011/12/23 20:00 9.56045709446
2011/12/24 20:00 9.56258416496

I'm multiplying by 24 to get hours rather than decimal days.

So the shortest day, at least here in the bay area, was actually yesterday, 2011/12/21. Not too surprising, since the solstice wasn't that long after sunset yesterday.

If you look at the actual sunrise and sunset times, you'll find that the latest sunrise and earliest sunset don't correspond to the solstice or the shortest day. But that's all tied up with the equation of time and the analemma ... and I'll cover that in a separate article.

Tags: astronomy, science, programming, python, writing
[ 11:28 Dec 22, 2011    More science/astro | permalink to this entry | ]

Convert patterns in only some lines to title case

A friend had a fun problem: she had some XML files she needed to import into GNUcash, but the program that produced them left names in all-caps and she wanted them more readable. So she'd have a file like this:

<STMTTRN>
   <TRNTYPE>DEBIT
   <DTPOSTED>20111125000000[-5:EST]
   <TRNAMT>-22.71
   <FITID>****

   <NAME>SOME    COMPANY
   <MEMO>SOME COMPANY    ANY TOWN   CA 11-25-11 330346
</STMTTRN>

and wanted to change the NAME and MEMO lines to read Some Company and Any Town. However, the tags, like <NAME>, all had to remain upper case, and presumably so did strings like DEBIT. How do you change just the NAME and MEMO lines from upper case to title case?

The obvious candidate to do string substitutes is sed. But there are several components to the problem.

Addresses

First, how do you ensure the replacement only happens on lines with NAME and MEMO?

sed lets you specify address ranges for just that purpose. If you say sed 's/xxx/yyy/' sed will change all xxx's to yyy; but if you say sed '/NAME/s/xxx/yyy/' then sed will only do that substitution on lines containing NAME.

But we need this to happen on lines that contain either NAME or MEMO. How do you do that? With \|, like this: sed '/\(NAME\|MEMO\)/s/xxx/yyy/'

Converting to title case

Next, how do you convert upper case to lower case? There's a sed command for that: \L. Run sed 's/.*/\L&/' and type some upper and lower case characters, and they'll all be converted to lower-case.

But here we want title case -- we want most of each word converted to lowercase, but the first letter should stay uppercase. That means we need to detect a word and figure out which is the first letter.

In the strings we're considering, a word is a set of letters A through Z with one of the following characteristics:

1. It's preceded by a space
2. It's preceded by a close-angle-bracket, >

So the pattern /[ >][A-Z]*/ will match anything we consider a word that might need conversion.

But we need to separate the first letter and the rest of the word, so we can treat them separately. sed's \( \) operators will let us do that. The pattern \([ >][A-Z]\) finds the first letter of a word (including the space or > preceding it), and saves that as its first matched pattern, \1. Then \([A-Z]*\) right after it will save the rest of the word as \2.

So, taking our \L case converter, we can convert to title case like this: sed 's/\([ >][A-Z]\)\([A-Z]*\)/\1\L\2/g

Starting to look long and scary, right? But it's not so bad if you build it up gradually from components. I added a g on the end to tell sed this is a global replace: do the operation on every word it finds in the line, otherwise it will only make the substitution once, on the first word it sees, then quit.

Putting it together

So we know how to seek out specific lines, and how to convert to title case. Put the two together, and you get the final command:

sed '/\(NAME\|MEMO\)/s/\([ >][A-Z]\)\([A-Z]*\)/\1\L\2/g'

I ran it on the test input, and it worked just fine.

For more information on sed, a good place to start is the sed regular expressions manual.

Tags: regexp, cmdline, sed
[ 14:13 Dec 18, 2011    More linux/cmdline | permalink to this entry | ]

Taming the loud system beep

One of the distros I'm trying on my Dell Latitude 2120 laptop is Arch Linux. I like a lot of things about Arch; I had to stop using it on my previous laptop because something broke in the wi-fi drivers, but I always regretted giving it up. And it seems to work quite well on the Dell, with one exception: the system beep (which other distros don't support at all) was ear-shatteringly loud, far too loud to consider being able to use this laptop in a public space. Clearly that needed to be fixed.

The usual approach to system beeps, unloading or blacklisting the pcspkr module, had no effect. xset b off turned the beep off in X; I could also set its pitch and duration to change it to a nice quiet click, though I wasn't able to change the volume that way. I actually do like having a system beep, as long as it's fairly quiet and won't disturb people nearby. Unfortunately, xset b only affects the bell while in X; it didn't have any effect on the deafening sound Arch gave upon shutdown or reboot.

It turns out that on some laptops, including this Dell, the system beep goes not through the old-style pcspkr driver, but through the normal sound card. And the sound card has a separate channel for the system beep, so even if you have your volume turned down, the beep may still be at 100%. All I needed to do was run alsamixer and find out what the channel was called: "Beep".

Given that, I could use the amixer program to ensure the beep volume will be sane when I log in. I added the following to .zlogin (for zsh; obviously, adjust for your own shell):

amixer -q -c 0 set Beep 5

That gave me a nice quiet beep. If I need to turn it off completely, amixer can do that too: amixer -q -c 0 set Beep mute (curiously, amixer -q -c 0 set Beep 0 doesn't actually set the volume to zero, just sets it very low).

That volume setting applies to the shutdown beep, too, fortunately. Though what I'd really like is to have quiet beeps while I'm running, but no shutdown beep at all. I don't understand the purpose of the shutdown beep; obviously I know when I've told my machine to shut down or reboot, so why do I need an audible reminder? But I've been unable to find anything explaining what's causing this beep. I tried adding a amixer -q -c 0 set Beep mute to /etc/rc.local.shutdown, but it didn't help; apparently the shutdown beep is called before that file is run. Which strongly suggests it is being run by Arch Linux, not by something in the BIOS. But nobody I've asked had any suggestions as to its source, or how to change it. Another enduring Linux mystery ...

Update: I mentioned xset b as a way to adjust beeps inside X -- in addition to turning beeps totally off, you can also set pitch, duration, and sometimes volume though the volume part didn't work for me. But outside X, you can make similar adjustments with setterm, e.g. setterm -bfreq 400 -blength 50. Thanks to Mikachu for the tip!

Tags: linux, archlinux, alsa
[ 12:05 Dec 13, 2011    More linux | permalink to this entry | ]

Set Ubuntu's system clock to use localtime, not UTC

Need your Ubuntu clock to stay in sync with a dual-boot Windows install? It seems to have changed over the years, and google wasn't finding any pages for me offering suggestions that still worked.

Turns out, in Ubuntu Oneiric Ocelot, it's controlled by the file /etc/default/rcS: set UTC=no.

Apparently it's also possible to get Windows to understand a UTC system clock using a registry tweak.

Ironically, that page, which I found by searching for windows system clock utc, also has the answer for setting Ubuntu to local time. So if I'd searched for Windows in the first place, I wouldn't have had to puzzle out the Ubuntu solution myself. Go figure!

Tags: linux, ubuntu, tip
[ 13:56 Dec 11, 2011    More linux | permalink to this entry | ]