def test_p_iasp91_geo_fallback_manual(self):
        """
        Manual test for P phase in IASP91 given geographical input.

        This version of the test checks that things still work when
        geographiclib is not installed.
        """
        has_geographiclib_real = geodetics.HAS_GEOGRAPHICLIB
        geodetics.HAS_GEOGRAPHICLIB = False
        m = TauPyModel(model="iasp91")
        arrivals = m.get_travel_times_geo(source_depth_in_km=10.0,
                                          source_latitude_in_deg=20.0,
                                          source_longitude_in_deg=33.0,
                                          receiver_latitude_in_deg=55.0,
                                          receiver_longitude_in_deg=33.0,
                                          phase_list=["P"])
        geodetics.HAS_GEOGRAPHICLIB = has_geographiclib_real
        self.assertEqual(len(arrivals), 1)
        p_arrival = arrivals[0]

        self.assertEqual(p_arrival.name, "P")
        self.assertAlmostEqual(p_arrival.time, 412.43, 2)
        self.assertAlmostEqual(p_arrival.ray_param_sec_degree, 8.613, 3)
        self.assertAlmostEqual(p_arrival.takeoff_angle, 26.74, 2)
        self.assertAlmostEqual(p_arrival.incident_angle, 26.70, 2)
        self.assertAlmostEqual(p_arrival.purist_distance, 35.00, 2)
        self.assertEqual(p_arrival.purist_name, "P")

    def test_pierce_p_iasp91_geo(self):
        """
        Test single pierce point against output from TauP using geo data.

        This version of the test is used when geographiclib is installed
        """
        m = TauPyModel(model="iasp91")
        arrivals = m.get_pierce_points_geo(source_depth_in_km=10.0,
                                           source_latitude_in_deg=-45.0,
                                           source_longitude_in_deg=-50.0,
                                           receiver_latitude_in_deg=-80.0,
                                           receiver_longitude_in_deg=-50.0,
                                           phase_list=["P"])
        self.assertEqual(len(arrivals), 1)
        p_arr = arrivals[0]

        # Open test file.
        filename = os.path.join(DATA,
                                "taup_pierce_-mod_isp91_ph_P_-h_10_-evt_" +
                                "-45_-50_-sta_-80_-50")

        expected = np.genfromtxt(filename, skip_header=1)

        np.testing.assert_almost_equal(expected[:, 0],
                                       np.degrees(p_arr.pierce['dist']), 2)
        np.testing.assert_almost_equal(expected[:, 1],
                                       p_arr.pierce['depth'], 1)
        np.testing.assert_almost_equal(expected[:, 2],
                                       p_arr.pierce['time'], 1)
        np.testing.assert_almost_equal(expected[:, 3],
                                       p_arr.pierce['lat'], 1)
        np.testing.assert_almost_equal(expected[:, 4],
                                       p_arr.pierce['lon'], 1)

    def test_regional_models(self):
        """
        Tests small regional models as this used to not work.

        Note: It looks like too much work to get a 1-layer model working.
        The problem is first in finding the moho, and second in coarsely-
        sampling slowness. Also, why bother.
        """
        model_names = ["2_layer_model", "5_layer_model"]
        expected_results = [
            [("p", 18.143), ("Pn", 19.202), ("PcP", 19.884), ("sP", 22.054),
             ("ScP", 23.029), ("PcS", 26.410), ("s", 31.509), ("Sn", 33.395),
             ("ScS", 34.533)],
            [("Pn", 17.358), ("P", 17.666), ("p", 17.804), ("P", 17.869),
             ("PcP", 18.039), ("ScP", 19.988), ("sP", 22.640), ("sP", 22.716),
             ("sP", 22.992), ("PcS", 23.051), ("sP", 24.039), ("sP", 24.042),
             ("Sn", 30.029), ("S", 30.563), ("s", 30.801), ("S", 30.913),
             ("ScS", 31.208)]]

        for model_name, expects in zip(model_names, expected_results):
            with TemporaryWorkingDirectory():
                folder = os.path.abspath(os.curdir)
                build_taup_model(
                    filename=os.path.join(DATA, os.path.pardir,
                                          model_name + ".tvel"),
                    output_folder=folder, verbose=False)
                model = TauPyModel(os.path.join(folder, model_name + ".npz"))

            arrvials = model.get_ray_paths(source_depth_in_km=18.0,
                                           distance_in_degree=1.0)

            self.assertEqual(len(arrvials), len(expects))
            for arrival, expect in zip(arrvials, expects):
                self.assertEqual(arrival.name, expect[0])
                self.assertAlmostEqual(arrival.time, expect[1], 3)

