# How do I get the velocity of the P-wave or S-wave without using the travel time graph?

I'm trying to write a program that would determine the magnitude and epicenter of an earthquake. From http://www.dartmouth.edu/~ears5/handouts/EQ_distance.html, the distance of the epicenter can be calculated with

$$d = (ts-tp)/(1/vs)-(1/vp)$$

However, I'm not really sure as to what the values Vs and Vp are. From the travel time graph Vs and Vp varies with the distance of the sensor from the earthquake. Given that I have 3 sensors with given latitude and longitude, from which I am able to get the lag time, is there anyway I can estimate the Vs and Vp from these variables?

Otherwise, is there any other way I could get the distance of the sensors from the epicenter?

I think that it is possible with fours stations yes, not with three.

With the assumption of (roughly) homogeneous velocities we get the following forward problem. Assuming that $$v_p$$ and $$v_s$$ are constants between the Earthquake and all your receivers yields the following system: \begin{align} d &= \sqrt{\Delta x^2 + \Delta y^2 + \Delta z^2}=\text{3 knowns, 3 unknowns},\quad \text{for} \left\{ \begin{array}{rl} \Delta x &= x_\text{earthquake}-x_\text{receiver} \\ \Delta y & = y_\text{earthquake} - y_\text{receiver} \\ \Delta z & = z_\text{earthquake} - z_\text{receiver} \end{array} \begin{split} \end{split}\right.,\\ \Delta t & = t_s - t_p = \text{measured},\\ p &= \left(\frac{1}{v_s} - \frac{1}{v_p} \right) = \text{1 unknown}. \end{align} Restating the equation of your question then gives the forward problem: $$\Delta t = p d.$$ For example, for $$v_p=2000$$ m/s and $$v_s=1200$$ m/s, we find $$p=1/3000$$ and can compute e.g. with $$d = 600$$ that $$\Delta t = 600/3000 =0.2$$ s.

Now, we must set=up the inverse problem. You generally need at least as many equations as unknowns. In this case we have 4 unknowns (the velocity and the $$(x,y,z)$$ earthquake location). That means that we need 4 (unique) equations. I think that the simplest way of obtaining 4 unique equations is to combine 4 measurements: $$p\begin{pmatrix} d_1 \\ d_2 \\ d_3 \\ d_4 \end{pmatrix} = \begin{pmatrix} \Delta t_1 \\ \Delta t_2 \\ \Delta t_3 \\ \Delta t_4 \end{pmatrix},$$ Denoting with $$x_{1,2,3,4}$$ and $$y_{1,2,3,4}$$ the latitude and longitude of the receivers, and assuming measurements at the surface, we obtain: $$p\begin{pmatrix} \sqrt{(x_\text{earthquake}-x_1)^2 + (y_\text{earthquake}-y_1)^2 + z_\text{earthquake}^2} \\ \sqrt{(x_\text{earthquake}-x_2)^2 + (y_\text{earthquake}-y_2)^2 + z_\text{earthquake}^2} \\ \sqrt{(x_\text{earthquake}-x_3)^2 + (y_\text{earthquake}-y_3)^2 + z_\text{earthquake}^2 } \\ \sqrt{(x_\text{earthquake}-x_4)^2 + (y_\text{earthquake}-y_4)^2 + z_\text{earthquake}^2} \end{pmatrix} = \begin{pmatrix} \Delta t_1 \\ \Delta t_2 \\ \Delta t_3 \\ \Delta t_4 \end{pmatrix}.$$

Finally, we must solve the inverse problem. Unfortunately, the system is non-linear (the unknowns are 'hidden' in the root of the square and can't be separated). You must therefore use some non-linear solver. For example, in MATLAB you can do this

% True earthquake location
x_earthquake=0;
y_earthquake=0;
z_earthquake=1000;

% True station locations
x_1 = 1000; y_1 = 1000;
x_2 =  500; y_2 = -300;
x_3 = -400; y_3 = -100;
x_4 =  -10; y_4 = 800;
X = [x_1;x_2;x_3;x_4];
Y = [y_1;y_2;y_3;y_4];

% True velocity structure
p0 = 1/3000;

% Forward equation (F(1)=x_eartquake, F(2)=y_earthquake, F(3)=z_earthquake, F(4)=p.
t =@(F) F(4) * sqrt( (F(1)-X).^2 + (F(2)-Y).^2 + F(3).^2);

% True recordings
measured_times = t([x_earthquake,y_earthquake,z_earthquake,p0]);

% Misfit function (=0 at optimum)
G =@(F) t(F) - measured_times;

% Invert
%options = optimoptions('fsolve','FiniteDifferenceType','central'); % using these will give better results in MATLAB
%F_inv = fsolve( G, [1000,1000,1000,0.1],options)
F_inv = fsolve( G, [1000,1000,1000,0.1])
x_earthquake_inv = F_inv(1)
y_earthquake_inv = F_inv(2)
z_earthquake_inv = F_inv(3)
p_inv            = F_inv(4)


If you copy this into https://octave-online.net/, for example, it will find an earthquake at (0,0,1000) and p=1/3000, with some small errors due to the non-linear solver.