У данной задачи есть 2 варианта (1.1, 1.2).
Ниже переставлен код выполняющий расчет скорости поперечных и продольных волн в породе состоящий из двух компонентов - минеральной (Кальцит) и пор наполненных водой. Расчет выполняется с использованием методов теории эффективных сред в зависимости от объемных долей компонентов, т.е. от пористости породы.
В данном коде представлено два метода: метод Фойгта-Ройсса-Хилла - Voigt_Reuss_Hill_method
и метод Хашина-Штрикмана для двух компонентов - Hashin_Shtrikman_method_2comp
. Оба метода основаны на усреднении некоторых граничных значений. В первом случае границы Фойгта - Voigt_bound
и Ройсса - Reuss_bound
, во втором Hashin_Shtrikman_bound_2comp
.
Задача
Построить на одном рисунке результаты работы одного из методов (два или четыре графика на одном figure):
в виде графиков и записать полученные результаты в CSV файл. На графиках должен быть заголовок, подписи по осям, координатная сетка.
Код выполняющий задание свойств компонентов среды и выполняющий построение графиков и запись результатов должен быть размещен в отдельном файле от методов расчета. Код не относящийся к используемому методу должен быть удален.
Базовый код
import numpy as np
import matplotlib.pyplot as plt
#%% Расчет
# Число компонентов
N_components = 2
# Число шагов
N_steps = 100
# Инициализируем массивы исходных данных и результатов
var_c = ['Vp', 'Vs', 'Rho', 'K', 'Mu', 'Vol']
c = np.zeros([N_components,N_steps], dtype=[(x, np.float64) for x in var_c])
var_r = ['K_V', 'Mu_V', 'Vp_V', 'Vs_V',
'K_R', 'Mu_R', 'Vp_R', 'Vs_R',
'K_VRH', 'Mu_VRH', 'Vp_VRH', 'Vs_VRH',
'K_HS_UP', 'Mu_HS_UP', 'K_HS_LO', 'Mu_HS_LO',
'K_HS', 'Mu_HS','Vp_HS','Vs_HS']
r = np.zeros([N_steps], dtype=[(x, np.float64) for x in var_r])
# Задаем исходные данные
# Calcite
c[0]['Vp'] = 6.54
c[0]['Vs'] = 3.35
c[0]['Rho'] = 2.7
# Water
c[1]['Vp'] = 1.6
c[1]['Vs'] = 1e-6
c[1]['Rho'] = 1.1
# Задаем пористость
poros = np.linspace(0.0, 1.0, num=N_steps)
# Определяем содержание компонентов
c[0]['Vol'] = 1 - poros
c[1]['Vol'] = poros
# Расчет
# Расчет K и Mu
c['K'] = c['Rho']*(c['Vp']**2-4./3.*c['Vs']**2)
c['Mu'] = c['Rho']*c['Vs']**2
# axis = 0 : Суммирование по компонентам N_components
# результат - одномерный массив из N_steps элементов
Rho_eff = np.sum(c['Vol']*c['Rho'], axis=0)
# Метод Фойгта
def Voigt_bound(c):
K = np.sum(c['Vol']*c['K'], axis=0)
Mu = np.sum(c['Vol']*c['Mu'], axis=0)
return K, Mu
r['K_V'], r['Mu_V'] = Voigt_bound(c)
# Метод Ройсса
def Reuss_bound(c):
K = 1.0 / np.sum(c['Vol']/c['K'], axis=0)
Mu = 1.0 / np.sum(c['Vol']/c['Mu'], axis=0)
return K, Mu
r['K_R'], r['Mu_R'] = Reuss_bound(c)
# Метод Фойгта-Ройсса-Хилла
def Voigt_Reuss_Hill_method(c):
K_V, Mu_V = Voigt_bound(c)
K_R, Mu_R = Reuss_bound(c)
K = (K_V+K_R) / 2.0
Mu = (Mu_V+Mu_R) / 2.0
return K, Mu
r['K_VRH'] , r['Mu_VRH'] = Voigt_Reuss_Hill_method(c)
# Метод Хашина-Штрикмана для двух компонентов
def Hashin_Shtrikman_bound_2comp(c0,c1):
K1, Mu1, f1 = c0['K'], c0['Mu'], c0['Vol']
K2, Mu2, f2 = c1['K'], c1['Mu'], c1['Vol']
Z0 = K1+4./3.*Mu1
YK = 1./(K2-K1)
YMu = 1./(Mu2-Mu1)
K_HS = K1+f2/(YK+f1/Z0)
Z1 = 2.*f1*(K1+2.*Mu1)/(5.*Mu1*Z0)
Mu_HS = Mu1+f2/(YMu+Z1)
return K_HS, Mu_HS
r['K_HS_UP'], r['Mu_HS_UP'] = Hashin_Shtrikman_bound_2comp(c[0],c[1])
r['K_HS_LO'], r['Mu_HS_LO'] = Hashin_Shtrikman_bound_2comp(c[1],c[0])
def Hashin_Shtrikman_method_2comp(c0,c1):
K_HS_UP, Mu_HS_UP = Hashin_Shtrikman_bound_2comp(c0,c1)
K_HS_LO, Mu_HS_LO = Hashin_Shtrikman_bound_2comp(c1,c0)
K = (K_HS_UP+K_HS_LO) / 2.0
Mu = (Mu_HS_UP+Mu_HS_LO) / 2.0
return K, Mu
r['K_HS'], r['Mu_HS'] = Hashin_Shtrikman_method_2comp(c[0],c[1])
# Для каждого метода рассчитываем Vp и Vs
for m in ['V', 'R', 'VRH', 'HS']:
r[f'Vp_{m}'] = np.sqrt((r[f'K_{m}']+4.0/3.0*r[f'Mu_{m}'])/Rho_eff)
r[f'Vs_{m}'] = np.sqrt(r[f'Mu_{m}']/Rho_eff)
У данной задачи есть 4 варианта (2.1, 2.2, 2.3, 2.4).
Задача
Зададимся координатами любого города. Предположим, что в этом городе находится сейсмостанция. Отметить на карте сам город и выбрать из каталога сейсмических событий события, отстоящие от него по координатам на не менее чем 15 градусов и не более чем на 60 градусов. Карту следует построить в стереографической проекции с центром в выбранном городе.
Данные
Каталог сейсмических событий с магнитудой больше 7 USGS
Базовый код
import numpy as np
import matplotlib.pyplot as plt
from cartopy import crs as ccrs
fig = plt.figure(figsize=(10,10))
central_point = [60,160] # сейсмостанция
ax = fig.add_subplot(projection=ccrs.Stereographic(
central_point[0],central_point[1],scale_factor=3))
ax.set_global()
ax.coastlines(resolution='110m')
Построить на карте мира теплокарту числа сейсмических событий из каталога с шагом 3 на 3 градуса. Использовать проекцию Робинсона с центром по линии перемены дат.
Данные
Каталог сейсмических событий с магнитудой больше 7 USGS
Базовый код
import numpy as np
import matplotlib.pyplot as plt
from cartopy import crs as ccrs
import matplotlib.ticker as mticker
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(projection=ccrs.Robinson(central_longitude=180))
ax.set_global()
ax.coastlines(resolution='110m')
Построить глобальную карту напряженности магнитного поля за 1980, 1990 и 2000 годы. Использовать проекцию Робинсона. Пример карты https://ngdc.noaa.gov/geomag/magfield-wist/ (Total Intensity: Main field) Использовать контуры с подписями.
Данные
Данные NGDC NOAA 1980 год 1980 год 1980 год
Базовый код
import numpy as np
import matplotlib.pyplot as plt
from cartopy import crs as ccrs
import matplotlib.ticker as mticker
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(projection=ccrs.Robinson(central_longitude=0))
ax.set_global()
ax.coastlines(resolution='110m')
В данном Shape-файле описаны точки трех разных видов (kind): base: опорные точки, profile1 и profile2 - результаты измерений смещений вдоль профилей. Построить карту используя подложку tile сервера OpenStreetMap на которую нанести: опорные точки, точки профилей 1 и 2 соединенные линиями различных цветов. Для каждой из точек профиля указать направление (az в градусах) и величину смещения (value) в виде стрелки.
Данные
Базовый код
import numpy as np
import matplotlib.pyplot as plt
from cartopy import crs as ccrs
import shapefile
from cartopy.io.img_tiles import OSM
shpfile = shapefile.Reader("/home/f0ma/tdata/tdata")
fields = [x[0] for x in shpfile.fields[1:] ]
for item in shpfile.iterShapeRecords():
print("Position:", item.shape.points[0])
print("Data:", dict(zip(fields, item.record)))
tiles_provider = OSM()
plt.figure(figsize=(16, 16))
ax = plt.axes(projection=tiles_provider.crs)
ax.set_extent([shpfile.bbox[0]-0.1, shpfile.bbox[2]+0.1,
shpfile.bbox[1]-0.1, shpfile.bbox[3]+0.1], ccrs.PlateCarree())
ax.add_image(tiles_provider, 11)