# Created by Yuk Tung Liu
#! /usr/bin/env python

from math import sin, cos, atan2, pi, acos, floor

# RA and DEC of the GW source:
alpha_hours = 8.0
delta_degrees = -70.0

# degrees to radian
deg_to_rad = pi/180

alpha = alpha_hours*15*deg_to_rad
delta = delta_degrees*deg_to_rad

# LIGO detector geometry from 
# http://www.ligo.org/scientists/GW100916/GW100916-geometry.html
# and also 
# http://www.ligo.org/scientists/GRB051103/GRB051103-geometry.php
# The two pages are consistent.

# LIGO Hanford (LHO) detector geometry 
# Latitude, Longitude and azimuth of the x-arm
# Note: Longitude is measured positively westward from Greenwich
# Azimuth is measured from the South point, turning positive to the West.
lat_LHO = (46 + 27.0/60 + 19.0/3600)*deg_to_rad
long_LHO = (119 + 24.0/60 + 28.0/3600)*deg_to_rad 
xazi_LHO = (180-36)*deg_to_rad

# LIGO Livingston (LLO) detector geometry
# Latitude, Longitude and azimuth of the x-arm
lat_LLO = (30 + 33.0/60 + 46.0/3600)*deg_to_rad
long_LLO = (90 + 46.0/60 + 27.0/3600)*deg_to_rad
xazi_LLO = (90-18)*deg_to_rad

# Time of detection of GW150914: 2015-09-14 09:50:45 UTC
# The Greenwich sidereal time theta_G is calculated by 
# the formula at https://en.wikipedia.org/wiki/Sidereal_time
# The result is theta_G = 9:26:56

# Days between 2015-09-14 0:00 and 2000-01-01 0:00
D = 5736.0 
# Days between 2015-09-14 09:50:45 and 2000-01-01 12:00
D = D + (9 + 50.0/60 + 45.0/3600 - 12)/24 
theta_G = 18.697374558 + 24.06570982441908*D
# theta_G mod 24h 
theta_G = theta_G - floor(theta_G/24.0)*24
theta_G = theta_G*15*deg_to_rad

# Equatorial -> horizontal using the formulae at
# https://en.wikipedia.org/wiki/Celestial_coordinate_system

# Hour angles
h_LHO = theta_G - long_LHO - alpha
h_LLO = theta_G - long_LLO - alpha

# Zenith distance theta and azimuth A
theta_LHO = acos( sin(lat_LHO)*sin(delta) + cos(lat_LHO)*cos(delta)*cos(h_LHO) )
theta_LLO = acos( sin(lat_LLO)*sin(delta) + cos(lat_LLO)*cos(delta)*cos(h_LLO) )
y = cos(delta)*sin(h_LHO)
x = cos(delta)*cos(h_LHO)*sin(lat_LHO) - sin(delta)*cos(lat_LHO)
A_LHO = atan2(y,x)
phi_LHO = xazi_LHO - A_LHO
y = cos(delta)*sin(h_LLO)
x = cos(delta)*cos(h_LLO)*sin(lat_LLO) - sin(delta)*cos(lat_LLO)
A_LLO = atan2(y,x)
phi_LLO = xazi_LLO - A_LLO

cth_LHO = cos(theta_LHO)
s2ph_LHO = sin(2*phi_LHO)
c2ph_LHO = cos(2*phi_LHO)
cth_LLO = cos(theta_LLO)
s2ph_LLO = sin(2*phi_LLO)
c2ph_LLO = cos(2*phi_LLO)

Fp_LHO = 0.5*(1+cth_LHO*cth_LHO)*c2ph_LHO
Fx_LHO = cth_LHO*s2ph_LHO
Fp_LLO = 0.5*(1+cth_LLO*cth_LLO)*c2ph_LLO
Fx_LLO = cth_LLO*s2ph_LLO

print "GW source RA (hours) = ",alpha_hours
print "GW source DEC (degrees) = ",delta_degrees
print "Time: 2015-09-14 09:50:45 UTC"
print " "
print "LIGO Hanford (LHO) detector:" 
print "F_+ and F_x: ",Fp_LHO, Fx_LHO
print " "
print "LIGO Livingston (LLO) detector:"
print "F_+ and F_x: ",Fp_LLO, Fx_LLO