Star Altitude / Azimuth calculation script

For the discussion of the sciences. Physics problems, chemistry equations, biology weirdness, it all goes here.

Moderators: gmalivuk, Moderators General, Prelates

User avatar
Minstrel
Posts: 55
Joined: Thu Jun 26, 2008 4:07 pm UTC

Star Altitude / Azimuth calculation script

Postby Minstrel » Fri Aug 11, 2017 6:37 pm UTC

Hi all, hoping for some feedback on a tiny library I wrote up to calculate the altitude and azimuth of a fixed star, given terrestrial coordinates.

The script is on github here, relevant parts pasted below.

Code: Select all

#!/usr/bin/ruby -w
# encoding: utf-8

require 'date'

class Numeric
  def to_rad
    self * Math::PI / 180
  end
  def to_deg
    self * ( 180 / Math::PI )
  end
end

class MyStars

  attr_accessor :jd, :jda, :t, :gmst, :gast, :last, :lat
 
  # Initialize method takes the local latitude and longitude (as decimal
  # degrees) as input and using the current time creates an object
  # containing the apparent sidereal time, both local and Greenwich.
  # This calculation is based on the USNO's calculations here:
  # http://aa.usno.navy.mil/faq/docs/GAST.php

  def initialize(local_lon, local_lat)
    # Local latitude setter, north is positive
    @lat = local_lat
    # Express local_lon as negative for west, positive for east
    @lon = local_lon
    # Current Julian Day, fractional
    @jd = DateTime.now.ajd.to_f
    # Julian Days since 1 Jan 2000 at 12 UTC
    @jda = @jd - 2451545.0
    # Julian centuries since 1 Jan 2000 at 12 UTC
    # Currently not using this in favor of a quick calculation, w/ loss of 0.1 sec / century.
    @t = @jda / 36525.0
    # Greenwhich Mean Sidereal Time
    # Quick and dirty version, with loss as mentioned above
    @gmst = 18.697374558 + ( 24.06570982441908 * @jda )
    # Reduce to a range of 0 - 24 h
    @gmst = @gmst % 24
    # Greenwich Apparent Sidereal Time
    omega = 125.04 - 0.052954 * @jda
    l = 280.47 + 0.98565 * @jda
    epsilon = 23.4393 - 0.0000004 * @jda
    deltapsi = ( -0.000319 * Math::sin(omega.to_rad) ) - ( 0.000024 * Math::sin( (2*l).to_rad ) )
    eqeq = deltapsi * Math.cos(epsilon.to_rad)
    @gast = @gmst + eqeq
    # Local Apparent Sidereal Time
    @last = @gast + ( local_lon / 15.0 )
  end

  # Altitude and Azimuth methods take the Right Ascension and Declination of a fixed
  # star as decimal hours for RA and decimal degrees for Dec and output the
  # Altitude and Azimuth as decimal degrees.
  # AA method just runs them both and sends a pretty output.
  # This is based on the USNO's calculations here:
  # http://aa.usno.navy.mil/faq/docs/Alt_Az.php
  # Results so far test out within a few minutes of Stellarium.

  def altitude(ra, dec)
    lha = ( @gast - ra ) * 15 + @lon
    a = Math::cos(lha.to_rad)
    b = Math::cos(dec.to_rad)
    c = Math::cos(@lat.to_rad)
    d = Math::sin(dec.to_rad)
    e = Math::sin(@lat.to_rad)
    Math::asin(a*b*c+d*e).to_deg
  end

  def azimuth(ra, dec)
    lha = ( @gast - ra ) * 15 + @lon
    a = -(Math::sin(lha.to_rad))
    b = Math::tan(dec.to_rad)
    c = Math::cos(@lat.to_rad)
    d = Math::sin(@lat.to_rad)
    e = Math::cos(lha.to_rad)
    Math::atan(a/(b*c-d*e)).to_deg
  end

  def aa(ra,dec)
    puts "Altitude is " + self.altitude(ra,dec).to_s
    puts "Azimuth is " + self.azimuth(ra,dec).to_s
  end

end


Example of use if you really want to run it (although I was mostly looking for feedback on the equations themselves used in the 'azimuth' and 'altitude' methods):

Code: Select all

require_relative 'mystars'
# Create a new MyStars object, passing in local longitude and latitude as
# decimal degrees
mystars = MyStars.new(-71, 43)
# Pass Right Ascension (as decimal hours) and Declination (as decimal degrees)
# to the methods altitude, azimuth and get back the corresponding value,
# as decimal degrees.
# The aa method passes them both back as pretty text.
# Calculating Vega's position below:
mystars.azimuth(18.616666,38.783333)
=> 43.96040907459006
mystars.altitude(18.616666,38.783333)
=> 8.97361521398292
mystars.aa(18.616666,38.783333)
Altitude is 8.97361521398292
Azimuth is 43.96040907459006


So far, it produces output within a few arcminutes of 'live' data from Stellarium, which is good enough for what I want. But I'm wondering if anyone who knows trig or astronomical calculations better than I do (read as: anyone who knows trig) can tell me if I'm missing out on anything, like periodicity that should be taken into account with the arcsin and arctan.

If anyone's curious, I'm hoping to use this as the basis for an ncurses based planetarium, because a) I want to learn ncurses (or the ruby accessor libraries thereof) and b) CLI all the things.

Edit: this is working pretty well now. I changed the azimuth calculation to use the double argument arctan, adding 360 to negative values, which now properly returns values of 0-360. Maybe I'll actually have something interesting to show off in a while.

