lab 9 Earthquakes - FouadM0426/comp.physics-spring-2025 GitHub Wiki

Goal

In this lab we will be analyzing the earthquake of Mexico City that erupted in 1985.We will be using Fourier Transform to try to figure out why buildings that were in the 6-15 story range were the majority of buildings destroyed in a crash, Aswell as finding the ideal height for the building to be able to withstand said earthquake if it would erupt again. The Fourier will be used to calculate the dominant frequencies, the soil near the area amplifies the waves.

Part 1

In this section of the lab, we will be taking all data from earthquakes from 1984-1986 that have been over a magnitude of 2.5 on the scale. For reference a 4 on the scale is what humans are able to physically sense. The code below will show importing the data and showing a visual to have a better understanding.

import numpy as np
import matplotlib.pyplot as plt
!wget https://facultyweb.siena.edu/~rfinn/teaching/phys250/earthquakes-1984-1986.csv
from astropy.table import Table
equakes = Table.read('earthquakes-1984-1986.csv')
plt.scatter(equakes['longitude'], equakes['latitude'], s=equakes['mag']*15)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Earthquakes from 1984-1986')

image

  • This portion will show how the magnitude is distributed, using a histogram
# plot histogram of magnitudes
mybins = np.arange(2,9,0.5)
plt.hist(equakes['mag'], bins=mybins)
plt.yscale('log')
plt.xlabel('Magnitude')
plt.ylabel('Number of Earthquakes')
plt.title('Distribution of Earthquake Magnitudes from 1984-1986')
fraction = len(equakes[equakes['mag'] == 8]) / len(equakes[equakes['mag'] > 2.5])
print(fraction)
print(f'(The mexico city magnitude of 8 is {fraction*100}that is the percent)')

image

Part 2

In this section of the lab we will be analyzing the seismic waves.The data we will be importing is broken down into 5 columns which will be used to plot different factors on a graph such as time vs displacement.

!wget https://facultyweb.siena.edu/~rfinn/teaching/phys250/mex1985_waveform.csv
# your code to read in the file use Table.read
seismic_data = Table.read('mex1985_waveform.csv')

# print the table so you can see the names of the columns
print(seismic_data)
plt.figure(figsize=(10, 10))
plt.subplot(3, 1, 1)
plt.plot(seismic_data['time'], seismic_data['accel'])
plt.xlabel('Time (s)')
plt.ylabel('Acceleration (cm/s^2)')
plt.title('Acceleration vs Time')
plt.subplot(3, 1, 2)
plt.plot(seismic_data['time'], seismic_data['vel'])
plt.xlabel('Time (s)')
plt.ylabel('Velocity (cm/s)')
plt.title('Velocity vs Time')
plt.subplot(3, 1, 3)
plt.plot(seismic_data['time'], seismic_data['displ'])
plt.xlabel('Time (s)')
plt.ylabel('Displacement (cm)')
plt.title('Displacement vs Time')
plt.tight_layout()
plt.show()

image

  • Determining the frequency spectrum for the waves
Y_accel = np.fft.rfft(seismic_data['accel'])
Y_vel = np.fft.rfft(seismic_data['vel'])
Y_displ = np.fft.rfft(seismic_data['displ'])
frequencies = np.fft.rfftfreq(len(seismic_data['accel']), d=seismic_data['time'][1] - seismic_data['time'][0])
plt.figure(figsize=(10, 10))
plt.subplot(3, 1, 1)
plt.plot(frequencies, np.abs(Y_accel))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.plot(frequencies, np.abs(Y_vel))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.plot(frequencies, np.abs(Y_displ))

image

plt.xlim(0, .75)
plt.plot(frequencies, np.abs(Y_accel))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')

image

Part 3 changes to construction

In this portion we will be using f0=10/Nstories, this being the relationship between resonant frequency/ period and the number of stories to help us find the ideal number of stories to prevent structure collapse

freq = (0.5)
N_stories = 10 / freq
print(N_stories)

This will resolve and show us that the number of stories we need are 20 for the optimal amount of stories to prevent collapse for high magnitude earthquakes.