Files
Earthquake/Program.m

152 lines
4.2 KiB
Matlab
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

%% By L_DelOff
% скрипт обеспечивает:
% 1. Закачку данных
% 2. Построение графиков
% 3. Дополнительный анализ
%% Начало
close all
clear all
% Дано
param.k=1.2; % коэффициент для расчёта доверительного интервала
%% Подгрузка сведений об катаклизме
% N = 1 - Япония [11.03.2011]
% N = 2 - Чили [16.09.2015]
% N = 3 - Перу [26.05.2019]
% N = 4 - Охотское море [24.05.2013]
% N = 5 - Залив Аляска [23.01.2018]
% N = 6 - Бейрут [04.08.2020]
N = 6;
prev=10;
next=5;
param.prev=prev;
load katalog
param.M=katalog(N).M; % магнитуда
% расположение эпицентра
param.fiA=katalog(N).fi; % широта
param.LA=katalog(N).L; % долгота
param.sattelite_available=katalog(N).sattelite; % доступные спутники
dd=katalog(N).day;
mm=katalog(N).month;
yyyy=katalog(N).year;
param.dates=day_array(dd,mm,yyyy,prev,next);
param.day_x=dd;
param.filename_report=[katalog(N).Location,'_', myformat(katalog(N).day,2)...
, myformat(katalog(N).month,2), myformat(katalog(N).year,4),'.mat'];
param.type_sensor='ssies';
%% Первичная информация
for i=param.sattelite_available
for j=[10^(0.43*param.M) 1000]
param.sattelite_range=i;%15:18; % Выбор с какого спутника требуются данные
param.R_quake=j; % радиус действия относительно магнитуды,км (еще тыщу ставил)
param=analysis(param);
grafik(param,katalog(N))
end
end
%% Корреляционные матрицы
param=analysis_corr(param);
%% Сохранение в файл
export_report(param.filename_report)
disp('Готово')
disp('--------')
function grafik(param,katalog)
win1=figure;
win1.Name=[katalog.Location,...
', ',myformat(katalog.day,2)...
,'.', myformat(katalog.month,2)...
,'.', myformat(katalog.year,4)...
,' [',num2str(katalog.fi),';',num2str(katalog.L),']',...
' M=',num2str(param.M),...
' R=',num2str(fix(param.R_quake)), 'км Sat: F-',num2str(param.sattelite_range),...
' K=',num2str(param.k)];
win1.Units='normalized';
win1.OuterPosition = [0 0 1 1];
t1=uicontrol(win1,'Style','text');
t1.Units='Normalized';
t1.Position = [0.1 0.98 0.8 0.02];
t1.String = win1.Name;
t1.FontSize = 12;
t1.BackgroundColor=[1 1 1];
XLim=[param.dates(1).day param.dates(end).day];
N_sat=findNsat(param.report,param.sattelite_range,param.R_quake);
ax1=subplot(2,2,1,'Parent',win1);
y=param.report.data(N_sat).RPA_HHeden;
explot(ax1,y,XLim,param)
title('H-He');
ax2=subplot(2,2,2,'Parent',win1);
y=param.report.data(N_sat).EP_Etemp;
explot(ax2,y,XLim,param)
title('E_{temp}');
ax3=subplot(2,2,3,'Parent',win1);
y=param.report.data(N_sat).RPA_O2den;
explot(ax3,y,XLim,param)
title('O_2');
ax4=subplot(2,2,4,'Parent',win1);
y=param.report.data(N_sat).EP_Eden;
explot(ax4,y,XLim,param)
title('E');
saveas(win1,[katalog.Location,...
'_', myformat(katalog.day,2)...
, myformat(katalog.month,2)...
, myformat(katalog.year,4)...
,'_R=', num2str(fix(param.R_quake)), '_F-', num2str(param.sattelite_range)],'bmp');
end
function N_sat=findNsat(report,f,R)
[~,b]=size(report.data);
N_sat=[];
for i=1:b
if report.data(i).N_sat==f
if report.data(i).R==R
N_sat=i;
end
end
end
end
function limit=explot(ax,y,setXLim,param)
%x=[];
%for i=1:length(param.dates)
% x(i)=param.dates(i).day;
%end
plot(ax,y,'Color','k','Marker','o');
xticks(1:length(param.dates))
N_sat=findNsat(param.report,param.sattelite_range,param.R_quake);
xtlabels={param.report.data(N_sat).date};
xticklabels(xtlabels{1})
xtickangle(45)
grid on
hold on
Q1=quantile(y,0.25);
Q2=quantile(y,0.5);
Q3=quantile(y,0.75);
L1=Q2-param.k*(Q3-Q1);
L2=Q2+param.k*(Q3-Q1);
XLim=ax.XLim;
plot(ax,XLim,[L1 L1],'Color','b','LineStyle','--');
plot(ax,XLim,[L2 L2],'Color','b','LineStyle','--');
plot(ax,XLim,[Q2 Q2],'Color','m','LineStyle','--');
YLim=ax.YLim;
plot(ax,[1+param.prev 1+param.prev],YLim,'Color','r');
ax.XLim=[1 length(param.dates)];
end