def get_first_arrival(st, model='ak135'):
    '''
    Returns first arrival information for a particular stream and a theoretical velocity model ak135 or iasp91.
    
    Output is phase name, arrival time (in s after origin) and slowness (in s / km).
        
    Parameters
    ----------
    st : ObsPy Stream object
        Stream of SAC format seismograms for the seismic array
    model : string
        Model to use for the travel times, either 'ak135' or 'iasp91'.

    Returns
    -------
    phase : Phase object
        Phase object containing the phase name, arrival time, and slowness of the first arrival
    '''

    # Read event depth and great circle distance from SAC header
    depth = st[0].stats.sac.evdp
    delta = st[0].stats.sac.gcarc

    taup = TauPyModel(model)
    first_arrival = taup.get_travel_times(depth, delta)[0]

    phase = Phase(first_arrival.name, first_arrival.time, first_arrival.ray_param_sec_degree / G_KM_DEG)

    return phase

def get_arrivals(st, model='ak135'):
    '''
    Returns complete arrival information for a particular stream and a theoretical velocity model ak135 or iasp91.
    
    Output is phase name, arrival time (in s after origin) and slowness (in s / km).
        
    Parameters
    ----------
    st : ObsPy Stream object
        Stream of SAC format seismograms for the seismic array
    model : string
        Model to use for the travel times, either 'ak135' or 'iasp91'.

    Returns
    -------
    phase_list : list
        List containing Phase objects containing the phase name, arrival time, and slowness of each arrival
    '''

    # Read event depth and great circle distance from SAC header
    depth = st[0].stats.sac.evdp
    delta = st[0].stats.sac.gcarc

    tau_model = TauPyModel(model)
    arrivals = tau_model.get_travel_times(depth, delta)
    
    phase_list = []
    
    for arrival in arrivals:
        phase = Phase(arrival.name, arrival.time, arrival.ray_param_sec_degree / G_KM_DEG)
        phase_list.append(phase)

    return phase_list

    def test_surface_wave_ttimes(self):
        """
        Tests the calculation of surface ttimes.

        Tested against a reference output from the Java TauP version.
        """
        for model, table in [("iasp91", "iasp91_surface_waves_table.txt"),
                             ("ak135", "ak135_surface_waves_table.txt")]:
            m = TauPyModel(model=model)
            filename = os.path.join(DATA, table)

            with open(filename, "rt") as fh:
                for line in fh:
                    _, distance, depth, phase, time, ray_param, _, _ \
                        = line.split()
                    distance, depth, time, ray_param = \
                        map(float, [distance, depth, time, ray_param])

                    arrivals = m.get_travel_times(
                        source_depth_in_km=depth, distance_in_degree=distance,
                        phase_list=[phase])

                    self.assertTrue(len(arrivals) > 0)

                    # Potentially multiple arrivals. Get the one closest in
                    # time and closest in ray parameter.
                    arrivals = sorted(
                        arrivals,
                        key=lambda x: (abs(x.time - time),
                                       abs(x.ray_param_sec_degree -
                                           ray_param)))
                    arrival = arrivals[0]
                    self.assertEqual(round(arrival.time, 2), round(time, 2))
                    self.assertEqual(round(arrival.ray_param_sec_degree, 2),
                                     round(ray_param, 2))

    def test_different_models(self):
        """
        Open all included models and make sure that they can produce
        reasonable travel times.
        """
        models = ["1066a", "1066b", "ak135", "herrin", "iasp91", "prem",
                  "sp6", "jb", "pwdk", "ak135f_no_mud"]
        for model in models:
            m = TauPyModel(model=model)

            # Get a p phase.
            arrivals = m.get_travel_times(
                source_depth_in_km=10.0, distance_in_degree=50.0,
                phase_list=["P"])
            # AK135 travel time.
            expected = 534.4
            self.assertTrue(abs(arrivals[0].time - expected) < 5)

            # Get an s phase.
            arrivals = m.get_travel_times(
                source_depth_in_km=10.0, distance_in_degree=50.0,
                phase_list=["S"])
            # AK135 travel time.
            expected = 965.1
            # Some models do produce s-waves but they are very far from the
            # AK135 value.
            self.assertTrue(abs(arrivals[0].time - expected) < 50)

