import numpy as np
import matplotlib.pyplot as plt

a = 5600
b = 2600
rho = 2200

def H(t):
    h = np.zeros(len(t))
    for i, it in enumerate(t):
        if it > 0:
            h[i]=1
        else:
            h[i]=0.0
    return h

def pwave(x,t,i,p):
    r = np.sqrt(x[0]**2 + x[1]**2 + x[2]**2)
    gi = x[i-1]/r
    gp = x[p-1]/r
    return gi*gp*H(t-r/a)*np.exp(-t)*np.exp(r/a)/(4.0*np.pi*a**2*r)

x = [100000.0, 100000.0, 100000.0]
t = np.linspace(0,100,128)

signal = pwave(x,t,1,3)

plt.plot(signal)
plt.show()