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!