In [1]:
import json
import matplotlib.pyplot as plt
import numpy as np
import pysra
%matplotlib inline
In [2]:
# Increased figure sizes
plt.rcParams['figure.dpi'] = 120
Load the WNA profile from Campbell (2003).
In [3]:
with open('../tests/data/qwl_tests.json') as fp:
data = json.load(fp)[1]
thickness = np.diff(data['site']['depth'])
profile = pysra.site.Profile()
for i, (thick, vel_shear, density) in enumerate(zip(
thickness, data['site']['velocity'], data['site']['density']
)):
profile.append(
pysra.site.Layer(
pysra.site.SoilType(f'{i}', density * pysra.motion.GRAVITY),
thick * 1000, vel_shear * 1000,
)
)
profile.update_layers(0)
Create simple point source motion
In [4]:
motion = pysra.motion.SourceTheoryRvtMotion(
magnitude=6.5, distance=20, region='cena'
)
motion.calc_fourier_amps(data['freqs'])
In [5]:
calc = pysra.propagation.QuarterWaveLenCalculator(site_atten=0.04)
input_loc = profile.location('outcrop', index=-1)
In [6]:
calc(motion, profile, input_loc)
In [7]:
fig, ax = plt.subplots()
ax.plot(motion.freqs, calc.crustal_amp, label='Crustal Amp.')
ax.plot(motion.freqs, calc.site_term, label='Site Term')
ax.set(
xlabel='Frequency (Hz)', xscale='log',
ylabel='Amplitude', yscale='linear',
)
ax.legend()
fig.tight_layout();
The quarter-wavelength calculation is tested against the WNA and CENA crustal amplification models provided by Campbell (2003). The test of the CENA model passes, but the WNA model fails. Below is a comparison of the two crustal amplifications.
In [8]:
fig, ax = plt.subplots()
ax.plot(motion.freqs, calc.crustal_amp, label='Calculated')
ax.plot(data['freqs'], data['crustal_amp'], label='Campbell (03)')
ax.set(
xlabel='Frequency (Hz)', xscale='log',
ylabel='Amplitude', yscale='linear',
)
ax.legend()
fig.tight_layout();
Adjust the profile to match the target crustal amplification -- no consideration of the site attenuation paramater although this can also be done. First, the adjustment is only performed on the velocity. Second set of plots adjusts velocity and thickness.
In [9]:
for adjust_thickness in [False, True]:
calc.fit(target_type='crustal_amp', target=data['crustal_amp'],
adjust_thickness=adjust_thickness)
fig, ax = plt.subplots()
ax.plot(motion.freqs, calc.crustal_amp, label='Calculated')
ax.plot(data['freqs'], data['crustal_amp'], label='Campbell (03)')
ax.set(
xlabel='Frequency (Hz)', xscale='log',
ylabel='Amplitude', yscale='linear',
)
ax.legend()
fig.tight_layout();
for yscale in ['log', 'linear']:
fig, ax = plt.subplots()
ax.plot(profile.initial_shear_vel, profile.depth,
label='Initial', drawstyle='steps-pre')
ax.plot(calc.profile.initial_shear_vel, calc.profile.depth,
label='Fit', drawstyle='steps-pre')
ax.legend()
ax.set(
xlabel='$V_s$ (m/s)', xlim=(0, 3500),
ylabel='Depth (m)', ylim=(8000, 0.1), yscale=yscale,
)
fig.tight_layout();
In [ ]: