Earthquake/Program.m

152 lines
4.2 KiB
Matlab
Raw Blame History

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

%% 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 = 1;
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