def main():
                                      #           MAIN PROGRAM BODY
    OT,stlat,stlon,evlat,evlon,depth = getoptions()
    origin_time = UTCDateTime(str(OT))
    result = client.distaz(stalat=stlat, stalon=stlon, evtlat=evlat,evtlon=evlon)
    model = TauPyModel(model="AK135")
    arrivals = model.get_travel_times(source_depth_in_km=depth,distance_in_degree=result['distance'])#,
                                     #phase_list = ['P','PcP','PP','PKiKP','S','SS','ScS','SKiKS'])
    print "Distance = {0:.1f} arc degrees.".format(result['distance'])
    print "{0:.0f} Km distance.".format(result['distance']*111.25)
    print "{0:.0f} deg back Azimuth.".format(result['backazimuth'])
    table = client.traveltime(evloc=(evlat,evlon),staloc=[(stlat,stlon)],evdepth=depth)
    print "Selected phase list:\n"
    print (table.decode())

                   #   Print the phases, travel time and forecasted arrival time.
    phasename = []
    phasetime = []
    arrivaltime = []
    print "For origin time {}, ".format(origin_time)
    print "TauP big list of phases and arrival times:"
    for i in range(0,len(arrivals)):
        phasename.append(arrivals[i].name)
        phasetime.append(arrivals[i].time)
        at = origin_time+(arrivals[i].time)
        arrivaltime.append(at)
        print 'Phase: {0} \t arrives in {1:.2f} sec. at time {2:02.0f}:{3:02.0f}:{4:02.0f}.{5:02.0f}' \
              .format(arrivals[i].name,arrivals[i].time,at.hour,at.minute,at.second,at.microsecond/10000)
    arrivalpaths = model.get_ray_paths(source_depth_in_km=depth,distance_in_degree=result['distance'])#,\
                                #        phase_list = ['P','PcP','PP','PKiKP','S','SS','ScS','SKiKS'])
    arrivalpaths.plot()

