diff --git a/bbp/comps/arias_duration.py b/bbp/comps/arias_duration.py index 77f566be..9e89c654 100644 --- a/bbp/comps/arias_duration.py +++ b/bbp/comps/arias_duration.py @@ -38,9 +38,17 @@ from __future__ import division, print_function # G2CMSS = 980.665 # Convert g to cm/s/s +import os +import sys import math import scipy.integrate +# Import Broadband modules +import bband_utils +import install_cfg +import bbp_formatter +from station_list import StationList + # Compatible with SciPy 1.4 try: cumulative_trapezoid = scipy.integrate.cumulative_trapezoid @@ -169,19 +177,31 @@ def ad_from_acc(a_in_peer_file, a_out_ad): else: ia_norm = [0.0 for i_acc in arias_intensity] - # Define the time for AI=5%, 75%, 95% + # Define the time for AI=5%, 20%, 75%, 80%, 95% time_ai5 = 0 for i in range(pts): if ia_norm[i] >= 5: break time_ai5 = dt * i + time_ai20 = 0 + for i in range(pts): + if ia_norm[i] >= 20: + break + time_ai20 = dt * i + time_ai75 = 0 for i in range(pts): if ia_norm[i] >= 75: break time_ai75 = dt * i + time_ai80 = 0 + for i in range(pts): + if ia_norm[i] >= 80: + break + time_ai80 = dt * i + time_ai95 = 0 for i in range(pts): if ia_norm[i] >= 95: @@ -191,6 +211,7 @@ def ad_from_acc(a_in_peer_file, a_out_ad): # Now, calculate the arias intervals 5% to 75% and 5% to 95% time_5_75 = time_ai75 - time_ai5 time_5_95 = time_ai95 - time_ai5 + time_20_80 = time_ai80 - time_ai20 #print("Arias Intervals: 5 to 75 % 6f (secs), 5 to 95 % 6f (secs) " % # (time_5_75, time_5_95)) @@ -220,8 +241,8 @@ def ad_from_acc(a_in_peer_file, a_out_ad): outfile.write("# Arias Intensities from input accel: %s\n" % a_in_peer_file) outfile.write("# Peak Arias Intensity (cm/sec): %f (secs): %f\n" % (arias_intensity_max, (arias_intensity_index * dt))) - outfile.write("# Arias Intervals: T5-75 %f (s), T5-95 %f (secs)\n" % - (time_5_75, time_5_95)) + outfile.write("# Arias Intervals: T5-75 %f (s), T5-95 %f (s), T20-80 %f (s)\n" % + (time_5_75, time_5_95, time_20_80)) outfile.write("# Seconds Accel (g) " "Arias Intensity (cm/s) " "ADNormalized (%)\n") @@ -230,4 +251,93 @@ def ad_from_acc(a_in_peer_file, a_out_ad): outfile.write("% 8f % 8f % 8f % 8f\n" % (dt_vals[i], acc[i], arias_intensity[i], ia_norm[i])) outfile.close() - return 0 + return arias_intensity_max, time_5_75, time_5_95, time_20_80; + +class AriasDuration(object): + """ + BBP module implementation of arias duration + """ + + def __init__(self, i_r_stations, sim_id=0): + """ + Initializes class variables + """ + self.sim_id = sim_id + self.r_stations = i_r_stations + + def run(self): + print("AriasDuration".center(80, '-')) + # + # convert input bbp acc files to peer format acc files + # + + install = install_cfg.InstallCfg.getInstance() + sim_id = self.sim_id + sta_base = os.path.basename(os.path.splitext(self.r_stations)[0]) + self.log = os.path.join(install.A_OUT_LOG_DIR, str(sim_id), + "%d.araisduration_%s.log" % (sim_id, sta_base)) + a_statfile = os.path.join(install.A_IN_DATA_DIR, + str(sim_id), + self.r_stations) + a_tmpdir = os.path.join(install.A_TMP_DATA_DIR, str(sim_id)) + a_outdir = os.path.join(install.A_OUT_DATA_DIR, str(sim_id)) + + # + # Make sure the tmp and out directories exist + # + bband_utils.mkdirs([a_tmpdir, a_outdir], print_cmd=False) + + slo = StationList(a_statfile) + site_list = slo.get_station_list() + + for site in site_list: + stat = site.scode + print("==> Processing station: %s" % (stat)) + bbpfile = os.path.join(a_outdir, + "%d.%s.acc.bbp" % (sim_id, stat)) + + # Now we need to convert to peer format + out_n_acc = os.path.join(a_tmpdir, + "%d.%s.peer_n.acc" % (sim_id, stat)) + out_e_acc = os.path.join(a_tmpdir, + "%d.%s.peer_e.acc" % (sim_id, stat)) + out_z_acc = os.path.join(a_tmpdir, + "%d.%s.peer_z.acc" % (sim_id, stat)) + bbp_formatter.bbp2peer(bbpfile, out_n_acc, out_e_acc, out_z_acc, accel=True) + + # Duration output files + out_n_arias = os.path.join(a_outdir, + "%d.%s_N.arias" % (sim_id, stat)) + out_e_arias = os.path.join(a_outdir, + "%d.%s_E.arias" % (sim_id, stat)) + out_z_arias = os.path.join(a_outdir, + "%d.%s_Z.arias" % (sim_id, stat)) + + # compute each one + n_max, n_time_5_75, n_time_5_95, n_time_20_80 = ad_from_acc(out_n_acc, out_n_arias) + e_max, e_time_5_75, e_time_5_95, e_time_20_80 = ad_from_acc(out_e_acc, out_e_arias) + z_max, z_time_5_75, z_time_5_95, z_time_20_80 = ad_from_acc(out_z_acc, out_z_arias) + geo_max = math.sqrt(n_max*e_max) + geo_time_5_75 = math.sqrt(n_time_5_75*e_time_5_75) + geo_time_5_95 = math.sqrt(n_time_5_95*e_time_5_95) + geo_time_20_80 = math.sqrt(n_time_20_80*e_time_20_80) + + # write summary file + out_duration = os.path.join(a_outdir, + "%d.%s.ard" % (sim_id, stat)) + outfile = open(out_duration, "w") + outfile.write("# Component Peak Arias (cm/s) T5-75 (s) T5-95 (s) T20-80 (s)\n") + outfile.write("N %f %f %f %f\n" % (n_max, n_time_5_75, n_time_5_95, n_time_20_80)) + outfile.write("E %f %f %f %f\n" % (e_max, e_time_5_75, e_time_5_95, e_time_20_80)) + outfile.write("Z %f %f %f %f\n" % (z_max, z_time_5_75, z_time_5_95, z_time_20_80)) + outfile.write("GEOM %f %f %f %f\n" % (geo_max, geo_time_5_75, geo_time_5_95, geo_time_20_80)) + outfile.close() + + # All done! + print("AriasDuration Completed".center(80, '-')) + +if __name__ == '__main__': + print("Testing Module: %s" % (os.path.basename(sys.argv[0]))) + ME = AriasDuration(sys.argv[1], sim_id=int(sys.argv[2])) + ME.run() + sys.exit(0) diff --git a/bbp/comps/bbp_formatter.py b/bbp/comps/bbp_formatter.py index d3d343d9..68bd2aa0 100644 --- a/bbp/comps/bbp_formatter.py +++ b/bbp/comps/bbp_formatter.py @@ -189,7 +189,7 @@ def peer2bbp(in_peer_n_file, in_peer_e_file, in_peer_z_file, out_bbp_file): # Lastly, close the file bbp_file.close() -def bbp2peer(in_bbp_file, out_peer_n_file, out_peer_e_file, out_peer_z_file): +def bbp2peer(in_bbp_file, out_peer_n_file, out_peer_e_file, out_peer_z_file, accel=True): """ Convert bbp file into three peer files for use by RotD50 and other programs that input PEER format seismograms @@ -249,9 +249,14 @@ def bbp2peer(in_bbp_file, out_peer_n_file, out_peer_e_file, out_peer_z_file): "Error in conversion.") else: dt_vals.append(dt) - n_vals.append(float(elems[1]) / bband_utils.G2CMSS) - e_vals.append(float(elems[2]) / bband_utils.G2CMSS) - z_vals.append(float(elems[3]) / bband_utils.G2CMSS) + if accel: + n_vals.append(float(elems[1]) / bband_utils.G2CMSS) + e_vals.append(float(elems[2]) / bband_utils.G2CMSS) + z_vals.append(float(elems[3]) / bband_utils.G2CMSS) + else: + n_vals.append(float(elems[1])) + e_vals.append(float(elems[2])) + z_vals.append(float(elems[3])) # Prepare to write 6 colume format n_file = open(out_peer_n_file, "w") @@ -271,11 +276,16 @@ def bbp2peer(in_bbp_file, out_peer_n_file, out_peer_e_file, out_peer_z_file): e_file.write(line) z_file.write(line) - n_file.write("Acceleration in g\n") + if accel: + n_file.write("Acceleration in g\n") + e_file.write("Acceleration in g\n") + z_file.write("Acceleration in g\n") + else: + n_file.write("Velicity in cm/sec\n") + e_file.write("Velicity in cm/sec\n") + z_file.write("Velicity in cm/sec\n") n_file.write(" %d %1.6f NPTS, DT\n" % (npts, dt)) - e_file.write("Acceleration in g\n") e_file.write(" %d %1.6f NPTS, DT\n" % (npts, dt)) - z_file.write("Acceleration in g\n") z_file.write(" %d %1.6f NPTS, DT\n" % (npts, dt)) cur_line = 0 diff --git a/bbp/comps/module.py b/bbp/comps/module.py index d710878c..ba747853 100644 --- a/bbp/comps/module.py +++ b/bbp/comps/module.py @@ -49,6 +49,7 @@ from uc_site import UCSite from wcc_siteamp import WccSiteamp from rotd50 import RotD50 +from rotd_vel import RotDVel from fas import FAS from obs_seismograms import ObsSeismograms from copy_seismograms import CopySeismograms @@ -75,6 +76,7 @@ from irikura_gen_srf import IrikuraGenSrf from irikura_hf import IrikuraHF from seismo_soil import SeismoSoil +from arias_duration import AriasDuration from lf_seismograms import LFSeismograms class Module(object): diff --git a/bbp/comps/rotd_vel.py b/bbp/comps/rotd_vel.py new file mode 100644 index 00000000..e6d1675e --- /dev/null +++ b/bbp/comps/rotd_vel.py @@ -0,0 +1,181 @@ +#!/usr/bin/env python3 +""" +BSD 3-Clause License + +Copyright (c) 2023, University of Southern California +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +""" +from __future__ import division, print_function + +# Import Python modules +import os +import sys + +# Import Broadband modules +import bband_utils +import install_cfg +import bbp_formatter +from station_list import StationList + +class RotDVel(object): + """ + BBP module implementation of rotd50 for velocity provided by UCB. + Rotd50 inputs seismograms and outputs response spectra + Largely copied from RotD50 with minor modifications by Kevin Milner for velocity. + """ + @staticmethod + def do_rotd_vel(workdir, peer_input_e_file, peer_input_n_file, + peer_input_z_file, output_rotd_file, + logfile): + """ + This function runs the rotd50 command inside workdir, using + the inputs and outputs specified + """ + install = install_cfg.InstallCfg.getInstance() + + # Make sure we don't have absolute path names + peer_input_e_file = os.path.basename(peer_input_e_file) + peer_input_n_file = os.path.basename(peer_input_n_file) + peer_input_z_file = os.path.basename(peer_input_z_file) + output_rotd_file = os.path.basename(output_rotd_file) + + # Save cwd, change back to it at the end + old_cwd = os.getcwd() + os.chdir(workdir) + + # Make sure we remove the output files first or Fortran will + # complain if they already exist + try: + os.unlink(output_rotd_file) + except OSError: + pass + + # + # write config file for rotd100 program + rd50_conf = open("rotd100_inp.cfg", 'w') + # This flag indicates inputs acceleration + rd50_conf.write("2 interp flag\n") + # This flag indicate we are processing two input files (horizontals) + rd50_conf.write("1 Npairs\n") + # Number of headers in the file + rd50_conf.write("6 Nhead\n") + rd50_conf.write("%s\n" % peer_input_e_file) + rd50_conf.write("%s\n" % peer_input_n_file) + rd50_conf.write("%s\n" % output_rotd_file) + # Close file + rd50_conf.close() + + progstring = ("%s/rotd100 >> %s 2>&1" % (install.A_UCB_BIN_DIR, logfile)) + bband_utils.runprog(progstring, abort_on_error=True, print_cmd=False) + + # Restore working directory + os.chdir(old_cwd) + + def __init__(self, i_r_stations, sim_id=0): + """ + Initializes class variables + """ + self.sim_id = sim_id + self.r_stations = i_r_stations + + def run(self): + print("RotDVel".center(80, '-')) + # + # convert input bbp acc files to peer format acc files + # + + install = install_cfg.InstallCfg.getInstance() + sim_id = self.sim_id + sta_base = os.path.basename(os.path.splitext(self.r_stations)[0]) + self.log = os.path.join(install.A_OUT_LOG_DIR, str(sim_id), + "%d.rotdvel_%s.log" % (sim_id, sta_base)) + a_statfile = os.path.join(install.A_IN_DATA_DIR, + str(sim_id), + self.r_stations) + a_tmpdir = os.path.join(install.A_TMP_DATA_DIR, str(sim_id)) + a_outdir = os.path.join(install.A_OUT_DATA_DIR, str(sim_id)) + + # + # Make sure the tmp and out directories exist + # + bband_utils.mkdirs([a_tmpdir, a_outdir], print_cmd=False) + + slo = StationList(a_statfile) + site_list = slo.get_station_list() + + for site in site_list: + stat = site.scode + print("==> Processing station: %s" % (stat)) + + # Create path names and check if their sizes are within bounds + nsfile = os.path.join(a_tmpdir, + "%d.%s.000" % (sim_id, stat)) + ewfile = os.path.join(a_tmpdir, + "%d.%s.090" % (sim_id, stat)) + udfile = os.path.join(a_tmpdir, + "%d.%s.ver" % (sim_id, stat)) + bbpfile = os.path.join(a_outdir, + "%d.%s.vel.bbp" % (sim_id, stat)) + + bband_utils.check_path_lengths([nsfile, ewfile, udfile], + bband_utils.GP_MAX_FILENAME) + + cmd = ("%s/wcc2bbp " % (install.A_GP_BIN_DIR) + + "wcc2bbp=0 nsfile=%s ewfile=%s udfile=%s < %s >> %s 2>&1" % + (nsfile, ewfile, udfile, bbpfile, self.log)) + bband_utils.runprog(cmd, abort_on_error=True, print_cmd=False) + + # Now we need to convert to peer format + out_n_vel = os.path.join(a_tmpdir, + "%d.%s.peer_n.vel" % (sim_id, stat)) + out_e_vel = os.path.join(a_tmpdir, + "%d.%s.peer_e.vel" % (sim_id, stat)) + out_z_vel = os.path.join(a_tmpdir, + "%d.%s.peer_z.vel" % (sim_id, stat)) + bbp_formatter.bbp2peer(bbpfile, out_n_vel, out_e_vel, out_z_vel, accel=False) + + # Let's have rotD50 create these output files + out_rotd_base = "%d.%s.rdvel" % (sim_id, stat) + tmp_rotd = os.path.join(a_tmpdir, out_rotd_base) + out_rotd = os.path.join(a_outdir, out_rotd_base) + + # Run the rotD50 program + self.do_rotd_vel(a_tmpdir, out_e_vel, out_n_vel, out_z_vel, + out_rotd, self.log) + + cmd = "cp %s %s" % (tmp_rotd, out_rotd) + bband_utils.runprog(cmd, abort_on_error=True, print_cmd=False) + + # All done! + print("RotDVel Completed".center(80, '-')) + +if __name__ == '__main__': + print("Testing Module: %s" % (os.path.basename(sys.argv[0]))) + ME = RotDVel(sys.argv[1], sim_id=int(sys.argv[2])) + ME.run() + sys.exit(0) diff --git a/bbp/src/ucb/rotd50/rotd100.f b/bbp/src/ucb/rotd50/rotd100.f index 25f60486..e712afca 100644 --- a/bbp/src/ucb/rotd50/rotd100.f +++ b/bbp/src/ucb/rotd50/rotd100.f @@ -30,22 +30,22 @@ program Calc_RotD50 real dt10 real x0(MAXPTS), y0(MAXPTS), u(MAXPTS), y2(MAXPTS) real rspTH1(MAXPTS), rspTH2(MAXPTS), rsp1(MAXPTS), rsp2(MAXPTS) - real x(MAXPTS), y(MAXPTS), rsp_Period(63), famp15(3) + real x(MAXPTS), y(MAXPTS), rsp_Period(65), famp15(3) real rotangle, w(200), sa(1000), workArray(1000), rotD50(3,200), rotD100(3,200), psa5E(200), psa5N(200) integer rD100ang(3,200), rD50ang(3,200) real saUnsort(1000) real damping complex cu1(MAXPTS) - data RSP_Period / 0.010, 0.011, 0.012, 0.013, 0.015, 0.017, 0.020, 0.022, 0.025, 0.029, - 1 0.032, 0.035, 0.040, 0.045, 0.050, 0.055, 0.060, 0.065, 0.075, 0.085, - 2 0.100, 0.110, 0.120, 0.130, 0.150, 0.170, 0.200, 0.220, 0.240, 0.260, - 3 0.280, 0.300, 0.350, 0.400, 0.450, 0.500, 0.550, 0.600, 0.650, 0.750, - 4 0.850, 1.000, 1.100, 1.200, 1.300, 1.500, 1.700, 2.000, 2.200, 2.400, - 5 2.600, 2.800, 3.000, 3.500, 4.000, 4.400, 5.000, 5.500, 6.000, 6.500, - 6 7.500, 8.500, 10.000 / + data RSP_Period / 0.0001, 0.001, 0.010, 0.011, 0.012, 0.013, 0.015, 0.017, 0.020, 0.022, + 1 0.025, 0.029, 0.032, 0.035, 0.040, 0.045, 0.050, 0.055, 0.060, 0.065, + 2 0.075, 0.085, 0.100, 0.110, 0.120, 0.130, 0.150, 0.170, 0.200, 0.220, + 3 0.240, 0.260, 0.280, 0.300, 0.350, 0.400, 0.450, 0.500, 0.550, 0.600, + 4 0.650, 0.750, 0.850, 1.000, 1.100, 1.200, 1.300, 1.500, 1.700, 2.000, + 5 2.200, 2.400, 2.600, 2.800, 3.000, 3.500, 4.000, 4.400, 5.000, 5.500, + 6 6.000, 6.500, 7.500, 8.500, 10.000 / - nFreq = 63 + nFreq = 65 damping = 0.05 dt_max = 0.001