Thanks!

User avatar
Jplus
Posts: 1692
Joined: Wed Apr 21, 2010 12:29 pm UTC
Location: Netherlands

Re: Star Altitude / Azimuth calculation script

Postby Jplus » Thu Aug 17, 2017 9:05 pm UTC

Hi Minstrel, thanks for sharing your little project. I think it looks fun and I starred it on GitHub.

To be honest, it was not clear to me what was your question until I read the edit. A little tip for the next time when you seek help with your code: contrast the expected output with the actual output. Whenever possible, also try to pinpoint a narrower chunk of your code where you suspect the mismatch originates.

Also, exploit the social web! GitHub has the issue list which is meant specifically to describe challenges and wishes related to your code. People who watch your repo (like me now) will be notified when you add a new issue ticket or when somebody comments to an existing one. They might even offer you a pull request with a fix. You can still post a link to your issue here on the forums, in order to attract attention. The more you use GitHub, the more your network of online coding acquaintances will grow. As you may already know there is also Stack Overflow, which is particularly great for the most challenging puzzles. You can link back and forth between GitHub and Stack Overflow to amplify the web goodness further.

Cheers!
"There are only two hard problems in computer science: cache coherence, naming things, and off-by-one errors." (Phil Karlton and Leon Bambrick)

coding and xkcd combined

(Julian/Julian's)

User avatar
Minstrel
Posts: 55
Joined: Thu Jun 26, 2008 4:07 pm UTC

Re: Star Altitude / Azimuth calculation script

Postby Minstrel » Mon Aug 21, 2017 12:35 pm UTC

Jplus, thanks for the advice and the star. It's my first time asking for feedback on code and I was too chicken to post on SO. I'll try to keep any future code more concise - actually I really could have just posted the equations rather than the code come to think of it.

If anyone's interested, I've gotten up to where it's displaying an ncurses representation of a given 10 degree FOV around any star visible in the user's sky. Just kind of a single snapshot right now but working on being able to pan around, magnify, filter by magnitude, ie. the usual planetarium stuff next.

Edit: you can now pan around with arrow keys or the keypad, and zoom in / out with + and -. Enter exits.

It's stretched vertically due to using a rectangular console font, but clear enough to get the point below (pink lines were drawn in after):

Lyra:
Image

Corona borealis:
Image

User avatar
Jplus
Posts: 1692
Joined: Wed Apr 21, 2010 12:29 pm UTC
Location: Netherlands

Re: Star Altitude / Azimuth calculation script

Postby Jplus » Wed Aug 30, 2017 7:04 pm UTC

I'm interested! Just replying late because I have so little time to lurk the forums.

Your screenshots look awesome.
"There are only two hard problems in computer science: cache coherence, naming things, and off-by-one errors." (Phil Karlton and Leon Bambrick)

coding and xkcd combined

(Julian/Julian's)

User avatar
Minstrel
Posts: 55
Joined: Thu Jun 26, 2008 4:07 pm UTC

Re: Star Altitude / Azimuth calculation script

Postby Minstrel » Thu Aug 31, 2017 3:27 pm UTC

Thanks for the kind words!

I've been chipping away at this daily and there is now a side panel with current settings, basic info on the currently selected star, sky is updated every 5 seconds and some other stuff..

What I'd really like to do is get my hands on a scope with a motor control so I can turn this into a bare bones, low resource, telescope controller with maybe some equally basic evening planning. But as my budget doesn't permit that and free telescopes with motor mounts aren't plentiful, I'm moving on to the display end.

Drawing the sky based on altitude and azimuth was pretty easy (basically just set the origin to be the top of the sky, and then calculate x/y points on a circle with radius = (90 - altitude) around it, then some fiddling to show the current field of view).

Shifting this to be a true point of view like big software like Stellarium turns out to be hard. Either hard, or I'm overthinking it like crazy. I feel like I should just be able to say "OK, I'm facing an azimuth angle of 135 degrees and an altitude of 15 degrees, so add or subtract those from the true values of the stars in the database when doing calculations", but spherical coordinates don't work that way (altitude of 80 and any given azimuth are pretty close together, but if you simply increase the effective altitude like I was planning, suddenly two nearby stars are now way across the sky from each other).

What it seems like I need to do is find the actual angular distance between the new center of view and each object in the database, then the angle to that object (essentially a new altitude and azimuth) and re-plot. I can either do that by doing those calculations directly, or maybe just inputting a fake longitude and latitude somewhere behind the scenes and re-use my existing calculations - I wonder if the two might actually turn out to be the same maths.

I think I'm going to do the latter and see how it goes. If all goes well, my biggest worry is then speed - replotting 4500 stars is working fine now, it happens each time the 5 second timer ticks with no delay and little CPU, but I'm not so sure it will continue to do so when I import more stars and deep sky objects and 4500 becomes 45,000+.

Edit - I accomplished it, in the fake_geo branch, by doing the fake latitude / longitude trick. You can now pan around both hemispheres, but the panning feels very awkward, since it spins you around the north celestial pole. I'm going to try to put in a fake ground layer and compass points to get a better feel for reference points, but not sure that will cut it.

What I really want is a 3d renderer for the terminal (there are some, but they all seem to rely on X or similar, which defeats the purpose for me). Maybe I need to write that first. Down the rabbit hole!


Return to “Science”

Who is online

Users browsing this forum: No registered users and 8 guests