Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
strdepth='%.4f' % zs
subfault=str(int(source[0])).rjust(4,'0')
if static==0 and tsunami==0: #Where to save dynamic waveforms
green_path=home+project_name+'/GFs/dynamic/'+model_name+"_"+strdepth+".sub"+subfault+"/"
if static==1 and tsunami==1: #Where to save dynamic waveforms
green_path=home+project_name+'/GFs/tsunami/'+model_name+"_"+strdepth+".sub"+subfault+"/"
if static==1 and tsunami==0: #Where to save statics
green_path=home+project_name+'/GFs/static/'+model_name+"_"+strdepth+".sub"+subfault+"/"
staname=genfromtxt(station_file,dtype="U",usecols=0)
if staname.shape==(): #Single staiton file
staname=array([staname])
#Compute distances and azimuths
d,az,lon_sta,lat_sta=src2sta(station_file,source,output_coordinates=True)
#Get moment corresponding to 1 meter of slip on subfault
mu=get_mu(structure,zs)
Mo=mu*ss_length*ds_length*1.0
Mw=(2./3)*(log10(Mo)-9.1)
#Move to output folder
os.chdir(green_path)
print('Processor '+str(rank)+' is working on subfault '+str(int(source[0]))+' and '+str(len(d))+' stations ')
for k in range(len(d)):
if static==0: #Compute full waveforms
diststr='%.6f' % d[k] #Need current distance in string form for external call
#Form the strings to be used for the system calls according to user desired options
if integrate==1: #Make displ.
#First Stike-Slip GFs
if custom_stf==None:
commandSS="syn -I -M"+str(Mw)+"/"+str(strike)+"/"+str(dip)+"/"+str(rakeSS)+" -D"+str(duration)+ \
"/"+str(rise)+" -A"+str(az[k])+" -O"+staname[k]+".subfault"+num+".SS.disp.x -G"+green_path+diststr+".grn.0"
commandSS=split(commandSS) #Split string into lexical components for system call
ds_length=source[9]
ss_length_in_km=ss_length/1000.
ds_length_in_km=ds_length/1000.
strdepth='%.4f' % zs
if static==0 and tsunami==0: #Where to save dynamic waveforms
green_path=green_path+'dynamic/'+model_name+"_"+strdepth+".sub"+subfault+"/"
if static==0 and tsunami==1: #Where to save dynamic waveforms
green_path=green_path+'tsunami/'+model_name+"_"+strdepth+".sub"+subfault+"/"
print("--> Computing synthetics at stations for the source at ("+str(xs)+" , "+str(ys)+")")
staname=genfromtxt(station_file,dtype="U",usecols=0)
if staname.shape==(): #Single staiton file
staname=array([staname])
#Compute distances and azimuths
d,az,lon_sta,lat_sta=src2sta(station_file,source,output_coordinates=True)
#Get moment corresponding to 1 meter of slip on subfault
mu=get_mu(structure,zs)
Mo=mu*ss_length*ds_length*1
Mw=(2./3)*(log10(Mo)-9.1)
#Move to output folder
log='' #Initalize log
os.chdir(green_path)
for k in range(len(d)):
if static==0: #Compute full waveforms
diststr='%.3f' % d[k] #Need current distance in string form for external call
#Form the strings to be used for the system calls according to suer desired options
if integrate==1: #Make displ.
#First Stike-Slip GFs
commandSS="syn -I -M"+str(Mw)+"/"+str(strike)+"/"+str(dip)+"/"+str(rakeSS)+" -D"+str(duration)+ \
"/"+str(rise)+" -A"+str(az[k])+" -O"+staname[k]+".subfault"+num+".SS.disp.x -G"+green_path+diststr+".grn.0"
print(commandSS) #Output to screen so I know we're underway
log=log+commandSS+'\n' #Append to log
commandSS=split(commandSS) #Split string into lexical components for system call
# are there forced onset times?
if onset_file is not None:
tdelay_constant = genfromtxt(home+project_name+'/data/model_info/'+onset_file)
#Get slip quantities
iss=2*arange(len(f)*num_windows)
ids=2*arange(len(f)*num_windows)+1
ss=sol[iss]
ds=sol[ids]
#Get rigidities
mu=zeros(len(ds))
trup=zeros(len(ds))
j=0
for krup in range(num_windows):
for k in range(len(f)):
mu[j]=get_mu(mod,f[k,3])
j+=1
#Get rupture start times
trupt=arange(0,num_windows)*trise/2 #Time delays fore ach sub-window
for krup in range(num_windows):
if onset_file==None: #No forced onset times
trup[krup*(len(ds)//num_windows):(krup+1)*(len(ds)//num_windows)]=epi2subfault(epicenter,f,rupture_speed,trupt[krup])
else: # use specified onset times plus window time delay
trup[krup*(len(ds)//num_windows):(krup+1)*(len(ds)//num_windows)]=tdelay_constant+trupt[krup]
#Prepare for output
out1=f[:,0:8]
out2=f[:,8:10]
for k in range(num_windows-1):
out1=r_[out1,f[:,0:8]]
out2=r_[out2,f[:,8:10]]
out=c_[out1,ss,ds,out2,trup,mu]
outdir=home+project_name+'/output/inverse_models/models/'+run_name+'.'+str(num).rjust(4,'0')+'.inv'
#if z_source>high_stress_depth:
# stress=stress_parameter*stress_multiplier
#else:
# stress=stress_parameter
# Frankel 95 scaling of corner frequency #verified this looks the same in GP
# Right now this applies the same factor to all faults
fc_scale=(M0)/(N*stress*dl**3*1e21) #Frankel scaling
small_event_M0 = stress*dl**3*1e21
#Get rho, alpha, beta at subfault depth
zs=fault[kfault,3]
mu,alpha,beta=get_mu(structure,zs,return_speeds=True)
rho=mu/beta**2
#Get radiation scale factor
Spartition=1/2**0.5
if component=='N' :
component_angle=0
elif component=='E':
component_angle=90
rho=rho/1000 #to g/cm**3
beta=(beta/1000)*1e5 #to cm/s
alpha=(alpha/1000)*1e5
#Verified this produces same value as in GP
CS=(2*Spartition)/(4*pi*(rho)*(beta**3))
CP=2/(4*pi*(rho)*(alpha**3))
#Open structure file
mod=loadtxt(home+project_name+'/structure/'+model_name,ndmin=2)
#Get slip quantities
iss=2*arange(len(sol)//2)
ids=2*arange(len(sol)//2)+1
ss=sol[iss]
ds=sol[ids]
#Get total slip
slip=(ss**2+ds**2)**0.5
#Get rigidities and areas
mu=zeros(len(ds))
A=zeros(len(ds))
i=0
for krupt in range((len(sol)//len(f))//2):
for k in range(len(f)):
mu[i]=get_mu(mod,f[k,3])
A[i]=f[k,8]*f[k,9]
i+=1
#Compute moments
try:
M0=mu*A*slip[:,0]
except:
M0=mu*A*slip
#Total up and copute magnitude
M0=M0.sum()
Mw=(2./3)*(log10(M0)-9.1)
return M0,Mw
def get_mean_slip(target_Mw,fault_array,vel_mod):
'''
Depending on the target magnitude calculate the necessary uniform slip
on the fault given the 1D layered Earth velocity model
'''
from numpy import genfromtxt,zeros,ones
from mudpy.forward import get_mu
vel=genfromtxt(vel_mod)
areas=fault_array[:,8]*fault_array[:,9]
mu=zeros(len(fault_array))
for k in range(len(mu)):
mu[k]=get_mu(vel,fault_array[k,3])
target_moment=10**(1.5*target_Mw+9.1)
mean_slip=ones(len(fault_array))*target_moment/sum(areas*mu)
return mean_slip,mu
#if z_source>high_stress_depth:
# stress=stress_parameter*stress_multiplier
#else:
# stress=stress_parameter
# Frankel 95 scaling of corner frequency #verified this looks the same in GP
# Right now this applies the same factor to all faults
fc_scale=(M0)/(N*stress*dl**3*1e21) #Frankel scaling
small_event_M0 = stress*dl**3*1e21
#Get rho, alpha, beta at subfault depth
zs=fault[kfault,3]
mu,alpha,beta=get_mu(structure,zs,return_speeds=True)
rho=mu/beta**2
#Get radiation scale factor
Spartition=1/2**0.5
if component=='N' :
component_angle=0
elif component=='E':
component_angle=90
rho=rho/1000 #to g/cm**3
beta=(beta/1000)*1e5 #to cm/s
alpha=(alpha/1000)*1e5
#Verified this produces same value as in GP
CS=(2*Spartition)/(4*pi*(rho)*(beta**3))
CP=2/(4*pi*(rho)*(alpha**3))
def get_mean_slip(target_Mw,fault_array,vel_mod):
'''
Depending on the target magnitude calculate the necessary uniform slip
on the fault given the 1D layered Earth velocity model
'''
from numpy import genfromtxt,zeros,ones
from mudpy.forward import get_mu
vel=genfromtxt(vel_mod)
areas=fault_array[:,8]*fault_array[:,9]
mu=zeros(len(fault_array))
for k in range(len(mu)):
mu[k]=get_mu(vel,fault_array[k,3])
target_moment=10**(1.5*target_Mw+9.1)
mean_slip=ones(len(fault_array))*target_moment/sum(areas*mu)
return mean_slip,mu