gad raifman
gad raifman

Reputation: 63

Solving an equation for several variables to create grid python3

after creating 4 sets of evenly spaced numbers (via linspace) and loading a set of coordinates and arrival times( relevant for the equation, there are 20 different coordinate sets ), i have created a four loops in order to get rms values for every set of variables (should have 24255 different sets), however as of now all i get is 20 rms values.

how would i get all the necessary values in a grid, that also shows the corresponding x,y,z,t0 values?

this is the code i used:

import numpy as np

#opening  files 

s_Coord = np.genfromtxt(fname = 's_Coord_Pset3.txt')

t_ob = np.genfromtxt(fname ='t_ob_Pset3File.txt')

#creating parameters

nt0 = 5
t0_vec = np.linspace(0,4,nt0)

nx = 21
x_vec = np.linspace(-20,20,nx)

ny = 21
y_vec = np.linspace(-20,20,ny)

nz = 11
z_vec = np.linspace(0,20,nz)

v = 6 #km

#calculating ti_p and rms
for x in x_vec:
    for y in y_vec:
        for z in z_vec:
            for t in t0_vec:
                rms_mat = np.zeros((nx,ny,nz,nt0))
                #the rms_mat zeros matrix is possibly for creating the grid to fill?
                g_range = np.arange(0,20,1) #for 20 stations from the file
                x1 = [s_Coord[g:g+1,0] for g in g_range]
                y1 = [s_Coord[g:g+1,1] for g in g_range]
                z1 = [s_Coord[g:g+1,2] for g in g_range]
                t0 = [t_ob[g:g+1] for g in g_range]
                ti_p = t+(((x1-x)**2+(y1-y)**2+(z1-z)**2)**0.5)/v
                rms = ((((t0-ti_p)**2)/(nx*ny*nz*nt0))**0.5)

any help with this would be much appreciated

as for the s_Coord file: this is sample the data:

   5.5096956e+00  -1.0756354e+01   0.0000000e+00
  -2.4259195e+00   1.4407220e+01   0.0000000e+00
   3.6707655e+00   5.3756022e+00   0.0000000e+00
   4.7193844e+00   9.2582846e+00   0.0000000e+00
  -1.8361331e+00  -1.3745789e+01   0.0000000e+00
     . . ........................................

as for the t_ob file: this is sample the data:

   6.0261696e+00
   5.1844492e+00
   4.7735687e+00
   5.1549808e+00
   6.0337739e+00

i tried appending the values in this way, but the code just got too heavy with a long run time without finishing running or just crashes(i believe there is a problem with how it was written

ti_p = []
rms_val = []
grid = []
for x in x_vec:
    for y in y_vec:
        for z in z_vec:
            for t in t0_vec:
                g_range = np.arange(0,20,1) #for 20 stations from the file
                x1 = [s_Coord[g:g+1,0] for g in g_range]
                y1 = [s_Coord[g:g+1,1] for g in g_range]
                z1 = [s_Coord[g:g+1,2] for g in g_range]

                ti_p.append(t+(((x1-x)**2+(y1-y)**2+(z1-z)**2)**0.5)/v)
                rms_val.append(((((t-ti_p)**2)/(nx*ny*nz*nt0))**0.5))
                grid.append((x,y,z,t,rms_val))

finally, here after the fix proposed below by Andrew H, this is the last block fixed to work

rms_mat = np.zeros((nx,ny,nz,nt0))
ti_p_mat =  np.zeros(rms_mat.shape , dtype=object) 
for i, x in enumerate(x_vec):
    for j, y in enumerate(y_vec):
        for k, z in enumerate(z_vec):
            for l, t in enumerate(t0_vec):
                g_range = np.arange(0,20,1) #for 20 stations from the file
                x1 = [s_Coord[g:g+1,0] for g in g_range]
                y1 = [s_Coord[g:g+1,1] for g in g_range]
                z1 = [s_Coord[g:g+1,2] for g in g_range]
                t0 = [t_ob[g:g+1] for g in g_range]
                ti_p = t+(((x1-x)**2+(y1-y)**2+(z1-z)**2)**0.5)/v
                rms = (((sum((t0-ti_p)**2))/(len(s_Coord)))**0.5)
                rms_mat[i, j, k, l] = rms
                ti_p_mat[i, j, k, l] = ti_p

Upvotes: 0

Views: 61

Answers (1)

Andrew Holmgren
Andrew Holmgren

Reputation: 1275

You are overwriting your rms array. You want to either preallocate rms and write the locally calculated array to a slice of rms, or you want to append.

You probably want something like this

rms_mat = np.zeros((nx,ny,nz,nt0))
ti_p_mat =  np.zeros(rms_mat.shape) 
for i, x in enumerate(x_vec):
    for j, y in enumerate(y_vec):
        for k, z in enumerate(z_vec):
            for l, t in enumerate(t0_vec):
                g_range = np.arange(0,20,1) #for 20 stations from the file
                x1 = [s_Coord[g:g+1,0] for g in g_range]
                y1 = [s_Coord[g:g+1,1] for g in g_range]
                z1 = [s_Coord[g:g+1,2] for g in g_range]
                t0 = [t_ob[g:g+1] for g in g_range]
                ti_p = t+(((x1-x)**2+(y1-y)**2+(z1-z)**2)**0.5)/v
                rms = ((((t0-ti_p)**2)/(nx*ny*nz*nt0))**0.5)
                rms_mat[i, j, k, l] = rms
                ti_pi_mat[i, j, k, l] = ti_p

Upvotes: 1

Related Questions