How to use mudpy - 10 common examples

top='free' ; bottom='locked' ; left='locked' ; right='locked' #'locked' or 'free'
bounds=(top,bottom,left,right)
################################################################################e=
########      Run-time modifications to the time series             ############
weight=True
decimate=None  #Decimate by constant (=None for NO decimation)
bandpass=None#np.array([0.5]) #Corner frequencies in Hz =None if no filter is desired
################################################################################

#Initalize project folders
if init==1:
    runslip.init(home,project_name)
    
# Run green functions          
if make_green==1 or make_synthetics==1:
    runslip.inversionGFs(home,project_name,GF_list,tgf_file,fault_name,model_name,
        dt,tsun_dt,NFFT,tsunNFFT,make_green,make_synthetics,dk,pmin,
        pmax,kmax,beta,time_epi,hot_start)   

#Run inversion
if invert==1:
    runslip.run_inversion(home,project_name,run_name,fault_name,model_name,GF_list,G_from_file,
            G_name,epicenter,rupture_speed,num_windows,reg_spatial,reg_temporal,
            nfaults,beta,decimate,bandpass,solver,bounds,weight,Ltype)

rupture_speed=3.0  #Only necessary if onset times are not identified in rupt file
stf_type='dreger'
################################################################################


#Initalize project folders
if init==1:
    runslip.init(home,project_name)

# Run green functions          
if make_green==1: 
    if ncpus<2:
        runslip.make_green(home,project_name,station_file,fault_name,model_name,
            dt,NFFT,static,tsunami,hot_start,dk,pmin,pmax,kmax)
    else:
        runslip.make_parallel_green(home,project_name,station_file,fault_name,
            model_name,dt,NFFT,static,tsunami,hot_start,dk,pmin,pmax,kmax,ncpus)


#Now make synthetics for source/station pairs
if make_synthetics==1:
    if ncpus<2: 
        runslip.make_synthetics(home,project_name,station_file,fault_name,model_name,integrate,
                static,tsunami,beta,hot_start,time_epi)
    else:
        runslip.make_parallel_synthetics(home,project_name,station_file,fault_name,model_name,integrate,
                static,tsunami,beta,hot_start,time_epi,ncpus,custom_stf=None)
    
