Earthquake as a triggering factor of an infinite slope mechanism

import os
import sys
import subprocess

if 'google.colab' in str(get_ipython()):
    print('Running on CoLab. Installing the required modules...')
    # subprocess.run('pip install ipympl', shell=True);
    subprocess.run('pip install pynewmarkdisp', shell=True);
    from google.colab import output, files
    output.enable_custom_widget_manager()
else:
    import tkinter as tk
    from tkinter.filedialog import askopenfilename

import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
from pynewmarkdisp.newmark import direct_newmark, plot_newmark_integration
from pynewmarkdisp.infslope import factor_of_safety, get_ky

# %matplotlib widget
%matplotlib inline
mpl.rcParams.update({
    "font.family": "serif",
    "font.serif": ["Computer Modern Roman", "cmr", "cmr10", "DejaVu Serif"],  # or 
    "mathtext.fontset": "cm",  # Use Computer Modern fonts for math
    "axes.formatter.use_mathtext": True,  # Use mathtext for axis labels
    "axes.unicode_minus": False,   # Use standard minus sign instead of a unicode character
})

# Create a folder where the records will be saved
workdir = os.getcwd()
records_dir = os.path.join(workdir, "records")
os.makedirs(f"{records_dir}", exist_ok=True)

def load_file(testing_data):
    if testing_data:
        url = "https://raw.githubusercontent.com/eamontoyaa/data4testing/main/pynewmarkdisp/"
        # record = np.loadtxt(f"{url}earthquake_data_friuli.csv", delimiter=";", skiprows=1)
        # record[:, 1] = - 2.0 * 9.81 * record[:, 1]  # Convert to SI units
        record = np.loadtxt(f"{url}earthquake_data_armenia1999.csv", delimiter=",")
    elif testing_data is False and 'google.colab' in str(get_ipython()):
        file = files.upload()
        # Move the file to the records folder
        os.rename(list(file.keys())[0], os.path.join(records_dir, list(file.keys())[0]))
        record = np.loadtxt(os.path.join(records_dir, list(file.keys())[0]))#, delimiter=",")
    else:  # GUI for file selection from local machine if not in CoLab
        tk.Tk().withdraw() # part of the import if you are not using other tkinter functions
        file = askopenfilename()
        # Move the file to the records folder
        os.rename(file, os.path.join(records_dir, os.path.basename(file)))
        record = np.loadtxt(os.path.join(records_dir, os.path.basename(file)))
    time, accel = record.T
    return time, accel
    
def get_arias_intensity(time, accel, g):
    return np.trapz(y=accel**2, x=time) * np.pi / (2 * g)

def get_pgv(time, accel, g):
    accel = 9.81 * accel / g  # Earthquake acceleration to SI units
    vel = np.cumsum(accel) * (time[1] - time[0])
    return np.max(np.abs(vel))

def get_pga(accel, g):
    return np.max(np.abs(accel / g))

testing_data = True  # Set to False to use the GUI to load the data from an external file

# Geometry
β = 35  # [°] - Slope angle - Plane inclination
d = 3.0  # [m] - Depth of the slip surface - Block height
d_w = 3.0  # [m] - Depth of the water table

# Material parameters
φ = 27  # [°] - Friction angle
c = 15  # [kPa] - Cohesion
γ = 19  # [kN/m³] - Unit weight

time, accel = load_file(testing_data)
time, accel

# Value of gravity acceleration in the same units as in the earthquake record
# g = 1.0   # → It means that acceleration units in the input file are given in [g]
g = 9.81  #  → It means that acceleration units in the input file are given in [m/s²]
# g = 981  #  → It means that acceleration units in the input file are given in [cm/s²]

PGA = get_pga(accel, g)
print(f"PGA: {PGA:.2f} g")

fs = factor_of_safety(d, d_w, β, φ, c, γ, ks=0)
print(f"Static factor of safety: {fs:.2f}")

ky = get_ky(d, d_w, β, φ, c, γ)
print(f"Critical seismic coefficient, ky: {ky:.2f}")

ks = 0.7 * ky
fs_ks = factor_of_safety(d, d_w, β, φ, c, γ, ks)
print(f"Seismic coefficient, ks: {ks:.3f}")
print(f"Pseudostatic factor of safety: {fs_ks:.2f}")

decrease_fs = (fs - fs_ks) / fs * 100
print(f"Decrease in factor of safety: {decrease_fs:.2f} %")

ks = 0.5 * PGA
fs_ks = factor_of_safety(d, d_w, β, φ, c, γ, ks)
print(f"Seismic coefficient, ks: {ks:.3f}")
print(f"Pseudostatic factor of safety: {fs_ks:.2f}")

decrease_fs = (fs - fs_ks) / fs * 100
print(f"Decrease in factor of safety: {decrease_fs:.2f} %")

permanent_disp = direct_newmark(time, accel, ky, g)
# permanent_disp = direct_newmark(time, accel, .135, g)
fig = plot_newmark_integration(permanent_disp)
fig.canvas.header_visible = False 
fig.canvas.toolbar_position = 'bottom' 
plt.show()