def get_station_delays(station_coords,sources,velmod='PREM',phase_list=['s','S']):
    '''
    Given an array of station corodiantes and sources calculate the travel time
    from each source toe ach site.
    
    velmod is the FULL path to the .npz file used by TauPy
    '''
    
    from obspy.taup import TauPyModel
    from numpy import zeros
    from obspy.geodetics.base import locations2degrees

    model = TauPyModel(model=velmod)
    
    #Initalize output variabe
    delay_time=zeros((len(sources),len(station_coords)))
    
    #Loop over sources
    for ksource in range(len(sources)):
        
        print('Working on source %d of %d ' % (ksource,len(sources)))
        
        #loop over sites
        for ksite in range(len(station_coords)):
        
            distance_in_degrees=locations2degrees(station_coords[ksite,1],station_coords[ksite,0],
                                                      sources[ksource,2],sources[ksource,1])
            
            arrivals = model.get_travel_times(source_depth_in_km=sources[ksource,3],
                                  distance_in_degree=distance_in_degrees,phase_list=phase_list)
            
            delay_time[ksource,ksite]=arrivals[0].time
            
    return delay_time

    def test_single_path_ak135(self):
        """
        Test the raypath for a single phase. This time for model AK135.
        """
        filename = os.path.join(
            DATA, "taup_path_-o_stdout_-h_10_-ph_P_-deg_35_-mod_ak135")
        expected = np.genfromtxt(filename, comments='>')

        m = TauPyModel(model="ak135")
        arrivals = m.get_ray_paths(source_depth_in_km=10.0,
                                   distance_in_degree=35.0, phase_list=["P"])
        self.assertEqual(len(arrivals), 1)

        # Interpolate both paths to 100 samples and make sure they are
        # approximately equal.
        sample_points = np.linspace(0, 35, 100)

        interpolated_expected = np.interp(
            sample_points,
            expected[:, 0],
            expected[:, 1])

        interpolated_actual = np.interp(
            sample_points,
            np.round(np.degrees(arrivals[0].path['dist']), 2),
            np.round(6371 - arrivals[0].path['depth'], 2))

        self.assertTrue(np.allclose(interpolated_actual,
                                    interpolated_expected, rtol=1E-4, atol=0))

    def test_single_path_geo_iasp91(self):
        """
        Test the raypath for a single phase given geographical input.

        This tests the case when geographiclib is installed.
        """
        filename = os.path.join(DATA,
                                "taup_path_-mod_iasp91_-o_stdout_-h_10_" +
                                "-ph_P_-sta_-45_-60_evt_-80_-60")
        expected = np.genfromtxt(filename, comments='>')

        m = TauPyModel(model="iasp91")
        arrivals = m.get_ray_paths_geo(source_depth_in_km=10.0,
                                       source_latitude_in_deg=-80.0,
                                       source_longitude_in_deg=-60.0,
                                       receiver_latitude_in_deg=-45.0,
                                       receiver_longitude_in_deg=-60.0,
                                       phase_list=["P"])
        self.assertEqual(len(arrivals), 1)

        # Interpolate both paths to 100 samples and make sure they are
        # approximately equal.
        sample_points = np.linspace(0, 35, 100)

        interpolated_expected_depth = np.interp(
            sample_points,
            expected[:, 0],
            expected[:, 1])
        interpolated_expected_lat = np.interp(
            sample_points,
            expected[:, 0],
            expected[:, 2])
        interpolated_expected_lon = np.interp(
            sample_points,
            expected[:, 0],
            expected[:, 3])

        interpolated_actual_depth = np.interp(
            sample_points,
            np.round(np.degrees(arrivals[0].path['dist']), 2),
            np.round(6371 - arrivals[0].path['depth'], 2))
        interpolated_actual_lat = np.interp(
            sample_points,
            np.round(np.degrees(arrivals[0].path['dist']), 2),
            np.round(arrivals[0].path['lat'], 2))
        interpolated_actual_lon = np.interp(
            sample_points,
            np.round(np.degrees(arrivals[0].path['dist']), 2),
            np.round(arrivals[0].path['lon'], 2))

        np.testing.assert_allclose(interpolated_actual_depth,
                                   interpolated_expected_depth,
                                   rtol=1E-4, atol=0)
        np.testing.assert_allclose(interpolated_actual_lat,
                                   interpolated_expected_lat,
                                   rtol=1E-4, atol=0)
        np.testing.assert_allclose(interpolated_actual_lon,
                                   interpolated_expected_lon,
                                   rtol=1E-4, atol=0)

   def migrate(self,plot=False):
      import geopy
      from geopy.distance import VincentyDistance

      '''
      This is a rewritten function that combines the functions find_pierce_coor and
      migrate_1d so that it's more efficient.  Still in testing stages. RM 2/6/16
      '''

      depth_range = np.arange(50,800,5)        #set range of pierce points
      value       = np.zeros((len(depth_range)))

      #geodetic info
      bearing     = self.az
      lon_s = self.ses3d_seismogram.sy
      lat_s = 90.0-self.ses3d_seismogram.sx
      lon_r = self.ses3d_seismogram.ry
      lat_r = 90.0-self.ses3d_seismogram.rx
      origin      = geopy.Point(lat_s, lon_s)

      #find how far away the pierce point is
      model  = TauPyModel(model='pyrolite_5km')

      for i in range(0,len(depth_range)):
         phase = 'P'+str(depth_range[i])+'s'
         pierce = model.get_pierce_points(self.eq_depth,self.delta_deg,phase_list=[phase])
         tt     = model.get_travel_times(self.eq_depth,self.delta_deg,phase_list=['P',phase])

         #in case there's duplicate phase arrivals
         for j in range(0,len(tt)):
            if tt[j].name == 'P':
               p_arr = tt[j].time           
            elif tt[j].name == phase:
               phase_arr = tt[j].time

         #determine value 
         Pds_time = phase_arr - p_arr
         i_start  = int((0.0 - self.window_start)/self.ses3d_seismogram.dt)
         i_t      = int(Pds_time/self.ses3d_seismogram.dt) + i_start
         value[i] = self.prf[i_t]

         points = pierce[0].pierce
         for j in range(0,len(points)):
            if points[j]['depth'] == depth_range[i] and points[j]['dist']*(180.0/np.pi) > 20.0:
               prc_dist = points[j]['dist']*(180.0/np.pi)
               d_km = prc_dist * ((2*np.pi*6371.0/360.0))
               destination = VincentyDistance(kilometers=d_km).destination(origin,bearing)
               lat = destination[0]
               lon = destination[1]
               row = {'depth':depth_range[i],'dist':prc_dist,'lat':lat,'lon':lon,'value':value[i]}
               self.pierce_dict.append(row)

      if plot == True:
         plt.plot(value,depth_range)
         plt.gca().invert_yaxis()
         plt.show()

      return value,depth_range

