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)  
    return total  
  
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  
  
# Station file read  
sta_X, sta_Y, sta_Z = [], [], []  
  
stafile = 'sta_list.txt'  
with open(stafile) as data:  
    lines = data.readlines()  
  
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.

Delete this entry (admin only).