# Demo entry 6786635

hypo_determine

Submitted by anonymous on Apr 01, 2019 at 13:20
Language: Python 3. Code size: 3.2 kB.

```import numpy as np
from scipy.sparse.linalg import cg

### Function
def get_sq_sum(list_val):
total = 0
for i in list_val:
total = total + pow(i, 2)

def get_residual(ori_x, ori_y, ori_z, ori_ot, x, y, z, ot, sta_x, sta_y, sta_z, vel):
d_obs, d_cal, residual = [], [], []
for i in range(0, 10, 1):
d_obs_t = (1 / vel) * np.sqrt(pow(ori_x - sta_x[i], 2) + pow(ori_y - sta_y[i], 2) + pow(ori_z - sta_z[i], 2)) + ori_ot
d_cal_t = (1 / vel) * np.sqrt(pow(x - sta_x[i], 2) + pow(y - sta_y[i], 2) + pow(z - sta_z[i], 2)) + ot
d_obs.append(d_obs_t)
d_cal.append(d_cal_t)

residual_t = d_obs[i] - d_cal[i]
residual.append(residual_t)
return d_obs, d_cal, residual

def get_jaco_mat(x, y, z, sta_x, sta_y, vel):
g_mat = []
for i in range(0, 10, 1):
g1 = ((x - sta_x[i]) / vel) * np.power(np.power(x - sta_x[i], 2) + np.power(y - sta_y[i], 2) + np.power(z, 2), -1/2)
g2 = ((y - sta_y[i]) / vel) * np.power(np.power(x - sta_x[i], 2) + np.power(y - sta_y[i], 2) + np.power(z, 2), -1/2)
g3 = (z / vel) * np.power(np.power(x - sta_x[i], 2) + np.power(y - sta_y[i], 2) + np.power(z, 2), -1/2)
g4 = 1

g_mat_t = [g1, g2, g3, g4]
g_mat.append(g_mat_t)
return g_mat

sta_X, sta_Y, sta_Z = [], [], []

stafile = 'sta_list.txt'
with open(stafile) as data:

for line in lines:
sta_X.append(float(line.split()[0]))
sta_Y.append(float(line.split()[1]))
sta_Z.append(float(line.split()[2]))

ori_X, ori_Y, ori_Z, ori_T = 0.0, 0.0, 10.0, 0
X, Y, Z, OT = 3.0, 4.0, 20.0, 2
vel = 5

residual = get_residual(ori_X, ori_Y, ori_Z, ori_T, X, Y, Z, OT, sta_X, sta_Y, sta_Z, vel)[2]
sq_sum_resi = get_sq_sum(residual)

print ("\nX: {0:.2f},    Y: {1:.2f},    Z: {2:.2f},    OT: {3:.2f}".format(X, Y, Z, OT))
for i in range(0, 10, 1):
print ("Residual of station {0:02d} : {1:.7f}".format(i+1, residual[i]))
print ("-> Residual sum of squares : %6.2f" %sq_sum_resi)
print ("===================================================\n")

while True:
G_Mat = get_jaco_mat(X, Y, Z, sta_X, sta_Y, vel)

G_Mat_T = np.transpose(G_Mat)

GTG = G_Mat_T.dot(G_Mat)
GTd = G_Mat_T.dot(residual)

dm = cg(GTG, GTd)

X = X + dm[0][0]
Y = Y + dm[0][1]
Z = Z + dm[0][2]
OT = OT + dm[0][3]

residual = get_residual(ori_X, ori_Y, ori_Z, ori_T, X, Y, Z, OT, sta_X, sta_Y, sta_Z, vel)[2]
sq_sum_resi = get_sq_sum(residual)

print("X: {0:.2f},    Y: {1:.2f},    Z: {2:.2f},    OT: {3:.2f}".format(X, Y, Z, OT))
for i in range(0, 10, 1):
print("Residual of station {0:02d} : {1:.7f}".format(i + 1, residual[i]))
print("-> Residual sum of squares : %6.2f" % sq_sum_resi)
print("===================================================\n")

if sq_sum_resi < 0.1:
break
```

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.