How to use the mudpy.forward.get_mu function in mudpy

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