#!/usr/bin/env python
# -*- coding: latin-1 -*-
import os, sys
import math
import matplotlib.pyplot as plt
#import travelmethods.py

# set constants
R = 6371.0
# assume average P velocity
vp = 11.
# rpd is radians per degree
rpd = math.pi/180.
# kpd is km per degree
kpd = R*rpd
# set earthquake origin time:
Oh=4
Om=49
Os=20

# initialize lists:
distances = []
t0s=[]
res3 = [ ]
# opens the date file ptimescopies, which has the header deleted all but P #and np columns
ptm = open('Ptimes.txt','r')
for line in ptm:
    fields = line.split()
    if fields[6]=='P' or fields[6]=='Pn':
        time = fields[7].split('Z')[0].split('T')[1].split(':')
        seconds = 3600*(float(time[0])-Oh)+60*(float(time[1])-Om)+float(time[2])-Os
        distances = distances + [float(fields[4])]
        t0s = t0s + [seconds] 
        res3 = res3 + [float(fields[9])]
        #print(“done”)

# Measuring along earths surface chord
t1s = []
for d in distances:
   b = d*kpd
# Then convert to distances, d
   
   t1 =  b/vp
   #print(t1s)
   #print('%g  %g' % (d,t1))
   #plt.plot(distances,t1,'blue')
   #plt.scatter(d,t1)
   #plt.show()
   #plt.plot(t0s, distances, 'o', color='black');
   t1s = t1s + [t1]
#print (t1s)

t2s = []
for d in distances:
   th = d*rpd
   hth = 0.5*th
# convert radians k
   a = 2*R*math.sin(hth)
   t2 = a/vp
   #print('%g  %g' % (d,t2s))
   #plt.plot(d,t2s, 'o', color='black')
t2s = t2s + [t2]
#plt.plot(t1s, t2s)

t3s = [t0-r3 for t0,r3 in zip(t0s,res3)]
pairedlist = zip(t0s,t1s)
pairedlist2 = zip(t0s,t2s)
pairedlist3 = zip(t0s,t3s)
res1 = [t0-t1 for t0,t1 in pairedlist]
res2 = [t0-t2 for t0,t2 in pairedlist2]
res3 = [t0-t3 for t0,t3 in pairedlist3]

fig = plt.figure()
#fig.subplots_adjust(bottom=0.2)
#ax1 = fig.add_subplot(211)
#line1=ax1.plot(t1s,'bo',label='list 1')
#line2 = ax1.plot(t2s,'go-',label='list 2')
#line3 = ax1.plot(t3s,'ro-',label='list 3')

#plt.plot(t3s,)
#plt.xlabel("Distance (degrees)")
plt.ylabel("P wave travel time (seconds)")
#lines = line1+line2+line3
#labels = [l.get_label() for l in lines]
#plt.show()
plt.hist(res1, 100, range=(-200,200),alpha=0.7)
plt.hist(res2, 100, range=(-200,200),alpha=0.7)
plt.hist(res3, 100, range=(-200,200),alpha=0.7)
plt.xlabel("Distance (degrees)")
plt.ylabel("P wave travel time (seconds)")
plt.show()