From 8f5b9d313861b69546497e33cc3896dff47ee322 Mon Sep 17 00:00:00 2001 From: schipp Date: Thu, 20 Aug 2020 11:38:23 +0200 Subject: [PATCH] fix --- logic.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/logic.py b/logic.py index 169fb9a..50730b8 100644 --- a/logic.py +++ b/logic.py @@ -39,10 +39,10 @@ def get_csdm(spectra:np.ndarray) -> np.ndarray: """ # einstein convention version - csdm = np.einsum('ik,jk->ijk', source_spectra, np.conj(source_spectra)) + csdm = np.einsum('ik,jk->ijk', spectra, np.conj(spectra)) # # old slow version - # csdm = np.zeros((n_stations, n_stations, len(freqs[freqs_of_interest])), dtype=np.complex) + # csdm = np.zeros((n_stations, n_stations, len(freqs[freqs_of_interest_idx])), dtype=np.complex) # for i, spec1 in enumerate(source_spectra): # for j, spec2 in enumerate(source_spectra): # csdm[i, j] = spec1 * np.conj(spec2) @@ -83,7 +83,7 @@ def bartlett_processor(csdm, gf_spectra): return beampower -def get_beampowers(csdm, gf_spectra, gp_dists): +def get_beampowers(csdm, gf_spectra, gp_dists, relevant_dists): """ Computes the Beampower for all grid-points. """ @@ -113,7 +113,7 @@ def get_beampowers(csdm, gf_spectra, gp_dists): voxel_confidence = .99 # geometry / structure - # TODO: determine lower bounds for grid geometry + # TODO: determine lower bounds for grid geometry that can be resolved grid_limits_x = (-2, 2) grid_limits_y = (-2, 2) grid_limits_z = (0, 4) @@ -123,8 +123,7 @@ def get_beampowers(csdm, gf_spectra, gp_dists): vel = .1 # decimals to round distances to - # TODO: Dynamically determine maximum reasonable accuracy - # should probably be dependent on spatial nyquist? + # TODO: Dynamically determine maximum reasonable accuracy, similar to above decimal_round = 2 # data @@ -134,6 +133,7 @@ def get_beampowers(csdm, gf_spectra, gp_dists): # synth example n_stations = 30 + # place source somewhere in study region and below depth = 2 synth_source_x = np.random.uniform(low=grid_limits_x[0], high=grid_limits_x[1]) synth_source_y = np.random.uniform(low=grid_limits_y[0], high=grid_limits_y[1]) synth_source_z = np.random.uniform(low=2, high=grid_limits_z[1]) @@ -148,11 +148,11 @@ def get_beampowers(csdm, gf_spectra, gp_dists): # compute fft frequencies from scipy.fftpack import fftfreq freqs = fftfreq(n_samples, sampling_rate) - freqs_of_interest = (freqs >= fmin) & (freqs <= fmax) + freqs_of_interest_idx = (freqs >= fmin) & (freqs <= fmax) # compute synthetic spectra dists_to_src = get_distances(list_of_locs=station_locations, point=synth_source) - source_spectra = np.array([get_gf_spectrum(freqs, dist, vel)[freqs_of_interest] for dist in dists_to_src]) + source_spectra = np.array([get_gf_spectrum(freqs, dist, vel)[freqs_of_interest_idx] for dist in dists_to_src]) grid_points, grid_x_coords, grid_y_coords, grid_z_coords = generate_grid( grid_limits=(grid_limits_x, grid_limits_y, grid_limits_z), @@ -173,7 +173,7 @@ def get_beampowers(csdm, gf_spectra, gp_dists): # extract relevant distances only relevant_dists = np.unique(gp_dists) - # compute Green's Functions (spectra) + # compute Green's Functions spectra gf_spectra = get_gf_spectra_for_dists( freqs=freqs, vel=vel, @@ -181,13 +181,14 @@ def get_beampowers(csdm, gf_spectra, gp_dists): ) # limit Green's Functions spectra to frequencies of interest - gf_spectra = gf_spectra[:, freqs_of_interest] + gf_spectra = gf_spectra[:, freqs_of_interest_idx] # compute beampowers beampowers = get_beampowers( csdm=csdm, gf_spectra=gf_spectra, - gp_dists=gp_dists + gp_dists=gp_dists, + relevant_dists=relevant_dists ) # np.argmax(beampowers)