# NJP Polar Alignment Tool # Jonathan Tomshine # December 20, 2009 from math import * from time import * from Tkinter import * from tkMessageBox import * # get Julian date from time.struct_time object def get_JD(now): # Algorithm from "Astronomical Formulae for Calculators" by Jean Meeus, 4th Ed. # parse time.struct_time object into integer year (Y), integer month (M) and fractional day (D) # assume that time.struct_time is provided as UTC Y = float(now.tm_year) M = float(now.tm_mon) D = float(now.tm_mday) D = D + float(now.tm_hour)/24 D = D + float(now.tm_min)/(24*60) D = D + float(now.tm_sec)/(24*60*60) if M <= 2: Y = Y -1 M = M + 12 # assume that we are looking at dates after Oct. 15, 1582 -- loss of some generality A = floor(Y/100) B = 2 - A + floor(A/4) JD = floor(365.25*Y) + floor(30.6001*(M+1)) + D + 1720994.5 + B return JD # compute Local Apparent Sidereal Time from JD and longitude def calc_lst(JD, Longitude): # algorithm taken from US Naval Observatory online, Dec. 20, 2009: # http://aa.usno.navy.mil/faq/docs/GAST.php # Let JD be the Julian date of the time of interest. # Let JD0 be the Julian date of the previous midnight (0h) UT (the value of # JD0 will end in .5 exactly), and let H be the hours of UT elapsed since # that time. Thus we have JD = JD0 + H/24. # find JD0 JD0 = floor(JD - 0.5) + 0.5 # find H, D, D0, T H = (JD - JD0)*24 D = JD - 2451545.0 D0 = JD0 - 2451545.0 T = D/36525 # the Greenwich mean sidereal time in hours is GMST = 6.697374558 + 0.06570982441908*D0 + 1.00273790935*H + 0.000026*T*T GMST = fmod(GMST, 24) # now compute some terms to correct the mean sidereal time to apparent # Omega, L, and epsilon (below) are given in degrees # the Longitude of the ascending node of the Moon, given as Omega = 125.04 - 0.052954*D # L, the Mean Longitude of the Sun, given as L = 280.47 + 0.98565*D # the obliquity is given as epsilon = 23.4393 - 0.0000004*D # delta_psi -- nutation in right ascension or the equation of the equinoxes delta_psi = -0.000319*sin(Omega/180*pi) - 0.000024*sin(2*L/180*pi) eqeq = delta_psi*cos(epsilon/180*pi) # correct GMST to GAST GAST = GMST + eqeq # translate GAST to LAST (from Greenwich to Local) LAST = GAST + (Longitude/360)*24 # be sure that we're in the range 0-24 if (LAST < 0): LAST = LAST + 24 else: if (LAST > 24): LAST = LAST - 24 return LAST # get hour angle of Polaris def polaris_ha(LAST, RA): # calcluate hour angle HA = LAST - RA if (HA < 0): HA = HA + 24.0 return HA # get theta for polaris: reticle angle -- last step! def get_theta(HA): # given an hour-angle, get the number of degrees between the right horizontal (x-axis) # and polaris, as measured counter-clockwise # i.e., for HA=6, theta=0, for HA=12, theta=90, for HA=18, theta=180, for HA=0/360, theta=270 theta = (HA/24)*360 - 90 if (theta < 0): theta = theta + 360 return theta # get radius for polaris -- figure out where to put it in the circle def get_radius(now): # uses a simple linear interpolation between the inner (2015, at r=0.9) and outer (1985, at r=1.1) rings # ignore leap-years -- will make no visible difference year = float(now.tm_year) + float(now.tm_yday)/365; m = (0.9-1.1)/(2015-1985) b = 1.1 - 1985*m radius = year*m + b return radius # verify a real prior to analysis def verify_float(float_in, min_val, max_val): error = False error_msg = "(no error)" try: # try the real conversion, see if it throws an error x = float(float_in) except ValueError: error = True error_msg = "Non-numerical input provided where a number is required." display_msg(error_msg) return error # x is a real -- is it within the proper range? if (x < min_val): error = True error_msg = "Value " + str(x) + " is below the minimum value of " + str(min_val) + "." if (x > max_val): error = True error_msg = "Value " + str(x) + " is above the maximum value of " + str(max_val) + "." # if we found an error, display a dialog box and exit if error: display_msg(error_msg) return error # display a message in a warning dialog box def display_msg(message_string): showwarning("Message", message_string) # destroy existing view, re-gather longitude, recalculate position of Polaris and draw a new view def update_figure(): global w global longitude # check if the supplied longitude is valid if (verify_float(longitude_entry.get(), -180.0, 180.0)): data_okay = False else: # destroy existing view window w.destroy() # draw a new one! draw_figure() return # gather longitude and draw the view-window def draw_figure(): global w global longitude # define some critical parameters win_sz = 350 radius = win_sz*0.70 star_radius = win_sz*(3.0/250.0) font_sz = 9 num_ticks = 72; # define right ascension of Polaris for/near current epoch RA = 2.6519444444 # get longitude from main window longitude = float(longitude_entry.get()) # find where to put Polaris now = gmtime() JD = get_JD(now) LAST = calc_lst(JD, longitude) HA = polaris_ha(LAST, RA) theta = get_theta(HA) pol_rad = get_radius(now) # prepare output text reporter_text = "Current UTC: " + asctime(now) + "\n" reporter_text = reporter_text + "Julian Date: " + str(JD) + "\n" reporter_text = reporter_text + "Local Longitude: " + str(longitude) + "\n" reporter_text = reporter_text + "Local Apparent Sidereal Time: " + str(LAST) + "\n" reporter_text = reporter_text + "Local Polaris Hour-Angle: " + str(HA) # draw finder view w = Canvas(root, width=2*win_sz, height=2*win_sz) w.grid(row=0, column=0, columnspan=3) # draw black background w.create_rectangle(0, 0, 2*win_sz, 2*win_sz, fill="black") # draw 3 concentric circles w.create_oval((win_sz-radius*1.1), (win_sz-radius*1.1), (win_sz+radius*1.1), (win_sz+radius*1.1), outline="red", fill="") w.create_oval((win_sz-radius*1.0), (win_sz-radius*1.0), (win_sz+radius*1.0), (win_sz+radius*1.0), outline="red", fill="") w.create_oval((win_sz-radius*0.9), (win_sz-radius*0.9), (win_sz+radius*0.9), (win_sz+radius*0.9), outline="red", fill="") # draw cross-hairs w.create_line(win_sz, (win_sz-radius*1.2), win_sz, (win_sz+radius*1.2), fill="red") w.create_line((win_sz-radius*1.2), win_sz, (win_sz+radius*1.2), win_sz, fill="red") # draw tick marks for i in range(1, num_ticks): tic_theta = (float(i))/num_ticks*360.0 if (fmod((i), 3) == 0): bonus = 0.03 else: bonus = 0.00 # create outer ticks start_x = win_sz + radius*(1.12+bonus)*cos(tic_theta/180*pi) start_y = win_sz + radius*(1.12+bonus)*sin(tic_theta/180*pi) end_x = win_sz + radius*1.10*cos(tic_theta/180*pi) end_y = win_sz + radius*1.10*sin(tic_theta/180*pi) w.create_line(start_x, start_y, end_x, end_y, fill="red") # create inner ticks start_x = win_sz + radius*0.90*cos(tic_theta/180*pi) start_y = win_sz + radius*0.90*sin(tic_theta/180*pi) end_x = win_sz + radius*(0.88-bonus)*cos(tic_theta/180*pi) end_y = win_sz + radius*(0.88-bonus)*sin(tic_theta/180*pi) w.create_line(start_x, start_y, end_x, end_y, fill="red") # draw polaris pol_x = win_sz + (pol_rad*cos(theta/180*pi)*radius) pol_y = win_sz - (pol_rad*sin(theta/180*pi)*radius) w.create_oval((pol_x-star_radius), (pol_y-star_radius), (pol_x+star_radius), (pol_y+star_radius), outline="white", fill="white") # add text w.create_text(0+star_radius, 0+star_radius, anchor="nw", fill="white", font=("courier",str(font_sz)), text=reporter_text) return # start the program! # declare the GUI root window root = Tk() root.title('Location of Polaris in NJP Finder - Viva La NJP!') # create a menu bar menu = Menu(root) root.config(menu=menu) # create the "File" menu filemenu = Menu(menu) menu.add_cascade(label="File", menu=filemenu) filemenu.add_command(label="Exit", command=root.quit) # add longitude input field, set default value longitude_label = Label(root, text="Local Longitude (+ East, - West): ") longitude_label.grid(row=1, column=0, sticky='nw') longitude_entry = Entry(root) longitude_entry.grid(row=1, column=1, sticky='nw') longitude_entry.insert(0, '-122.0650') # add refresh button below view window go_button = Button(root, text="Refresh View", command=update_figure) go_button.grid(row=1, column=2) # create the initial polaris view (at default longitude value) draw_figure() # start up the mainloop mainloop()