# -*- coding: utf-8 -*-
"""
Created on Mon Nov 5 15:21:14 2018
@author: Alejandro
"""
#%%
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
#datos = np.loadtxt("3/09141_step3171.dat",skiprows=9)
#datos = np.loadtxt("3/09141_step3187.dat",skiprows=9)
#datos = np.loadtxt("3/09141_step3229.dat",skiprows=9)
#datos = np.loadtxt("3/09141_step3247.dat",skiprows=9)
#datos = np.loadtxt("3/09141_step3264.dat",skiprows=9)
datos = np.loadtxt("537/00537_step2104.dat",skiprows=9)
datos2 = np.loadtxt("537/moment005372.dat",skiprows=9)
grid = 2 ; ds = 4*4**(grid)
dt = ds/2/6000
x = datos[:,0]; y = datos[:,1];
slip_vel = datos[:,2]; # m/s
slip = datos[:,3] # mm
shear_stres = datos[:,4] # MPa
t = datos2[:,0] # s
moment_rate = datos2[:,1] # Nm/s
moment = datos2[:,2] # Nm
Mw = datos2[:,3]
#%%
#% Scatter 3D
function = slip*4**(grid)
fig = plt.figure(figsize=(10,6))
ax = fig.gca(projection='3d')
#tt = list(range(-1,8))
sct = ax.scatter(x*ds, y*ds, function,vmin=0, vmax=np.max(function), s = 50, alpha
= 1 , c=function, cmap='jet', linewidth=0.1, marker = "s"); #velocidad
deslizamientos
#sct = ax.scatter(x*ds, y*ds, function,vmin=0, vmax=4, s = 40, alpha = 1 ,
c=function, cmap='jet', linewidth=0.1, marker = "s"); #deslizamientos
#sct = ax.scatter(x*4*4**(grid), y*4*4**(grid), function, s = 50, c=function,
cmap='jet', linewidth=0.1);
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_zticklabels([]); #ax.set_zlabel('MPa')
#ax.set_zlim(-1, 7);
ax.set_xlim(np.min(x*4*4**(grid)), np.max(x*4*4**(grid)));
ax.set_ylim(np.min(y*4*4**(grid)), np.max(y*4*4**(grid)))
fig.colorbar(sct, shrink=0.5, aspect=5, label='Deslizamiento [mm]');#
fig.set_label('label')
ax.view_init(elev=90., azim=0) # Vista en planta
#ax.view_init(elev=20., azim=30) # Vista 3D
plt.title('Deslizamiento terremoto Mw = 3.77943, t = 0.554667 [s]')
plt.show()
#%%
"""
LEER TODOS LOS ARCHIVOS PRODUCIDOS POR HIDEO AOCHI
"""
import glob
import numpy as np
import matplotlib.pyplot as plt
import obspy as ob
filenames = sorted(glob.glob('537/00537_step2*.dat'))
# 0: x, 1: y, 2: slip vel [m/s], 3: slip [mm], 4: shear stress [MPa]
col = 2
grid = 2 ; ds = 4*4**(grid)
ds = 4*4**(grid)
dt = ds/2/6000
fun_0 = np.loadtxt('537/00537_step0001.dat')
X = fun_0[:,0]
Y = fun_0[:,1]
fun = fun_0[:,col]
fun2 = fun_0[:,4]
fun3 = fun_0[:,3]
for f in filenames:
data = np.loadtxt(fname=f)
fun = np.c_[fun, data[:,col]]
fun2 = np.c_[fun2, data[:,4]]
fun3 = np.c_[fun3, data[:,3]]
t2 = dt*np.linspace(0,len(fun[0,:]),len(fun[0,:]))
#%%
Datos_finales = np.loadtxt("537/output00537f.dat",skiprows=9)
time_step = Datos_finales[:,0]
time = Datos_finales[:,1] # [s]
moment2 = Datos_finales[:,2] # [Nm/s]
moment_rate2 = Datos_finales[:,3] # [Nm]
#%%
from scipy.stats.stats import pearsonr
ss = np.where(fun == np.max(fun[:,:]))
int = np.int(ss[0])
plt.figure()
plt.subplot(221)
plt.plot(t2,fun[int,:],'k',linewidth=1.5)
plt.xlabel('Tiempo [s]'); plt.ylabel('Velocidad de deslizamiento [m/s]')
plt.xlim(t2[0], t2[len(t2)-1])
#plt.figure()
plt.subplot(222)
plt.plot(t2,fun2[int,:],'k',linewidth=1.5)
plt.xlabel('Tiempo [s]'); plt.ylabel('Esfuerzo de Corte [MPa]')
plt.xlim(t2[0], t2[len(t2)-1])
#plt.figure()
plt.subplot(223)
plt.plot(t2,fun3[int,:]*4**(grid),'k',linewidth=1.5)
plt.xlabel('Tiempo [s]'); plt.ylabel('Deslizamientos [mm]')
plt.xlim(t2[0], t2[len(t2)-1])
plt.subplot(224)
plt.plot(time, moment_rate2,'k',linewidth=1.5)
plt.xlabel('Tiempo [s]'); plt.ylabel('Tasa de momento sísmico [Nm/s] \n de toda la
ruptura')
plt.xlim(time[0], time[len(time)-1])
plt.suptitle('Parámetros de ruptura \n simulación ruptura dinámica terremoto Mw =
3.77943 ',size=16)
p = pearsonr(fun[int,:],fun2[int,:]) # Correlation and p-value
print('La correlación es', p[0])
print('El punto está en', (x[int])*4*4**(grid),(y[int])*4*4**(grid))
#%% realizar filtro pasabajo 5 Hz
#%% calcular Gc (energía de fractura)
## fun2[int,:] # Esfuerzo de corte
## fun3[int,:] # Deslizamientos
##import matplotlib.pyplot as plt
##
##fun2[int,0:65] = fun2[int,0:65]
##fun2[int,65:] = 0
plt.figure()
plt.plot((fun3[int,:])*(10**(-3)),fun2[int,:],'o-k')
plt.xlabel('Deslizamientos [m]'); plt.ylabel('Esfuerzo de corete [MPa]')
#plt.plot(fun3[500,:],fun2[500,:],'k')
index = np.where( fun2[int,:] == 0 ) # indices para encontrar el valor cuando la
tensión cae a cero.
index = np.asmatrix(index)
dc1 = (fun3[int,index[0,0]])*(10**(-3)) # Distancia a la cual se alcanza por
primera vez el esfuerzo cero.
tau1 = np.max(fun2[int,:])
G_c1 = (tau1*dc1)/2
print('La energía de fractura para ese punto es:', G_c1,'M J/m/m')
#%% Energía de Fractura para toda la falla
n = len(fun2[:,0])
dc = (np.zeros((n)))
tau = (np.zeros((n)))
G_c = (np.zeros((n)))
for k in range(n):
if (( max(fun3[k,:]) > 0 )):
if (min(fun2[k,:]) == 0 ):
index = np.where( fun2[k,:] == 0 ) # indices para encontrar el valor
cuando la tensión cae a cero.
index = np.asmatrix(index)
dc[k] = (fun3[k,index[0,0]])*(10**(-3)) # en metros
tau[k] = (np.max(fun2[k,:])) # MPa
G_c[k] = (((tau[k])*(dc[k]))/2 )*(10**(6)) # J/m/m
# else:
# dc[k] = 0
# tau[k] = 0
# G_c[k] = (0.16)
else:
dc[k] = 0
tau[k] = 0
G_c[k] = 3* (1.25*(10**(3))) # J/m/m
# pass
#%% Graficación de energía de fractura
fig = plt.figure(figsize=(10,6))
ax = fig.gca(projection='3d')
sct = ax.scatter(X*ds, Y*ds, G_c, s = 50, alpha = 1 , c=G_c, cmap='jet',
linewidth=0.1, marker = "s"); #velocidad deslizamientos
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_zticklabels([]); #ax.set_zlabel('MPa')
#ax.set_zlim(-1, 7);
ax.set_xlim(np.min(X*4*4**(grid)), np.max(X*4*4**(grid)));
ax.set_ylim(np.min(Y*4*4**(grid)), np.max(Y*4*4**(grid)))
fig.colorbar(sct, shrink=0.5, aspect=5, label='Energía de Fractura [J/m**2]');#
fig.set_label('label')
ax.view_init(elev=90., azim=0) # Vista en planta
#ax.view_init(elev=20., azim=30) # Vista 3D
plt.title('Deslizamiento terremoto Mw = 3.77943, t = 0.554667 [s]')
plt.show()
#%% Relación de caída de esfuerzo propemdio y ponderado
# slip : fun3
# shear stres : fun2
I_1 = np.zeros((len(fun3)))
I_2 = np.zeros((len(fun3)))
I_3 = np.zeros((len(fun3)))
#dsigma = np.zeros((len(fun3)))
du = np.zeros((len(fun3)))
for i in range(len(I_1)):
dsigma = max(fun2[i,:]) - fun2[i,len(fun2[0,:])-1]
du[i] = fun3[i,len(fun3[0,:])-1]
I_1[i] = dsigma*(du[i])*ds
I_2[i] = (du[i])*ds
I_3[i] = dsigma*ds
ddddd = np.where(du > 0)
A = (len(ddddd[0]))*(4*4**(grid))
dsigmaE = ((np.sum(I_1)) / (np.sum(I_2)))
dsigmaA = np.sum(I_3) /A
% Programa de correlaciones TwoPoint Statistic
clear all; clc;% close all;
datos = load('1534/01534_step3066.dat');
Z = 3 ;
grid = 4*(4^(Z)) ;
x = (datos(:,1))*grid ; % metros
y = (datos(:,2))*grid ; % metros
slip_vel = datos(:,3) ; % m/s
slip = datos(:,4)*grid/4 ; % mm
shear_stress = datos(:,5) ; % MPa
H = vec2mat(slip,64);
%
X = 1:64; X = X*grid;
Y = 1:64; Y = Y*grid;
[X,Y] = meshgrid(X,Y);
figure()
s = surf(X,Y,H);colormap jet; axis xy;
s.EdgeColor = 'none';
c = colorbar;
c.Label.String = 'Esfuerzo de corte [MPa]';
xlim( [x(1) , x(end)] ); ylim( [y(1) , y(end)] );
title('Esfuerzo de corte final Terremoto Mw4.62181')
xlabel({'Distancia [m]'});ylabel({'Distancia [m]'})
view(2)
%% Two point statistic
%[ corrfun r rw] = twopointcorr( x,y,dr,blksize,verbose)
[kx , ky ] = find( 0 < H & H <= 2000 ) ;
%k = (0 < H) & (H < 200);
%ky = (0< H) & (H < 220);
%dr = grid*2
dr =550
[ corrfun, r, rw] = twopointcorr( kx*grid,ky*grid,dr,1000);
%[ corrfun, r, rw] = tw opointcorr( x(k),y(k),dr,1000);
%[ corrfun, r, rw] = twopointcorr( kx,ky,dr);
%figure(); plot(r, corrfun/max(corrfun), '-k'); xlabel('Distancia [m]');
ylabel('Funci�n de correlaci�n')
figure()
hold on
plot(r, corrfun); xlabel('Distancia [m]'); ylabel('Funci�n de correlaci�n')
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 4 15:53:58 2018
@author: Alejandro
"""
#%% Datos de caida de esfuerzos
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
#from mpl_toolkits.mplot3d import Axes3D
datos = np.loadtxt("datos.txt")
numero = datos[:,0]
Mw = datos[:,1]
S_A = datos[:,2]
S_E = datos[:,3]
function = Mw
#plt.scatter(S_A, S_E, function,vmin=np.max(function), vmax=np.max(function), s =
50, alpha = 1 , c=function, cmap='jet', linewidth=0.1, marker = "s"); #velocidad
deslizamientos
sc = plt.scatter(S_E, S_A,s = 100 ,vmin=np.min(function), vmax=np.max(function),
alpha = 1 , c=function, cmap='jet', linewidth=0.1);
plt.title('Caída de esfuerzos')
plt.xlabel('$\Delta \sigma_E$ [MPa]'); plt.ylabel('$\Delta \sigma_A$ [MPa]')
#plt.xlim(0, n);plt.ylim(0, n)
plt.colorbar(sc,label='Magnitud de momento Mw')
plt.show()