def get_taupy_points(
    center_lat, center_lon, ev_lat, ev_lon, ev_depth, stime, etime, mini, maxi, ev_otime, phase_shift, sll, slm
):

    distance = locations2degrees(center_lat, center_lon, ev_lat, ev_lon)
    # print(distance)

    model = TauPyModel(model="ak135")
    arrivals = model.get_pierce_points(ev_depth, distance)
    # arrivals = earthmodel.get_pierce_points(ev_depth,distance,phase_list=('PP','P^410P'))

    # compute the vespagram window
    start_vespa = stime - mini
    end_vespa = etime - maxi

    # compare the arrival times with the time window
    count = 0
    k = 0
    phase_name_info = []
    phase_slowness_info = []
    phase_time_info = []

    for i_elem in arrivals:
        # print(i_elem)
        dummy_phase = arrivals[count]
        # print(dummy_phase)
        # phase time in seconds
        taup_phase_time = dummy_phase.time
        # print(taup_phase_time)
        # slowness of the phase
        taup_phase_slowness = dummy_phase.ray_param_sec_degree
        # compute the UTC travel phase time
        taup_phase_time2 = ev_otime + taup_phase_time + phase_shift

        # print(start_vespa)
        # print(end_vespa)
        # print(taup_phase_time2)

        if start_vespa <= taup_phase_time2 <= end_vespa:  # time window
            if sll <= taup_phase_slowness <= slm:  # slowness window

                # seconds inside the vespagram
                taup_mark = taup_phase_time2 - start_vespa
                # store the information
                phase_name_info.append(dummy_phase.name)
                phase_slowness_info.append(dummy_phase.ray_param_sec_degree)
                phase_time_info.append(taup_mark)
                # print(phases_info[k])
                k += 1

        count += 1

    # print(phase_name_info)

    phase_slowness_info = np.array(phase_slowness_info)
    phase_time_info = np.array(phase_time_info)

    return phase_name_info, phase_slowness_info, phase_time_info

 def test_p_ak135(self):
     """
     Test P phase arrival against TauP output in in model AK135.
     """
     m = TauPyModel(model="ak135")
     arrivals = m.get_travel_times(source_depth_in_km=10.0,
                                   distance_in_degree=35.0,
                                   phase_list=["P"])
     self._compare_arrivals_with_file(
         arrivals, "taup_time_-h_10_-ph_ttall_-deg_35_-mod_ak135")

    def test_buried_receiver(self):
        """
        Simple test for a buried receiver.
        """
        m = TauPyModel(model="iasp91")
        arrivals = m.get_travel_times(
            source_depth_in_km=10.0, distance_in_degree=90.0,
            receiver_depth_in_km=50,
            phase_list=["P", "PP", "S"])

        self._compare_arrivals_with_file(arrivals, "buried_receivers.txt")

 def test_ak135(self):
     """
     Test travel times for lots of phases against output from TauP in model
     AK135.
     """
     m = TauPyModel(model="ak135")
     arrivals = m.get_travel_times(source_depth_in_km=10.0,
                                   distance_in_degree=35.0,
                                   phase_list=["ttall"])
     self._compare_arrivals_with_file(
         arrivals, "taup_time_-h_10_-ph_ttall_-deg_35_-mod_ak135")

 def test_p_pwdk(self):
     """
     Test P phase arrival against TauP output in model pwdk
     with different cache values to test `TauModel.load_from_depth_cache`
     """
     for cache in self.caches:
         m = TauPyModel(model="pwdk", cache=cache)
         arrivals = m.get_travel_times(source_depth_in_km=10.0,
                                       distance_in_degree=35.0,
                                       phase_list=["P"])
         self._compare_arrivals_with_file(
             arrivals, "taup_time_-h_10_-ph_P_-deg_35_-mod_pwdk")

 def test_ppointvsobspytaup_S2P(self):
     slowness = 12.33
     evdep = 12.4
     evdist = 67.7
     pp1 = self.model.ppoint_distance(200, slowness, phase="P")
     model = TauPyModel(model="iasp91")
     arrivals = model.get_ray_paths(evdep, evdist, ("S250p",))
     arrival = arrivals[0]
     index = np.searchsorted(arrival.path["depth"][::-1], 200)
     pdist = arrival.path["dist"]
     pp2 = degrees2kilometers((pdist[-1] - pdist[-index - 1]) * 180 / np.pi)
     self.assertLess(abs(pp1 - pp2) / pp2, 0.2)

