diff --git a/bbp/comps/anderson_gof.py b/bbp/comps/anderson_gof.py index ffb03272..e35f085d 100755 --- a/bbp/comps/anderson_gof.py +++ b/bbp/comps/anderson_gof.py @@ -2,7 +2,7 @@ """ BSD 3-Clause License -Copyright (c) 2023, University of Southern California +Copyright (c) 2024, University of Southern California All rights reserved. Redistribution and use in source and binary forms, with or without @@ -51,6 +51,16 @@ from scipy import signal, stats from scipy.signal import butter, filtfilt +# Compatible with SciPy 1.4 +try: + cumulative_trapezoid = integrate.cumulative_trapezoid +except AttributeError: + cumulative_trapezoid = integrate.cumtrapz +try: + simpson = integrate.simpson +except AttributeError: + simpson = integrate.simps + # Import BBP modules import bband_utils import plot_config @@ -135,15 +145,15 @@ def eval_func2(p1, p2): @staticmethod def integral(x, idt): - return integrate.simps(x, x=None, dx=idt) + return simpson(x, x=None, dx=idt) @staticmethod def integral2(x, idt): - return integrate.simps(x * x, x=None, dx=idt) + return simpson(x * x, x=None, dx=idt) @staticmethod def integ(array_in, idt): - return integrate.cumtrapz(array_in, dx=idt) + return cumulative_trapezoid(array_in, dx=idt) @staticmethod def shift_bit_length(x): diff --git a/bbp/comps/arias_duration.py b/bbp/comps/arias_duration.py index fad42b23..77f566be 100644 --- a/bbp/comps/arias_duration.py +++ b/bbp/comps/arias_duration.py @@ -1,7 +1,8 @@ +#!/usr/bin/env python3 """ BSD 3-Clause License -Copyright (c) 2021, University of Southern California +Copyright (c) 2024, University of Southern California All rights reserved. Redistribution and use in source and binary forms, with or without @@ -40,6 +41,12 @@ import math import scipy.integrate +# Compatible with SciPy 1.4 +try: + cumulative_trapezoid = scipy.integrate.cumulative_trapezoid +except AttributeError: + cumulative_trapezoid = scipy.integrate.cumtrapz + # Converting to cm units. Use approximation to g G_TO_CMS = 981.0 # %(cm/s) @@ -97,9 +104,9 @@ def ad_from_acc(a_in_peer_file, a_out_ad): # Integrate and calculate vel and disp arrays # acc_cms = [value * G_TO_CMS for value in acc] - velo = list(dt * scipy.integrate.cumtrapz(acc_cms)) + velo = list(dt * cumulative_trapezoid(acc_cms)) velo.insert(0, 0) - disp = list(dt * scipy.integrate.cumtrapz(velo)) + disp = list(dt * cumulative_trapezoid(velo)) disp.insert(0, 0) # @@ -138,7 +145,7 @@ def ad_from_acc(a_in_peer_file, a_out_ad): # Arias Intensities # Using the trapezoidal integration arias_intensity = [pow((value * G_TO_CMS), 2) for value in acc] - tsum = list(scipy.integrate.cumtrapz(arias_intensity) * dt) + tsum = list(cumulative_trapezoid(arias_intensity) * dt) tsum.insert(0, 0) arias_intensity = [num * math.pi / (2 * G_TO_CMS) for num in tsum] diff --git a/bbp/comps/jbsim.py b/bbp/comps/jbsim.py index 6733b422..74878e5f 100755 --- a/bbp/comps/jbsim.py +++ b/bbp/comps/jbsim.py @@ -2,7 +2,7 @@ """ BSD 3-Clause License -Copyright (c) 2023, University of Southern California +Copyright (c) 2024, University of Southern California All rights reserved. Redistribution and use in source and binary forms, with or without @@ -78,9 +78,18 @@ def run(self): a_statfile = os.path.join(install.A_IN_DATA_DIR, str(sim_id), self.r_stations) - a_srffile = os.path.join(install.A_IN_DATA_DIR, - str(sim_id), - self.r_srffile) + a_rupture_file = os.path.join(install.A_IN_DATA_DIR, + str(sim_id), + self.r_srffile) + # RWG 20241025: jbsim-v3.0.0 is backward compatible with jbsim-v2.0.0 + # it can use rupture files in both SRF and MRF formats + if os.path.splitext(a_rupture_file)[-1].lower() == ".srf": + rupmod_type = "SRF" + elif os.path.splitext(a_rupture_file)[-1].lower() == ".mrf": + rupmod_type = "MRF" + else: + print("[ERROR]: Unknown rupture format in file %s" % (a_rupture_file)) + sys.exit(-1) # Set directories, and make sure they exist a_indir = os.path.join(install.A_IN_DATA_DIR, str(sim_id)) @@ -106,9 +115,11 @@ def run(self): # progstring = ("%s latloncoords=1 slon=%f slat=%f " % (os.path.join(install.A_GP_BIN_DIR, - "jbsim-v2.0.0"), slon, slat) + + "jbsim-v3.0.0"), slon, slat) + "tshift_timedomain=1 use_closest_gf=1 " + - "rupmodtype=SRF rupmodfile=%s " % a_srffile + + "cgf_flag=%d " % (config.GF_CGFFLAG) + + "rupmodtype=%s " % (rupmod_type) + + "rupmodfile=%s " % (a_rupture_file) + "moment=-1 outdir=%s stat=%s " % (a_veldir, site) + "min_taper_range=0.0 max_taper_range=0.0 " + "gftype=fk gflocs=%s gftimes=%s gf_swap_bytes=%d " % diff --git a/bbp/comps/jbsim_cfg.py b/bbp/comps/jbsim_cfg.py index f3f116c8..736a6894 100755 --- a/bbp/comps/jbsim_cfg.py +++ b/bbp/comps/jbsim_cfg.py @@ -1,8 +1,8 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 """ BSD 3-Clause License -Copyright (c) 2021, University of Southern California +Copyright (c) 2024, University of Southern California All rights reserved. Redistribution and use in source and binary forms, with or without @@ -93,6 +93,11 @@ def __init__(self, vmodel_name): self.NTOUT = int(vmodel_params['NTOUT']) else: self.NTOUT = 4096 + # RWG 20241025: Added parameter GF_CGFFLAG + if 'GF_CGFFLAG' in vmodel_params: + self.GF_CGFFLAG = int(vmodel_params['GF_CGFFLAG']) + else: + self.GF_CGFFLAG = 0 # Now, configure other paramteres self.GF_SWAP_BYTES = 0 diff --git a/bbp/comps/match.py b/bbp/comps/match.py index a22d448f..a4682695 100755 --- a/bbp/comps/match.py +++ b/bbp/comps/match.py @@ -2,7 +2,7 @@ """ BSD 3-Clause License -Copyright (c) 2023, University of Southern California +Copyright (c) 2024, University of Southern California All rights reserved. Redistribution and use in source and binary forms, with or without @@ -87,11 +87,25 @@ def run(self): dirs = [a_tmpdir] bband_utils.mkdirs(dirs, print_cmd=False) + # Get pointer to the velocity model object + vel_obj = velocity_models.get_velocity_model_by_name(self.vmodel_name) + if vel_obj is None: + raise bband_utils.ParameterError("Cannot find velocity model: %s" % + (self.vmodel_name)) + + # Check for velocity model-specific parameters + vmodel_params = vel_obj.get_codebase_params('gp') + # Start with defaults self.phase = config.PHASE self.hf_fhi = config.HF_FHI self.lf_flo = config.LF_FLO + # Check if we have a different merging frequency + if 'MATCH_MERGING_FREQUENCY' in vmodel_params: + self.hf_fhi = float(vmodel_params['MATCH_MERGING_FREQUENCY']) + self.lf_flo = float(vmodel_params['MATCH_MERGING_FREQUENCY']) + # Set match method if config.MATCH_METHOD == 1: self.phase = 1 @@ -109,15 +123,6 @@ def run(self): slo = StationList(a_statfile) site_list = slo.get_station_list() - # Get pointer to the velocity model object - vel_obj = velocity_models.get_velocity_model_by_name(self.vmodel_name) - if vel_obj is None: - raise bband_utils.ParameterError("Cannot find velocity model: %s" % - (self.vmodel_name)) - - # Check for velocity model-specific parameters - vmodel_params = vel_obj.get_codebase_params('gp') - # Figure out what DT we should use when resampling # Figure out the LF DT value diff --git a/bbp/comps/match_cfg.py b/bbp/comps/match_cfg.py index f1ee9453..ee32db9c 100755 --- a/bbp/comps/match_cfg.py +++ b/bbp/comps/match_cfg.py @@ -1,8 +1,8 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 """ BSD 3-Clause License -Copyright (c) 2021, University of Southern California +Copyright (c) 2024, University of Southern California All rights reserved. Redistribution and use in source and binary forms, with or without @@ -45,13 +45,16 @@ class MatchCfg(object): Define the configuration parameters for the Match program """ def __init__(self): - self.HF_FHI = 1.0 + # Default merging frequency set + self.DEFAULT_MERGING_FREQUENCY = 1.0 + + self.HF_FHI = self.DEFAULT_MERGING_FREQUENCY self.HF_FLO = "1.0e+15" self.HF_ORD = 4 self.HF_TSTART = 0.0 self.LF_FHI = 0.0 - self.LF_FLO = 1.0 + self.LF_FLO = self.DEFAULT_MERGING_FREQUENCY self.LF_ORD = 4 self.LF_TSTART = 0.0 diff --git a/bbp/comps/rzz2015.py b/bbp/comps/rzz2015.py index 00187bad..5e2c935c 100755 --- a/bbp/comps/rzz2015.py +++ b/bbp/comps/rzz2015.py @@ -2,7 +2,7 @@ """ BSD 3-Clause License -Copyright (c) 2023, University of Southern California +Copyright (c) 2024, University of Southern California All rights reserved. Redistribution and use in source and binary forms, with or without @@ -45,6 +45,12 @@ mpl.use('Agg') # Disable use of Tk/X11 import pylab +# Compatible with SciPy 1.4 +try: + cumulative_trapezoid = integrate.cumulative_trapezoid +except AttributeError: + cumulative_trapezoid = integrate.cumtrapz + # Import BBP modules import bband_utils import plot_config @@ -106,7 +112,7 @@ def arias(self, F, dt, percent): n = len(F) a_i = [pow(value, 2) for value in F] - I = integrate.cumtrapz(a_i) * dt + I = cumulative_trapezoid(a_i) * dt # Arias Intensity Ia = (F[0]**2) * dt / 2.0 + I[n-2] + (F[n-1]**2) * dt / 2.0 It = (percent / 100.0) * Ia diff --git a/bbp/src/gp/JordanBailey/.gitignore b/bbp/src/gp/JordanBailey/.gitignore index 52e5d527..6b6d3bde 100644 --- a/bbp/src/gp/JordanBailey/.gitignore +++ b/bbp/src/gp/JordanBailey/.gitignore @@ -2,4 +2,5 @@ jbsim jbsim3d ray_stimes -jbsim-v2.0.0 \ No newline at end of file +jbsim-v2.0.0 +jbsim-v3.0.0 \ No newline at end of file diff --git a/bbp/src/gp/JordanBailey/defs.h b/bbp/src/gp/JordanBailey/defs.h index f7e8e630..9fb4af7c 100644 --- a/bbp/src/gp/JordanBailey/defs.h +++ b/bbp/src/gp/JordanBailey/defs.h @@ -11,6 +11,7 @@ #define U1FLAG 1 #define U2FLAG 2 #define U3FLAG 3 +#define MOMTEN_FLAG 4 #define RDONLY_FLAGS O_RDONLY #define RDWR_FLAGS O_RDWR diff --git a/bbp/src/gp/JordanBailey/function.h b/bbp/src/gp/JordanBailey/function.h index a70d8840..6a4ccb64 100644 --- a/bbp/src/gp/JordanBailey/function.h +++ b/bbp/src/gp/JordanBailey/function.h @@ -36,6 +36,8 @@ void sum_4gf_v2(float *,int,float *,struct gfheader *,int,int,float *,float *,fl void timeshift_gf_v2(float *,int,float *,struct gfheader *,int,float *,float *,int); void mech_4gf_v2(float *,float *,struct gfheader *,struct gfparam gfp,int,struct mechparam,float *,float *); +void momten_4gf(struct standrupformat *,int,int,float *,float *,struct gfheader *,struct gfparam,int,float *); + void read_beroza(struct beroza *,char *rfile,float *); void get_brmpars(struct beroza *,int,int,float *,float *,float *,float *); void beroza_stf(struct beroza *,int,int,float *,float *,float *,int,float *,float *); @@ -57,6 +59,9 @@ void get_srfpars(struct standrupformat *,int,int,float *,float *,float *,float * void get_srfpars_v2(struct standrupformat *,int,int,float *,float *,struct mechparam *); void get_srfparsOLD(struct standrupformat *,int,int,float *,float *,float *,int,int,struct mechparam *); void srf_stf(struct standrupformat *,int,int,float *,float *,float *,int,float *,struct mechparam,float *); +void srf_stf_v2(struct standrupformat *,int,int,float *,float *,float *,int,float *,struct mechparam,float *); + +void get_indx_pow(float *,float *,float *,struct sgtindex *,struct geoprojection *,int); char *skipval(int,char *); void do_cnvlv(float *,float *,int,float *,int); diff --git a/bbp/src/gp/JordanBailey/greenfunc.c b/bbp/src/gp/JordanBailey/greenfunc.c index 75ffb403..d58cee2b 100644 --- a/bbp/src/gp/JordanBailey/greenfunc.c +++ b/bbp/src/gp/JordanBailey/greenfunc.c @@ -318,24 +318,32 @@ if(gfpar.flag3d == 0) fdr = opfile_ro(str); + if(gfpar.cgf_flag == 1) + reed(fdr,&gfh->nc,sizeof(int)); reed(fdr,&gfh->nt,sizeof(int)); reed(fdr,&gfh->dt,sizeof(float)); reed(fdr,&gfh->dep,sizeof(float)); reed(fdr,&gfh->rng,sizeof(float)); reed(fdr,&gfh->tst,sizeof(float)); reed(fdr,&gfh->mom,sizeof(float)); + if(gfpar.cgf_flag == 1) + reed(fdr,&gfh->lam,sizeof(float)); reed(fdr,&gfh->mu,sizeof(float)); if(gfpar.swap_flag == 0) /* byte-order is OK, no swapping needed */ swap_flag = 0; else /* byte-order not OK, swapping needed */ { + if(gfpar.cgf_flag == 1) + swap_in_place(1,(char *)(&gfh->nc)); swap_in_place(1,(char *)(&gfh->nt)); swap_in_place(1,(char *)(&gfh->dt)); swap_in_place(1,(char *)(&gfh->dep)); swap_in_place(1,(char *)(&gfh->rng)); swap_in_place(1,(char *)(&gfh->tst)); swap_in_place(1,(char *)(&gfh->mom)); + if(gfpar.cgf_flag == 1) + swap_in_place(1,(char *)(&gfh->lam)); swap_in_place(1,(char *)(&gfh->mu)); swap_flag = 1; @@ -1146,7 +1154,7 @@ for(kg=0;kg<4;kg++) /* resample if needed */ if((*dtout)/gfh[kg].dt > pratio_tol || (*dtout)/gfh[kg].dt < mratio_tol) { - fprintf(stderr,"*** RESAMPLED diff= %13.5e ratio= %13.5e\n",(*dtout)-gfh[kg].dt,(*dtout)/gfh[kg].dt); + fprintf(stderr,"*** RESAMPLED dtout= %.5e gfdt= %.5e diff= %13.5e ratio= %13.5e\n",*dtout,gfh[kg].dt,(*dtout)-gfh[kg].dt,(*dtout)/gfh[kg].dt); ntpad = 2*gfh[kg].nt; fnt = ntpad*gfh[kg].dt/(*dtout); @@ -1497,6 +1505,10 @@ else while(ntp2 < ntsum) ntp2 = ntp2*2; + maxnt = ntout; + if(ntp2 < maxnt) + maxnt = ntp2; + nf2 = ntp2/2; for(ig=0;ig<4;ig++) @@ -1561,53 +1573,22 @@ else gc2[i].im = gc2[i].re*sinA + gc2[i].im*cosA; gc2[i].re = tmpre; } - } - - /* reset pointers to first group of 8 GFs */ - - gf0 = gf; - gf1 = gf + ntsum; - gf2 = gf + 2*ntsum; - - gc0 = (struct complex *) gf0; - gc1 = (struct complex *) gf1; - gc2 = (struct complex *) gf2; - - /* sum over 4 GFs to interpolated response */ - - nts3 = 3*nf2; - nts6 = 6*nf2; - nts9 = 9*nf2; - for(i=0;i= *minr) } } } + +/* RWG 20241018: Added full moment tensor option using CGF and MRF input. + +The cgf's are used to convolve with moment tensor source representations. Conversion of cgf +to full 18 component strain Green's tensor is outlined below. + +Using A&R convention for Green's Tensor: + + n = positive North + e = positive East + d = positive Down + + Gnn,Gee,Gdd,Gne,Gnd,Ged for each component n,e,d + +Transformation from Lupei Zhu (fk2mt in radiat.c): + +Note that in Lupei's z,r,t system z=Up. To explicitly reflect this below, +I have replaced Lupei's "z" with "u" for the ground motion component. The "z" +for the moment tensor is still "down", consistent with A&R. + +This means all of Lupei's .u components are the opposite polarity of my .d components, +and thus my conversions for the .d comps are multiplied by -1 relative to Lupei's formulas. +Additionally, I believe these formula also include a factor of 2 for the off-diagonal components, +so that symmetry of the tensors can be used to reduce computations. + +In these relations: cosAz=cos(azi), cos2Az=cos(2*azi), sinAz=sin(azi), sin2Az=sin(2*azi). + + Gnn_r = rex/3 - rdd/6 - rss*cos2Az/2 + Gnn_t = -tss*sin2Az/2 + Gnn_d = -(uex/3 - udd/6 - uss*cos2Az/2) + + Gee_r = rex/3 - rdd/6 + rss*cos2Az/2 + Gee_t = -Gnn_t + Gee_d = -(uex/3 - udd/6 + uss*cos2Az/2) + + Gdd_r = rex/3 + rdd/3 + Gdd_t = 0 + Gdd_d = -(uex/3 + udd/3) + + Gne_r = -rss*sin2Az + Gne_t = tss*cos2Az + Gne_d = uss*sin2Az + + Gnd_r = -rds*cosAz + Gnd_t = -tds*sinAz + Gnd_d = uds*cosAz + + Ged_r = -rds*sinAz + Ged_t = tds*cosAz + Ged_d = uds*sinAz + +*/ + +void momten_4gf(struct standrupformat *srf,int off,int ip,float *gfmech,float *gf,struct gfheader *gfh,struct gfparam gfp,int nts,float *azi) +{ +float *udd, *rdd, *uds, *rds, *tds, *uss, *rss, *tss, *uex, *rex; +float *gfn, *gfe, *gfu; +double Mnn, Mee, Mdd, Mne, Mnd, Med; +double Gnn_r, Gee_r, Gdd_r, Gne_r, Gnd_r, Ged_r; +double Gnn_t, Gee_t, Gdd_t, Gne_t, Gnd_t, Ged_t; +double Gnn_d, Gee_d, Gdd_d, Gne_d, Gnd_d, Ged_d; +double scale, cosAz, sinAz, cos2Az, sin2Az; +double rad, tan; +int it, ig; + +struct srf_planerectangle *prect_ptr; +struct srf_prectsegments *prseg_ptr; +struct srf_allpoints *apnts_ptr; +struct srf_apointvalues *apval_ptr; + +float half, third, sixth; + +half = 1.0/2.0; +third = 1.0/3.0; +sixth = 1.0/6.0; + +prect_ptr = &(srf->srf_prect); +prseg_ptr = prect_ptr->prectseg; +apnts_ptr = &(srf->srf_apnts); +apval_ptr = apnts_ptr->apntvals + off; + +Mnn = apval_ptr[ip].mnn; +Mee = apval_ptr[ip].mee; +Mdd = apval_ptr[ip].mdd; +Mne = apval_ptr[ip].mne; +Mnd = apval_ptr[ip].mnd; +Med = apval_ptr[ip].med; + +zapit(gfmech,12*nts); + +cosAz = cos(*azi); +sinAz = sin(*azi); + +cos2Az = cosAz*cosAz - sinAz*sinAz; +sin2Az = 2.0*sinAz*cosAz; + +for(ig=0;ig<4;ig++) + { + if(gfh[ig].wt > (float)(0.0)) + { + gfu = gfmech + 3*ig*nts; + gfn = gfmech + 3*ig*nts + nts; + gfe = gfmech + 3*ig*nts + 2*nts; + + udd = gf + gfp.nc*ig*nts; + rdd = gf + gfp.nc*ig*nts + nts; + uds = gf + gfp.nc*ig*nts + 2*nts; + rds = gf + gfp.nc*ig*nts + 3*nts; + tds = gf + gfp.nc*ig*nts + 4*nts; + uss = gf + gfp.nc*ig*nts + 5*nts; + rss = gf + gfp.nc*ig*nts + 6*nts; + tss = gf + gfp.nc*ig*nts + 7*nts; + uex = gf + gfp.nc*ig*nts + 8*nts; + rex = gf + gfp.nc*ig*nts + 9*nts; + + scale = (gfh[ig].wt)/(gfh[ig].mom); + + for(it=0;it cm^2 */ + +float half = 0.5; +float two = 2.0; + +int latloncoords = 0; +int geoproj = 1; + +char statfile[2048]; +int nstat, istat; +int nbuf = 1000; +char *stat_buf, *stat_ptr; +float *slon_buf, *slat_buf; + +float tstart = 0.0; + +sprintf(rupmodtype,"SRF"); +statfile[0] = '\0'; + +sprintf(gfpar.gftype,"fk"); + +setpar(ac, av); + +getpar("statfile","s",statfile); +getpar("latloncoords","d",&latloncoords); +getpar("geoproj","d",&geoproj); + +if(statfile[0] != '\0') + latloncoords = 0; +else + latloncoords = 1; + +if(latloncoords == 1) + { + mstpar("slat","f",&slat); + mstpar("slon","f",&slon); + } + +getpar("rupmodtype","s",rupmodtype); +if(strcmp(rupmodtype,"SRF") != 0 && strcmp(rupmodtype,"srf") != 0 + && strcmp(rupmodtype,"MRF") != 0 && strcmp(rupmodtype,"mrf") != 0) + { + sprintf(rupmodtype,"SRF"); + } + +gfpar.cgf_flag = 0; +getpar("cgf_flag","d",&gfpar.cgf_flag); + +gfpar.nc = 8; +if(gfpar.cgf_flag > 0) + gfpar.nc = 10; + +srf_flag = 0; +mrf_flag = 0; +if(strcmp(rupmodtype,"SRF") == 0 || strcmp(rupmodtype,"srf") == 0) + { + srf_flag = 1; + mrf_flag = 0; + } +if(strcmp(rupmodtype,"MRF") == 0 || strcmp(rupmodtype,"mrf") == 0) + { + srf_flag = 0; + mrf_flag = 1; + + if(gfpar.cgf_flag == 0) + { + fprintf(stderr,"cgf_flag= %d: input GFs must be Complete GFs for moment tensor source, exiting...\n",gfpar.cgf_flag); + exit(-1); + } + } + +if(srf_flag == 1 || mrf_flag == 1) + { + mstpar("rupmodfile","s",rupmodfile); + + getpar("nseg","d",&nseg); + getpar("inbin","d",&inbin); + + mstpar("outdir","s",outdir); + if(latloncoords == 1) + mstpar("stat","s",stat); + } + +if((strncmp(gfpar.gftype,"fk",2) == 0) || (strncmp(gfpar.gftype,"FK",2) == 0)) + { + gfpar.flag3d = 0; + + mstpar("gflocs","s",gfpar.gflocs); + mstpar("gftimes","s",gfpar.gftimes); + + gfpar.swap_flag = 0; + getpar("gf_swap_bytes","d",&gfpar.swap_flag); + } +else + { + fprintf(stderr,"gftype= %s invalid option, exiting...\n",gfpar.gftype); + exit(-1); + } + +mstpar("gfpath","s",gfpath); +mstpar("gfname","s",gfname); + +getpar("use_closest_gf","d",&use_closest_gf); +getpar("min_taper_range","f",&min_taper_range); +getpar("max_taper_range","f",&max_taper_range); + +mstpar("maxnt","d",&maxnt); +mstpar("mindt","f",&mindt); +getpar("ntout","d",&ntout); +getpar("dtout","f",&dtout); + +getpar("moment","f",&moment); +getpar("tstart","f",&tstart); + +getpar("tshift_timedomain","d",&tshift_timedomain); + +endpar(); + +maxmech = 1; +mechpar.nmech = 1; +mechpar.flag[0] = U1FLAG; +mechpar.flag[1] = 0; +mechpar.flag[2] = 0; + +if(srf_flag == 1 || mrf_flag == 1) + { + if(srf_flag == 1) + maxmech = 3; + if(mrf_flag == 1) + mechpar.flag[0] = MOMTEN_FLAG; + + read_srf(&srf,rupmodfile,inbin); + prect_ptr = &srf.srf_prect; + prseg_ptr = prect_ptr->prectseg; + apnts_ptr = &srf.srf_apnts; + apval_ptr = apnts_ptr->apntvals; + +/* +03/03/06 + + Default is to use all POINTS specified in SRF file. This is + done with nseg=-1. Values for 'shypo', 'dhypo', 'len', 'wid' + are not important. + + -or- + + Only use one segment from standard rupture model format; + specified with 'nseg'. +*/ + + if(nseg < 0) /* do all POINTS */ + { + np = srf.srf_apnts.np; + + apv_off = 0; + } + else + { + np = prseg_ptr[nseg].nstk*prseg_ptr[nseg].ndip; + + /* reset POINTS pointer to appropriate segment */ + + apv_off = 0; + for(i=0;i 1.0001) + tstart = tstart - apval_ptr[0].dt*((int)(10.0*dtout/(apval_ptr[0].dt) + 0.5)); + +gf = (float *) check_malloc (4*gfpar.nc*ntsum*sizeof(float)); +gfmech = (float *) check_malloc (maxmech*12*ntsum*sizeof(float)); +space = (float *) check_malloc (2*ntsum*sizeof(float)); + +seis = (float *) check_malloc (3*ntout*sizeof(float)); +subseis = (float *) check_malloc (maxmech*3*ntout*sizeof(float)); +stf = (float *) check_malloc (ntout*sizeof(float)); + +if(statfile[0] != '\0') + { + stat_buf = (char *) check_malloc (NCHAR_STAT*nbuf*sizeof(char)); + slon_buf = (float *) check_malloc (nbuf*sizeof(float)); + slat_buf = (float *) check_malloc (nbuf*sizeof(float)); + + fpr = fopfile(statfile,"r"); + + nstat = 0; + while(fgets(string,2048,fpr) != NULL) + { + if(strncmp(string,"#",1) != 0) + { + stat_ptr = stat_buf + NCHAR_STAT*nstat; + sscanf(string,"%f %f %s",&slon_buf[nstat],&slat_buf[nstat],stat_ptr); + nstat++; + + if(nstat % nbuf == 0) + { + stat_buf = (char *) check_realloc (stat_buf,NCHAR_STAT*(nstat+nbuf)*sizeof(char)); + slon_buf = (float *) check_realloc (slon_buf,(nstat+nbuf)*sizeof(float)); + slat_buf = (float *) check_realloc (slat_buf,(nstat+nbuf)*sizeof(float)); + } + } + } + + stat_buf = (char *) check_realloc (stat_buf,NCHAR_STAT*nstat*sizeof(char)); + slon_buf = (float *) check_realloc (slon_buf,nstat*sizeof(float)); + slat_buf = (float *) check_realloc (slat_buf,nstat*sizeof(float)); + } +else if(latloncoords == 1) + { + nstat = 1; + stat_buf = (char *) check_malloc (NCHAR_STAT*sizeof(char)); + slon_buf = (float *) check_malloc (sizeof(float)); + slat_buf = (float *) check_malloc (sizeof(float)); + + strcpy(stat_buf,stat); + slon_buf[0] = slon; + slat_buf[0] = slat; + } + +for(istat=0;istat (float)(0.0)) + { + get_ard_srf(&srf,apv_off,ip,&azi,&rng,&z0,&deast,&dnorth,&slon_buf[istat],&slat_buf[istat],&geop); + + find_4gf_v2(gfpar,gfhead,&rng,&z0,&deast,&dnorth,&maxgft); + + if(use_closest_gf == 1 && rng > min_taper_range) + reset_wt_4gf_v2(gfhead,&rng,&min_taper_range,&max_taper_range,&tadjust); + + read_4gf_v2(gfpath,gfname,gf,ntsum,gfhead,gfpar,&maxnt,&dtout,space); + + if(mrf_flag == 1) /* moment-tensor source */ + { + scale = 1.0; + momten_4gf(&srf,apv_off,ip,gfmech,gf,gfhead,gfpar,ntsum,&azi); + } + else /* double-couple source specified by slip */ + { + scale = normf*apval_ptr[ip].area; + mech_4gf_v2(gfmech,gf,gfhead,gfpar,ntsum,mechpar,&azi,&scale); + } + tmom = tmom + src_amp*scale; + + sum_4gf_v2(subseis,ntout,gfmech,gfhead,ntsum,maxnt,&rt,&maxgft,&tstart,tshift_timedomain,mechpar,&tadjust); + + srf_stf_v2(&srf,apv_off,ip,seis,subseis,stf,ntout,&dtout,mechpar,space); + } + } + + sv = seis; + sn = seis + ntout; + se = seis + 2*ntout; + + strncpy(sname,stat_ptr,7); + sname[7] = '\0'; + + fprintf(stderr,"Total summed moment= %13.5e\n",tmom); + + if(moment > 0.0) + { + sfac = moment/tmom; + for(it=0;itzsgt = (int)((double)((*dep))*invh + 1.5); indx->indx = (long long)(indx->xsgt)*(long long)(100000000) + (long long)(indx->ysgt)*(long long)(10000) + (long long)(indx->zsgt); } +void get_indx_pow(float *lon,float *lat,float *dep,struct sgtindex *indx,struct geoprojection *gp,int ipow) +{ +float xs, ys, xr, yr; +double invh; +long long indx_shft; + +indx_shft = 1; +while(ipow > 0) + { + indx_shft = indx_shft*10; + ipow--; + } + +invh = 1.0/(double)(indx->h); + +if(gp->geoproj == 0) + { + xs = ((*lon) - gp->modellon)*gp->kmlon; + ys = (gp->modellat - (*lat))*gp->kmlat; + + xr = xs*gp->cosR + ys*gp->sinR - gp->xshift; + yr = -xs*gp->sinR + ys*gp->cosR - gp->yshift; + } +else if(gp->geoproj == 1) + gcproj(&xr,&yr,lon,lat,&gp->erad,&gp->g0,&gp->b0,gp->amat,gp->ainv,1); + +indx->xsgt = (int)((double)(xr)*invh + 0.5); +indx->ysgt = (int)((double)(yr)*invh + 0.5); +indx->zsgt = (int)((double)((*dep))*invh + 1.5); +indx->indx = (long long)(indx->xsgt)*(long long)(indx_shft)*(long long)(indx_shft) + (long long)(indx->ysgt)*(long long)(indx_shft) + (long long)(indx->zsgt); +} + void get_ard_srf(struct standrupformat *srf,int off,int ip,float *az,float *rg,float *z0,float *de,float *dn,float *slon,float *slat,struct geoprojection *gp) { struct srf_planerectangle *prect_ptr; @@ -133,9 +165,7 @@ else if(gp->geoproj == 1) *z0 = apval_ptr[ip].dep; } -void resample(s,nt,dt,isamp,ntpad,ntrsmp,newdt,p) -float *s, *p, *dt, *newdt; -int nt, isamp, ntpad, ntrsmp; +void resample(float *s,int nt,float *dt,int isamp,int ntpad,int ntrsmp,float *newdt,float *p) { float df, f, f0, fl, fl2, fac; int i, j; @@ -210,9 +240,7 @@ while(n--) } } -void taper_norm(g,dt,nt) -float *g, *dt; -int nt; +void taper_norm(float *g,float *dt,int nt) { float fac, df, arg; int i; @@ -232,9 +260,7 @@ for(i=nt-ntap;isrf_prect); prseg_ptr = prect_ptr->prectseg; apnts_ptr = &(srf->srf_apnts); @@ -173,25 +175,44 @@ mpar->flag[0] = 0; mpar->flag[1] = 0; mpar->flag[2] = 0; -if(apval_ptr[ip].nt1 > 0) - { - mpar->flag[mpar->nmech] = U1FLAG; - mpar->nmech = mpar->nmech + 1; - } -if(apval_ptr[ip].nt2 > 0) +mrf_flag = 0; +if(strcmp(srf[0].src_format,"MOMENT") == 0) + mrf_flag = 1; + +if(mrf_flag == 1) { - mpar->flag[mpar->nmech] = U2FLAG; - mpar->nmech = mpar->nmech + 1; + mpar->flag[0] = MOMTEN_FLAG; + mpar->nmech = 1; + + *samp = sqrt(0.5*apval_ptr[ip].mnn*apval_ptr[ip].mnn + + 0.5*apval_ptr[ip].mee*apval_ptr[ip].mee + + 0.5*apval_ptr[ip].mdd*apval_ptr[ip].mdd + + apval_ptr[ip].mne*apval_ptr[ip].mne + + apval_ptr[ip].mnd*apval_ptr[ip].mnd + + apval_ptr[ip].med*apval_ptr[ip].med); } -if(apval_ptr[ip].nt3 > 0) +else { - mpar->flag[mpar->nmech] = U3FLAG; - mpar->nmech = mpar->nmech + 1; - } + if(apval_ptr[ip].nt1 > 0) + { + mpar->flag[mpar->nmech] = U1FLAG; + mpar->nmech = mpar->nmech + 1; + } + if(apval_ptr[ip].nt2 > 0) + { + mpar->flag[mpar->nmech] = U2FLAG; + mpar->nmech = mpar->nmech + 1; + } + if(apval_ptr[ip].nt3 > 0) + { + mpar->flag[mpar->nmech] = U3FLAG; + mpar->nmech = mpar->nmech + 1; + } -*vs = sqrt(apval_ptr[ip].slip1*apval_ptr[ip].slip1 + *samp = sqrt(apval_ptr[ip].slip1*apval_ptr[ip].slip1 + apval_ptr[ip].slip2*apval_ptr[ip].slip2 + apval_ptr[ip].slip3*apval_ptr[ip].slip3); + } mpar->stk = apval_ptr[ip].stk; mpar->dip = apval_ptr[ip].dip; @@ -199,3 +220,222 @@ mpar->rak = apval_ptr[ip].rake; *rt = apval_ptr[ip].tinit; } + +void srf_stf_v2(struct standrupformat *srf,int off,int ip,float *s,float *u,float *stf,int nt,float *dt,struct mechparam mp,float *unused) +{ +FILE *fpw; +int it, nstf, im; +float sum, *sptr1, *uptr; +float *space, *sptr2; + +struct srf_planerectangle *prect_ptr; +struct srf_prectsegments *prseg_ptr; +struct srf_allpoints *apnts_ptr; +struct srf_apointvalues *apval_ptr; + +float fnt; +int resamp, ntpad, ntrsmp, gnt; + +float slip, sfac, dt_stf, dt_tmp, ds_dt; +int nt_fac, nt_tmp, jj, kk; + +float tol = 1.0e-02; + +float pratio_tol, mratio_tol; +float ratio_tol = 0.001; + +pratio_tol = 1.0 + ratio_tol; +mratio_tol = 1.0 - ratio_tol; + +prect_ptr = &(srf->srf_prect); +prseg_ptr = prect_ptr->prectseg; +apnts_ptr = &(srf->srf_apnts); +apval_ptr = apnts_ptr->apntvals + off; + +zapit(stf,nt); + +sptr1 = NULL; +space = NULL; +sptr2 = NULL; + +for(im=0;im 0) + { + dt_stf = apval_ptr[ip].dt; + +/* resample if needed */ + + if((*dt)/dt_stf > pratio_tol || (*dt)/dt_stf < mratio_tol) + { + if(nstf == 1) /* add zero points in front and back for this case */ + { + nstf = nstf + 2; + sptr1 = (float *)check_realloc((void *)sptr1,nstf*sizeof(float)); + + sptr1[2] = 0.0; + sptr1[1] = sptr1[0]; + sptr1[0] = 0.0; + } + else /* add zero point in back for this case */ + { + nstf++; + sptr1 = (float *)check_realloc((void *)sptr1,nstf*sizeof(float)); + + sptr1[nstf-1] = 0.0; + } + +/* first try linear interpolation as this best mimics target slip-rate */ + + nt_fac = (int)(dt_stf/(*dt) + 0.5); + if(nt_fac >= 2) /* otherwise, skip linear interpolation */ + { + nt_tmp = nt_fac*(nstf-1) + 1; + dt_tmp = dt_stf/nt_fac; + sptr2 = (float *)check_realloc((void *)sptr2,nt_tmp*sizeof(float)); + + for(it=0;it pratio_tol) /* add some initial padding to prevent artifacts */ + { + ntpad = (int)(10.0*(*dt)/dt_stf + 0.5); + + sptr1 = (float *) check_realloc((void *)sptr1,(2*ntpad+nstf)*sizeof(float)); + + for(it=(nstf)-1;it>=0;it--) + sptr1[it+ntpad] = sptr1[it]; + + for(it=0;it pratio_tol || (*dt)/dt_stf < mratio_tol) + { + ntpad = 2*nstf; + fnt = ntpad*dt_stf/(*dt); + gnt = (int)(fnt + 0.5); + + while(nt_tol(fnt,gnt) > tol) + { + ntpad++; + fnt = ntpad*dt_stf/(*dt); + gnt = (int)(fnt + 0.5); + } + + ntrsmp = (int)(fnt); + + if((*dt) < dt_stf) + { + space = (float *) check_realloc ((void *)space,2*ntrsmp*sizeof(float)); + sptr2 = (float *)check_realloc((void *)sptr2,2*ntrsmp*sizeof(float)); + + for(it=0;it<2*ntrsmp;it++) + sptr2[it] = 0.0; + + resamp = 1; + } + else + { + space = (float *) check_realloc ((void *)space,2*ntpad*sizeof(float)); + sptr2 = (float *)check_realloc((void *)sptr2,2*ntpad*sizeof(float)); + + for(it=0;it<2*ntpad;it++) + sptr2[it] = 0.0; + + resamp = -1; + } + + for(it=0;it 0.0) + { + sfac = 0.0; + for(it=0;it nt) + nstf = nt; + +/* add factor of dt to prenormalize convolution */ + for(it=0;it= 3.0) + { + ip = sscanf(str,"%*s %s",srf[0].src_format); + + if(ip != 1 || (strcmp(srf[0].src_format,"SLIP") != 0 && strcmp(srf[0].src_format,"MOMENT") != 0)) + sprintf(srf[0].src_format,"SLIP"); + + if(strcmp(srf[0].src_format,"MOMENT") == 0) + mrf_flag = 1; + } + /* reserve first 3 lines for writing out command that generated SRF */ srf[0].srf_hcmnt.nline = 3; srf[0].srf_hcmnt.cbuf = (char *)check_malloc((srf[0].srf_hcmnt.nline)*MAXLINE*sizeof(char)); @@ -219,10 +234,11 @@ else &(apval_ptr->tinit), &(apval_ptr->dt)); + apval_ptr->vp = -1; apval_ptr->vs = -1; apval_ptr->den = -1; } - else if(atof(srf->version) >= 2.0) + else if(atof(srf->version) < 3.0) { sscanf(str,"%f %f %f %f %f %f %f %f %f %f",&(apval_ptr->lon), &(apval_ptr->lat), @@ -234,50 +250,96 @@ else &(apval_ptr->dt), &(apval_ptr->vs), &(apval_ptr->den)); + apval_ptr->vp = -1; } + else if(atof(srf->version) >= 2.0) + { + sscanf(str,"%f %f %f %f %f %f %f %f %f %f %f",&(apval_ptr->lon), + &(apval_ptr->lat), + &(apval_ptr->dep), + &(apval_ptr->stk), + &(apval_ptr->dip), + &(apval_ptr->area), + &(apval_ptr->tinit), + &(apval_ptr->dt), + &(apval_ptr->vp), + &(apval_ptr->vs), + &(apval_ptr->den)); + } + + if(mrf_flag != 1) /* expecting slip */ + { + fgets(str,MAXLINE,fpr); + sscanf(str,"%f %f %d %f %d %f %d",&(apval_ptr->rake), + &(apval_ptr->slip1), + &(apval_ptr->nt1), + &(apval_ptr->slip2), + &(apval_ptr->nt2), + &(apval_ptr->slip3), + &(apval_ptr->nt3)); - fgets(str,MAXLINE,fpr); - sscanf(str,"%f %f %d %f %d %f %d",&(apval_ptr->rake), - &(apval_ptr->slip1), - &(apval_ptr->nt1), - &(apval_ptr->slip2), - &(apval_ptr->nt2), - &(apval_ptr->slip3), - &(apval_ptr->nt3)); + if(apval_ptr->nt1) + apval_ptr->stf1 = (float *)check_malloc((apval_ptr->nt1)*sizeof(float)); + else + apval_ptr->stf1 = NULL; - if(apval_ptr->nt1) - apval_ptr->stf1 = (float *)check_malloc((apval_ptr->nt1)*sizeof(float)); - else - apval_ptr->stf1 = NULL; + stf = apval_ptr->stf1; - stf = apval_ptr->stf1; + for(it=0;it<(apval_ptr->nt1);it++) + fscanf(fpr,"%f",&stf[it]); - for(it=0;it<(apval_ptr->nt1);it++) - fscanf(fpr,"%f",&stf[it]); + if(apval_ptr->nt2) + apval_ptr->stf2 = (float *)check_malloc((apval_ptr->nt2)*sizeof(float)); + else + apval_ptr->stf2 = NULL; - if(apval_ptr->nt2) - apval_ptr->stf2 = (float *)check_malloc((apval_ptr->nt2)*sizeof(float)); - else - apval_ptr->stf2 = NULL; + stf = apval_ptr->stf2; - stf = apval_ptr->stf2; + for(it=0;it<(apval_ptr->nt2);it++) + fscanf(fpr,"%f",&stf[it]); - for(it=0;it<(apval_ptr->nt2);it++) - fscanf(fpr,"%f",&stf[it]); + if(apval_ptr->nt3) + apval_ptr->stf3 = (float *)check_malloc((apval_ptr->nt3)*sizeof(float)); + else + apval_ptr->stf3 = NULL; - if(apval_ptr->nt3) - apval_ptr->stf3 = (float *)check_malloc((apval_ptr->nt3)*sizeof(float)); - else - apval_ptr->stf3 = NULL; + stf = apval_ptr->stf3; - stf = apval_ptr->stf3; + for(it=0;it<(apval_ptr->nt3);it++) + fscanf(fpr,"%f",&stf[it]); - for(it=0;it<(apval_ptr->nt3);it++) - fscanf(fpr,"%f",&stf[it]); + /* get rouge newline character */ + if((apval_ptr->nt1) || (apval_ptr->nt2) || (apval_ptr->nt3)) + fgets(str,MAXLINE,fpr); + } - /* get rouge newline character */ - if((apval_ptr->nt1) || (apval_ptr->nt2) || (apval_ptr->nt3)) + else if(mrf_flag == 1) /* expecting moment */ + { fgets(str,MAXLINE,fpr); + sscanf(str,"%lf %lf %lf %lf %lf %lf %d", + &(apval_ptr->mnn), + &(apval_ptr->mee), + &(apval_ptr->mdd), + &(apval_ptr->mne), + &(apval_ptr->mnd), + &(apval_ptr->med), + &(apval_ptr->ntmr)); + + if(apval_ptr->ntmr) + apval_ptr->mrf = (float *)check_malloc((apval_ptr->ntmr)*sizeof(float)); + else + apval_ptr->mrf = NULL; + + stf = apval_ptr->mrf; + + for(it=0;it<(apval_ptr->ntmr);it++) + fscanf(fpr,"%f",&stf[it]); + + /* get rouge newline character */ + if((apval_ptr->ntmr)) + fgets(str,MAXLINE,fpr); + } + } } @@ -317,8 +379,10 @@ void write_srf(struct standrupformat *srf,char *file,int bflag) { if(atof(srf->version) < 2.0) write_srf1(srf,file,bflag); -else if(atof(srf->version) >= 2.0) +else if(atof(srf->version) < 3.0) write_srf2(srf,file,bflag); +else if(atof(srf->version) >= 3.0) + write_srf3(srf,file,bflag); } void write_srf1(struct standrupformat *srf,char *file,int bflag) @@ -693,6 +757,245 @@ else } } +void write_srf3(struct standrupformat *srf,char *file,int bflag) +{ +FILE *fpw, *fopfile(); +struct srf_planerectangle *prect_ptr; +struct srf_prectsegments *prseg_ptr; +struct srf_allpoints *apnts_ptr; +struct srf_apointvalues *apval_ptr; + +float area; +float *stf; +int i, j, k, nt6, it, ig; + +char *sptr, pword[32]; +int fdw; + +int ip, nprite; + +int mrf_flag; + +prect_ptr = &(srf->srf_prect); +prseg_ptr = prect_ptr->prectseg; +apnts_ptr = &(srf->srf_apnts); +apval_ptr = apnts_ptr->apntvals; + +if(bflag) + { + if(strcmp(file,"stdout") == 0) + fdw = STDOUT_FILENO; + else + fdw = croptrfile(file); + + rite(fdw,srf->version,sizeof(srf->version)); + + if(strcmp(srf->type,"PLANE") == 0) + { + rite(fdw,srf->type,sizeof(srf->type)); + rite(fdw,&(prect_ptr->nseg),sizeof(prect_ptr->nseg)); + rite(fdw,prseg_ptr,(prect_ptr->nseg)*sizeof(struct srf_prectsegments)); + } + + sprintf(pword,"POINTS"); + rite(fdw,pword,sizeof(pword)); + rite(fdw,&(apnts_ptr->np),sizeof(apnts_ptr->np)); + for(i=0;inp;i++) + { + rite(fdw,&(apval_ptr[i].lon),sizeof(float)); + rite(fdw,&(apval_ptr[i].lat),sizeof(float)); + rite(fdw,&(apval_ptr[i].dep),sizeof(float)); + rite(fdw,&(apval_ptr[i].stk),sizeof(float)); + rite(fdw,&(apval_ptr[i].dip),sizeof(float)); + rite(fdw,&(apval_ptr[i].area),sizeof(float)); + rite(fdw,&(apval_ptr[i].tinit),sizeof(float)); + rite(fdw,&(apval_ptr[i].dt),sizeof(float)); + rite(fdw,&(apval_ptr[i].rake),sizeof(float)); + rite(fdw,&(apval_ptr[i].slip1),sizeof(float)); + rite(fdw,&(apval_ptr[i].nt1),sizeof(int)); + rite(fdw,&(apval_ptr[i].slip2),sizeof(float)); + rite(fdw,&(apval_ptr[i].nt2),sizeof(int)); + rite(fdw,&(apval_ptr[i].slip3),sizeof(float)); + rite(fdw,&(apval_ptr[i].nt3),sizeof(int)); + + rite(fdw,apval_ptr[i].stf1,(apval_ptr[i].nt1)*sizeof(float)); + rite(fdw,apval_ptr[i].stf2,(apval_ptr[i].nt2)*sizeof(float)); + rite(fdw,apval_ptr[i].stf3,(apval_ptr[i].nt3)*sizeof(float)); + } + close(fdw); + } +else + { + if(strcmp(file,"stdout") == 0) + fpw = stdout; + else + fpw = fopfile(file,"w"); + + fprintf(fpw,"%s %s\n",srf->version,srf->src_format); + + for(i=0;isrf_hcmnt.nline;i++) + { + sptr = (srf->srf_hcmnt.cbuf) + i*MAXLINE; + if(sptr[0] != '\0') + fprintf(fpw,"%s\n",sptr); + } + + if(strcmp(srf->type,"PLANE") == 0) + { + fprintf(fpw,"%s %d\n",srf->type,prect_ptr->nseg); + for(ig=0;ignseg;ig++) + { + fprintf(fpw,"%12.6f %11.6f %5d %5d %10.4f %10.4f\n",prseg_ptr[ig].elon, + prseg_ptr[ig].elat, + prseg_ptr[ig].nstk, + prseg_ptr[ig].ndip, + prseg_ptr[ig].flen, + prseg_ptr[ig].fwid); + fprintf(fpw,"%4.0f %4.0f %10.4f %10.4f %10.4f\n",prseg_ptr[ig].stk, + prseg_ptr[ig].dip, + prseg_ptr[ig].dtop, + prseg_ptr[ig].shyp, + prseg_ptr[ig].dhyp); + } + } + + mrf_flag = 0; + if(strcmp(srf->src_format,"MOMENT") == 0) + mrf_flag = 1; + + nprite = 0; + for(ig=0;ignseg;ig++) + { + fprintf(fpw,"POINTS %d\n",srf->np_seg[ig]); + for(ip=0;ipnp_seg[ig];ip++) + { + i = ip + nprite; + + fprintf(fpw,"%12.6f %11.6f %12.5e %4.0f %4.0f %12.5e %13.6e %12.5e %13.5e %13.5e %13.5e\n", + apval_ptr[i].lon, + apval_ptr[i].lat, + apval_ptr[i].dep, + apval_ptr[i].stk, + apval_ptr[i].dip, + apval_ptr[i].area, + apval_ptr[i].tinit, + apval_ptr[i].dt, + apval_ptr[i].vp, + apval_ptr[i].vs, + apval_ptr[i].den); + + if(mrf_flag != 1) /* expecting slip */ + { + fprintf(fpw,"%4.0f %10.4f %6d %10.4f %6d %10.4f %6d\n", + apval_ptr[i].rake, + apval_ptr[i].slip1, + apval_ptr[i].nt1, + apval_ptr[i].slip2, + apval_ptr[i].nt2, + apval_ptr[i].slip3, + apval_ptr[i].nt3); + + stf = apval_ptr[i].stf1; + nt6 = (apval_ptr[i].nt1)/6; + for(k=0;knp_seg[ig]; + } + + fclose(fpw); + } +} + void free_srf_stf(struct standrupformat *srf) { struct srf_allpoints *apnts_ptr; @@ -1648,3 +1951,215 @@ nb = sprintf(sptr1,"# ending_seed= %lld",ending_seed); sptr1 = sptr1 + MAXLINE; sprintf(sptr1,"#"); } + +void srf_to_mrf(struct standrupformat *srf,struct standrupformat *mrf,struct velmodel *vm,int use_srf_lame,int pflag,int ac,char **av) +{ +struct srf_prectsegments *prseg_in, *prseg_out; +struct srf_apointvalues *apval_in, *apval_out; +float *stfin, *stfout; +int i, j, k, it, ig; + +double s_mom, m_mom; + +float u1, u2, u3, sum; +double lam, l2m, mu; +double ux, uy, uz, vx, vy, vz; + +double arg; +double cosS, sinS; +double cosD, sinD; +double cosL, sinL; + +double rperd = 0.017453292519943; + +if(atof(srf[0].version) < 2.0) + { + fprintf(stderr,"srf version= %s < 2.0, exiting ... \n",srf[0].version); + exit(-1); + } + +s_mom = 0.0; +m_mom = 0.0; + +sprintf(mrf[0].version,"3.0"); +sprintf(mrf[0].src_format,"MOMENT"); + +copy_hcmnt(mrf,srf); + +if(pflag && atof(mrf[0].version) >= 2.0) + load_command_srf(mrf,ac,av); + +mrf[0].type[0] = '\0'; +if(strncmp(srf[0].type,"PLANE",5) == 0) + { + strcpy(mrf[0].type,srf[0].type); + + mrf[0].srf_prect.nseg = srf[0].srf_prect.nseg; + mrf[0].srf_prect.prectseg = (struct srf_prectsegments *)check_malloc(mrf[0].srf_prect.nseg*sizeof(struct srf_prectsegments)); + + prseg_in = srf[0].srf_prect.prectseg; + prseg_out = mrf[0].srf_prect.prectseg; + for(ig=0;ig apval_out[i].ntmr) + apval_out[i].ntmr = apval_in[i].nt2; + if(apval_in[i].nt3 > apval_out[i].ntmr) + apval_out[i].ntmr = apval_in[i].nt3; + + apval_out[i].mrf = (float *)check_realloc(apval_out[i].mrf,(apval_out[i].ntmr)*sizeof(float)); + stfout = apval_out[i].mrf; + + if(apval_in[i].nt1) + { + stfin = apval_in[i].stf1; + for(it=0;it<(apval_in[i].nt1);it++) + stfout[it] = stfin[it]*stfin[it]; + } + + if(apval_in[i].nt2) + { + stfin = apval_in[i].stf2; + for(it=0;it<(apval_in[i].nt2);it++) + stfout[it] = stfout[it] + stfin[it]*stfin[it]; + } + + if(apval_in[i].nt3) + { + stfin = apval_in[i].stf3; + for(it=0;it<(apval_in[i].nt3);it++) + stfout[it] = stfout[it] + stfin[it]*stfin[it]; + } + + sum = 0.0; + for(it=0;it<(apval_out[i].ntmr);it++) + { + stfout[it] = sqrt(stfout[it]); + sum = sum + (apval_out[i].dt)*stfout[it]; + } + + sum = 1.0/sum; + for(it=0;it<(apval_out[i].ntmr);it++) + stfout[it] = sum*stfout[it]; + + if(use_srf_lame != 0 && apval_out[i].vp > 0.0 && apval_out[i].vs > 0.0 && apval_out[i].den > 0.0) + { + l2m = apval_out[i].vp*apval_out[i].vp*apval_out[i].den; + mu = apval_out[i].vs*apval_out[i].vs*apval_out[i].den; + } + else + { + j = 0; + while(vm->dep[j] < apval_out[i].dep) + j++; + + l2m = 1.0e+10*vm->vp[j]*vm->vp[j]*vm->den[j]; + mu = 1.0e+10*vm->vs[j]*vm->vs[j]*vm->den[j]; + + if(apval_out[i].vp < 0.0) + apval_out[i].vp = 1.0e+05*vm->vp[j]; + + if(apval_out[i].vs < 0.0) + apval_out[i].vs = 1.0e+05*vm->vs[j]; + + if(apval_out[i].den < 0.0) + apval_out[i].den = vm->den[j]; + } + lam = l2m - 2.0*mu; + + u1 = apval_in[i].slip1; + u2 = apval_in[i].slip2; + u3 = apval_in[i].slip3; + + arg = apval_in[i].stk*rperd; + cosS = cos(arg); + sinS = sin(arg); + + arg = apval_in[i].dip*rperd; + cosD = cos(arg); + sinD = sin(arg); + + arg = apval_in[i].rake*rperd; + cosL = cos(arg); + sinL = sin(arg); + + vx = -sinD*sinS; + vy = sinD*cosS; + vz = -cosD; + + ux = -(u3*sinD - cosD*(u1*sinL + u2*cosL))*sinS + (u1*cosL - u2*sinL)*cosS; + uy = (u3*sinD - cosD*(u1*sinL + u2*cosL))*cosS + (u1*cosL - u2*sinL)*sinS; + uz = -u3*cosD - (u1*sinL + u2*cosL)*sinD; + +/* XXXXX +fprintf(stderr,"vx= %.5e vy= %.5e vz= %.5e\n",vx,vy,vz); +fprintf(stderr,"cosL= %.5e sinL= %.5e cosD= %.5e sinD= %.5e\n",cosL,sinL,cosD,sinD); +fprintf(stderr,"u1= %.5e u2= %.5e u3= %.5e\n",u1,u2,u3); +fprintf(stderr,"ux= %.5e uy= %.5e uz= %.5e\n",ux,uy,uz); +*/ + + apval_out[i].mnn = (l2m*vx*ux + lam*(vy*uy + vz*uz))*apval_out[i].area; + apval_out[i].mee = (l2m*vy*uy + lam*(vx*ux + vz*uz))*apval_out[i].area; + apval_out[i].mdd = (l2m*vz*uz + lam*(vx*ux + vy*uy))*apval_out[i].area; + + apval_out[i].mne = mu*(vx*uy + vy*ux)*apval_out[i].area; + apval_out[i].mnd = mu*(vx*uz + vz*ux)*apval_out[i].area; + apval_out[i].med = mu*(vy*uz + vz*uy)*apval_out[i].area; + + s_mom = s_mom + sqrt(u1*u1 + u2*u2 + u3*u3)*mu*apval_in[i].area; + m_mom = m_mom + sqrt(0.5*apval_out[i].mnn*apval_out[i].mnn + + 0.5*apval_out[i].mee*apval_out[i].mee + + 0.5*apval_out[i].mdd*apval_out[i].mdd + + apval_out[i].mne*apval_out[i].mne + + apval_out[i].mnd*apval_out[i].mnd + + apval_out[i].med*apval_out[i].med); + + } +fprintf(stderr,"slip moment= %.5e\n",s_mom); +fprintf(stderr,"MT moment= %.5e\n",m_mom); +} diff --git a/bbp/src/gp/StandRupFormat/structure.h b/bbp/src/gp/StandRupFormat/structure.h index d022a323..2e85ccb6 100644 --- a/bbp/src/gp/StandRupFormat/structure.h +++ b/bbp/src/gp/StandRupFormat/structure.h @@ -83,7 +83,7 @@ struct srf_apointvalues20141109 float *stf3; }; -struct srf_apointvalues +struct srf_apointvalues20241010 { float lon; float lat; @@ -107,6 +107,39 @@ struct srf_apointvalues float *stf3; }; +struct srf_apointvalues + { + float lon; + float lat; + float dep; + float stk; + float dip; + float area; + float tinit; + float dt; + float vp; /* added for V3.0 */ + float vs; + float den; + float rake; + float slip1; + int nt1; + float slip2; + int nt2; + float slip3; + int nt3; + float *stf1; + float *stf2; + float *stf3; + double mnn; /* moment-rate stuff added for V3.0 */ + double mee; + double mdd; + double mne; + double mnd; + double med; + int ntmr; + float *mrf; + }; + struct srf_allpoints { int np; @@ -150,7 +183,7 @@ struct standrupformat20141109 struct srf_allpoints srf_apnts; }; -struct standrupformat +struct standrupformat20241010 { char version[32]; char type[32]; @@ -161,6 +194,18 @@ struct standrupformat struct srf_allpoints srf_apnts; }; +struct standrupformat + { + char version[32]; + char type[32]; + char src_format[64]; /* added for V3.0, options are "SLIP" or "MOMENT", default is "SLIP" */ + int nseg; + int *np_seg; + struct srf_headercomment srf_hcmnt; + struct srf_planerectangle srf_prect; + struct srf_allpoints srf_apnts; + }; + struct slippars20171024 { float lon; @@ -197,6 +242,8 @@ struct slippars float den; float rt_s; float rt_e; + float aseis; + float rvf; }; struct segpars