====== Задачи ======

===== Задача 1 =====

У данной задачи есть 2 варианта (1.1, 1.2).

Ниже переставлен код выполняющий расчет скорости поперечных и продольных волн в породе состоящий из двух компонентов - минеральной (Кальцит) и пор наполненных водой. Расчет выполняется с использованием методов теории эффективных сред в зависимости от объемных долей компонентов, т.е. от пористости породы.

В данном коде представлено два метода: метод Фойгта-Ройсса-Хилла - ''Voigt_Reuss_Hill_method'' и метод Хашина-Штрикмана для двух компонентов - ''Hashin_Shtrikman_method_2comp''. Оба метода основаны на усреднении некоторых граничных значений. В первом случае границы Фойгта - ''Voigt_bound'' и  Ройсса - ''Reuss_bound'', во втором ''Hashin_Shtrikman_bound_2comp''.

**Задача**
  * задача 1.1 для метода Фойгта-Ройсса-Хилла
  * задача 1.2 для метода Хашина-Штрикмана

Построить на одном рисунке результаты работы **одного** из методов (два или четыре графика на одном figure):
  * Модуль упругости сжатия K, ГПа
  * Модуль упругости сдвига Mu, ГПа
  * Скорость продольных волн Vp, км/с
  * Скорость поперечных волн Vs, км/с

в виде графиков и записать полученные результаты в CSV файл.
На графиках должен быть заголовок, подписи по осям, координатная сетка.

Код выполняющий задание свойств компонентов среды и выполняющий построение графиков и запись результатов должен быть размещен в **отдельном файле** от методов расчета. Код не относящийся к используемому методу должен быть **удален**.

** Базовый код **

<sxh python>
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)    
</sxh>


===== Задача 2 =====

У данной задачи есть 4 варианта (2.1, 2.2, 2.3, 2.4).

==== 2.1 Сейсмические события на удалении от станции ====

**Задача**

Зададимся координатами любого города. Предположим, что в этом городе находится сейсмостанция.
Отметить на карте сам город и выбрать из каталога сейсмических событий события, отстоящие от него по координатам на не менее чем 15 градусов и не более чем на 60 градусов.
Карту следует построить в стереографической проекции с центром в выбранном городе.

**Данные**

{{ :7pluseventsusgs.csv | Каталог сейсмических событий с магнитудой больше 7 USGS}}

** Базовый код **
<sxh python>
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')
</sxh>


==== 2.2 Количество сейсмических событий на сетке ====

Построить на карте мира теплокарту числа сейсмических событий из каталога с шагом 3 на 3 градуса.
Использовать проекцию Робинсона с центром по линии перемены дат.

**Данные**

{{ :7pluseventsusgs.csv | Каталог сейсмических событий с магнитудой больше 7 USGS}}

** Базовый код **
<sxh python>
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')
</sxh>

==== 2.3 Карта напряженности магнитного поля ====

Построить глобальную карту напряженности магнитного поля за 1980, 1990 и 2000 годы.
Использовать проекцию Робинсона. Пример карты [[https://ngdc.noaa.gov/geomag/magfield-wist/]] (Total Intensity: Main field)
Использовать контуры с подписями.

**Данные**

Данные NGDC NOAA {{ ::f_map_mf_1980.zip | 1980 год }} {{ :f_map_mf_1990.zip | 1980 год }} {{ :f_map_mf_2000.zip | 1980 год }}

** Базовый код **
<sxh python>
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')
</sxh>

==== 2.4 Построение смещений на карте ====

В данном Shape-файле описаны точки трех разных видов (kind): base: опорные точки, profile1 и profile2 - результаты измерений смещений вдоль профилей.
Построить карту используя подложку tile сервера OpenStreetMap на которую нанести: опорные точки, точки профилей 1 и 2 соединенные линиями различных цветов.
Для каждой из точек профиля указать направление (az в градусах) и величину смещения (value) в виде стрелки.

**Данные**

{{ ::tdata.zip | Синтетические данные}}

** Базовый код **
<sxh python>
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)
</sxh>