def zero_phase(tr,phase,rf_window=[-10,120],**kwargs):
##############################################################################
   '''
   Finds predicted PP arrival and zeros a window centered on the arrival
   args--------------------------------------------------------------------------
   tr: obspy trace
   phase: phase to zero
   time_start: starting point of receiver function time window

   kwargs---------------------------------------------------------------------
   window_half_dur : half duration of zero window (default = 2.5 s)
   taup_model
   '''
   window_half_dur = kwargs.get('window_half_dur',2.5)
   taup_model = kwargs.get('taup_model','none')
   if taup_model == 'none':
      taup_model = TauPyModel('prem')

   arrs = taup_model.get_travel_times(source_depth_in_km=tr.stats.evdp,
                                      distance_in_degree=tr.stats.gcarc,
                                      phase_list=['P',phase])
   print 'arrs = ', arrs
   P_arr  = 'none'
   phase_arr = 'none'

   for arr in arrs:
      if arr.name == 'P':
         P_arr = arr
      elif arr.name == phase:
         phase_arr = arr

   if P_arr == 'none' or phase_arr == 'none':
      raise ValueError('problem occured in function "zero_phase", no matching arrivals found')
   else:
      P_time = P_arr.time
      phase_time = phase_arr.time
      delay_time = phase_time - P_time
      zero_window_center = -1.0*rf_window[0] + delay_time
      zero_window_start  = zero_window_center - window_half_dur
      zero_window_end    = zero_window_center + window_half_dur
      zero_window_start_index = int(zero_window_start/tr.stats.delta)
      zero_window_end_index   = int(zero_window_end/tr.stats.delta)

      #case 1: entire window is in range
      if zero_window_start_index >= 0 and zero_window_end_index <= len(tr.data):
         tr.data[zero_window_start_index:zero_window_end_index] = 0.0
      #case 2: end of window is out of range
      if zero_window_start_index >= 0 and zero_window_end_index >= len(tr.data):
         tr.data[zero_window_start_index:] = 0.0
      #case 3: entire window is out of range
      if zero_window_start_index >= len(tr.data):
         print "PP arrives outside the receiver function window"

 def test_ak135(self):
     """
     Test travel times for lots of phases against output from TauP in model
     AK135 with different cache values to test
     `TauModel.load_from_depth_cache`
     """
     for cache in self.caches:
         m = TauPyModel(model="ak135", cache=cache)
         arrivals = m.get_travel_times(source_depth_in_km=10.0,
                                       distance_in_degree=35.0,
                                       phase_list=["ttall"])
         self._compare_arrivals_with_file(
             arrivals, "taup_time_-h_10_-ph_ttall_-deg_35_-mod_ak135")