extremely simplified eq location

surface-wave-eq-location.py

import numpy as np
import matplotlib.pyplot as plt
import datetime
import geopy.distance

T0 = datetime.datetime(2023,9,16,12,31,0).timestamp() # test origin time
v = 1 # surface wave speed

def G(lon,lat,m,vv):
    #m = (x0,y0,t0) is the source parameter vector
#     return np.sqrt( (xx-m[0])**2 + (yy-m[1])**2 ) / vv
    out = []
    for lo,la in zip(lon,lat):
        this = geopy.distance.geodesic( (la,lo) , (m[1],m[0]) ).km
        out.append(this)
    return out

# predict travel times at sites: IU.KBS,IU.KONO,IU.SFJD,NS.JNW,DK.SCO
# (the following vectors are arranged in that order)
lon_obs = np.array([11.9385,9.5946,-50.62076,-8.42367,-21.949699])
lat_obs = np.array([78.9154,59.6521,66.9961,71.02734,70.485603])
tobs = np.array([datetime.datetime(2023,9,16,12,40,0).timestamp(),
                 datetime.datetime(2023,9,16,12,45,0).timestamp(),
                 datetime.datetime(2023,9,16,12,40,0).timestamp(),
                 datetime.datetime(2023,9,16,12,38,0).timestamp(),
                 datetime.datetime(2023,9,16,12,36,0).timestamp()])

# now loop over test event locations
nlon=20
nlat=20
lons = np.linspace(-30,-20,nlon)
lats = np.linspace(70,82,nlat)
residuals = np.zeros((nlon,nlat))
for i,lon_test in enumerate(lons):
    for j,lat_test in enumerate(lats):
        tt = G(lon_obs,lat_obs,[lon_test,lat_test],v)
        residuals[i,j] = np.mean(np.abs(tt - (tobs - T0)))
        
fig,ax=plt.subplots()
c=plt.contourf(lons,lats,residuals)
plt.colorbar(c,label='Misfit (seconds)')
ax.axis('equal')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()