#Run forward comptuation or solve for inverse problem
if solve==1:
    if static==0: #Forward problem (full waveforms)
        forward.waveforms_fakequakes(home,project_name,fault_name,None,GF_list,

# [0.5] is a low pass filter
 # [0.02,0.5] is a bdan pass filter
 # [0.02,np.inf] is a high pass filter
displacement_bandpass=np.array([0.5]) 
velocity_bandpass=None
tsunami_bandpass=None
bandpass=[displacement_bandpass,velocity_bandpass,tsunami_bandpass]
################################################################################

#Initalize project folders
if init==1:
    runslip.init(home,project_name)
    
# Run green functions          
if make_green==1 or make_synthetics==1:
    runslip.inversionGFs(home,project_name,GF_list,tgf_file,fault_name,model_name,
        dt,tsun_dt,NFFT,tsunNFFT,make_green,make_synthetics,dk,pmin,
        pmax,kmax,beta,time_epi,hot_start,ncpus,custom_stf)   

#Run inversion
if invert==1:
    runslip.run_inversion(home,project_name,run_name,fault_name,model_name,GF_list,G_from_file,
            G_name,epicenter,rupture_speed,num_windows,reg_spatial,reg_temporal,
            nfaults,beta,decimate,bandpass,solver,bounds,weight,Ltype)

nstrike=297 ; ndip=1 ; nfaults=(nstrike,ndip) #set nstrike to total no. of faults and ndip to 1 if using Tikh
beta=45 #Rotational offset (in degrees) applied to rake (0 for normal)
Ltype=0 # 0 for Tikhonov and 2 for Laplacian
solver='nnls' # 'lstsq','nnls'
top='free' ; bottom='locked' ; left='locked' ; right='locked' #'locked' or 'free'
bounds=(top,bottom,left,right)
################################################################################e=
########      Run-time modifications to the time series             ############
weight=True
decimate=None  #Decimate by constant (=None for NO decimation)
bandpass=None#np.array([0.5]) #Corner frequencies in Hz =None if no filter is desired
################################################################################

#Initalize project folders
if init==1:
    runslip.init(home,project_name)
    
# Run green functions          
if make_green==1 or make_synthetics==1:
    runslip.inversionGFs(home,project_name,GF_list,tgf_file,fault_name,model_name,
        dt,tsun_dt,NFFT,tsunNFFT,make_green,make_synthetics,dk,pmin,
        pmax,kmax,beta,time_epi,hot_start)   

#Run inversion
if invert==1:
    runslip.run_inversion(home,project_name,run_name,fault_name,model_name,GF_list,G_from_file,
            G_name,epicenter,rupture_speed,num_windows,reg_spatial,reg_temporal,
            nfaults,beta,decimate,bandpass,solver,bounds,weight,Ltype)

epicenter=array([84.708,28.147,15]) 
resample=1 #Resample synthetics to this rate (in Hz)
integrate=1 #=0 produces velocities, =1 makes displacements
beta=0 #Rake offset, usually a good idea to keep at zero
rupture_speed=3.0 #Fastest rupture allowed in km/s
num_windows=1
################################################################################


#Initalize project folders
if init==1:
    runslip.init(home,project_name)

# Run green functions          
if make_green==1: 
    runslip.make_green(home,project_name,station_file,fault_name,model_name,dt,NFFT,static,
                tsunami,hot_start,dk,pmin,pmax,kmax)  

#Now make synthetics for source/station pairs
if make_synthetics==1:
    #runslip.rupt2fault(home,project_name,rupture_name) 
    runslip.make_synthetics(home,project_name,station_file,fault_name,model_name,integrate,
                static,tsunami,beta,hot_start,time_epi)
    
#Run forward comptuation or solve for inverse problem
if solve==1:
    if static==0: #Forward problem (full waveforms)
        forward.waveforms_matrix(home,project_name,fault_name,rupture_name,station_file,GF_list,model_name,
                run_name,epicenter,time_epi,integrate,tsunami,hot_start,resample,beta,rupture_speed,
                num_windows,dt,NFFT)
    if static==1: #Forward problem (coseismics)
        forward.coseismics(home,project_name,rupture_name,station_file)

runslip.init(home,project_name)

# Run green functions          
if make_green==1: 
    if ncpus<2:
        runslip.make_green(home,project_name,station_file,fault_name,model_name,
            dt,NFFT,static,tsunami,hot_start,dk,pmin,pmax,kmax)
    else:
        runslip.make_parallel_green(home,project_name,station_file,fault_name,
            model_name,dt,NFFT,static,tsunami,hot_start,dk,pmin,pmax,kmax,ncpus)


#Now make synthetics for source/station pairs
if make_synthetics==1:
    if ncpus<2: 
        runslip.make_synthetics(home,project_name,station_file,fault_name,model_name,integrate,
                static,tsunami,beta,hot_start,time_epi)
    else:
        runslip.make_parallel_synthetics(home,project_name,station_file,fault_name,model_name,integrate,
                static,tsunami,beta,hot_start,time_epi,ncpus,custom_stf=None)
    
#Run forward comptuation or solve for inverse problem
if solve==1:
    if static==0: #Forward problem (full waveforms)
        forward.waveforms_fakequakes(home,project_name,fault_name,None,GF_list,
                model_name,run_name,dt,NFFT,G_from_file,G_name,source_time_function=stf_type,
                stf_falloff_rate=4.0,rupture_name=rupture_name,epicenter=epicenter,time_epi=time_epi)
    if static==1: #Forward problem (coseismics)
        forward.coseismics_matrix(home,project_name,rupture_name,station_file,G_from_file,G_name)

integrate=0 #=0 produces velocities, =1 makes displacements
beta=0 #Rake offset, usually a good idea to keep at zero
num_windows=1
rupture_speed=3.0  #Only necessary if onset times are not identified in rupt file
stf_type='dreger'
################################################################################


#Initalize project folders
if init==1:
    runslip.init(home,project_name)

# Run green functions          
if make_green==1: 
    if ncpus<2:
        runslip.make_green(home,project_name,station_file,fault_name,model_name,
            dt,NFFT,static,tsunami,hot_start,dk,pmin,pmax,kmax)
    else:
        runslip.make_parallel_green(home,project_name,station_file,fault_name,
            model_name,dt,NFFT,static,tsunami,hot_start,dk,pmin,pmax,kmax,ncpus)


#Now make synthetics for source/station pairs
if make_synthetics==1:
    if ncpus<2: 
        runslip.make_synthetics(home,project_name,station_file,fault_name,model_name,integrate,
                static,tsunami,beta,hot_start,time_epi)
    else:
        runslip.make_parallel_synthetics(home,project_name,station_file,fault_name,model_name,integrate,
                static,tsunami,beta,hot_start,time_epi,ncpus,custom_stf=None)
    
#Run forward comptuation or solve for inverse problem

################################################################################

############                 Synthetics parameters               ###############
time_epi=UTCDateTime('2015-04-25T06:11:26')
epicenter=array([84.708,28.147,15]) 
resample=1 #Resample synthetics to this rate (in Hz)
integrate=1 #=0 produces velocities, =1 makes displacements
beta=0 #Rake offset, usually a good idea to keep at zero
rupture_speed=3.0 #Fastest rupture allowed in km/s
num_windows=1
################################################################################


#Initalize project folders
if init==1:
    runslip.init(home,project_name)

# Run green functions          
if make_green==1: 
    runslip.make_green(home,project_name,station_file,fault_name,model_name,dt,NFFT,static,
                tsunami,hot_start,dk,pmin,pmax,kmax)  

#Now make synthetics for source/station pairs
if make_synthetics==1:
    #runslip.rupt2fault(home,project_name,rupture_name) 
    runslip.make_synthetics(home,project_name,station_file,fault_name,model_name,integrate,
                static,tsunami,beta,hot_start,time_epi)
    
#Run forward comptuation or solve for inverse problem
if solve==1:
    if static==0: #Forward problem (full waveforms)
        forward.waveforms_matrix(home,project_name,fault_name,rupture_name,station_file,GF_list,model_name,

################################################################################


#Initalize project folders
if init==1:
    runslip.init(home,project_name)

# Run green functions          
if make_green==1: 
    runslip.make_green(home,project_name,station_file,fault_name,model_name,dt,NFFT,static,
                tsunami,hot_start,dk,pmin,pmax,kmax)  

#Now make synthetics for source/station pairs
if make_synthetics==1:
    #runslip.rupt2fault(home,project_name,rupture_name) 
    runslip.make_synthetics(home,project_name,station_file,fault_name,model_name,integrate,
                static,tsunami,beta,hot_start,time_epi)
    
#Run forward comptuation or solve for inverse problem
if solve==1:
    if static==0: #Forward problem (full waveforms)
        forward.waveforms_matrix(home,project_name,fault_name,rupture_name,station_file,GF_list,model_name,
                run_name,epicenter,time_epi,integrate,tsunami,hot_start,resample,beta,rupture_speed,
                num_windows,dt,NFFT)
    if static==1: #Forward problem (coseismics)
        forward.coseismics(home,project_name,rupture_name,station_file)

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