143 lines
3.9 KiB
Matlab
143 lines
3.9 KiB
Matlab
%% By L_DelOff
|
||
% скрипт обеспечивает:
|
||
% 1. Закачку данных
|
||
% 2. Построение графиков
|
||
% 3. Дополнительный анализ
|
||
|
||
%% Начало
|
||
% Дано
|
||
param.k=1.2; % коэффициент для расчёта доверительного интервала
|
||
%% Подгрузка сведений об катаклизме
|
||
% N = 1 - Япония [11,3,2011]
|
||
% N = 2 - Чили [16,9,2015]
|
||
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);
|
||
%% Превышение порога
|
||
%param=analysis_limit(param);
|
||
|
||
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)];
|
||
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 |