Earthquake/analysis.m

function param=analysis(param)
%% By L_DelOff
% собираю данные по землетрясению
%% Поехали
type_sensor=param.type_sensor;
DATA.DATA=[];
%% Объединяю данные в диапазоне
if 1
    for i=1:length(param.dates)
        year=param.dates(i).year;
        month=param.dates(i).month;
        day=param.dates(i).day;
        DATA_temp=load_file(param.sattelite_range, year, month, day,...
            type_sensor); %загрузка данных
        DATA.DATA=[DATA.DATA DATA_temp.DATA];
    end
end
%% Сортировка данных
R_quake=param.R_quake;    % радиус действия относительно магнитуды,км (еще тыщу ставил)
% расположение эпицентра
fiA=param.fiA;             % широта
LA=param.LA;             % долгота

fname_rep=param.filename_report;

if exist(fname_rep,'file')
    load(fname_rep,'report');              % тут будут все данные для графика
else
    %report=[];
    %report.data=[];
    report.data.N_sat=[];
    report.data.R=[]; 
end


for i=1:length(param.dates)
    HHeden=[];
    O2den=[];
    Eden=[];
    Etemp=[];
    
    for j=1:length(DATA.DATA)
        date_true=1;
        if DATA.DATA(j).year~=param.dates(i).year % сортировка по году
            date_true=0;
        end
        if DATA.DATA(j).month~=param.dates(i).month % по месяцу
            date_true=0;
        end
        if (DATA.DATA(j).day~=param.dates(i).day)% по числам
            date_true=0;
        end
        if date_true
            for k=1:3 % координаты спутника только 3 раза в минуту(кадр) записываются
                fiB=DATA.DATA(j).GeoLat(k);        % 1. Географическая широта (градусы, север)(Geographic latitude (degrees, north))
                LB=DATA.DATA(j).GeoLong(k);        % 2. Географическая долгота (градусы, восток)(Geographic longitude (degrees, east))
                R=mydistance(fiA,LA,fiB,LB);
                if R<=R_quake
                    HHeden(end+1)= mymean(DATA.DATA(j).RPA_HHeden((k-1)*5+1:(k-1)*5+5),1e10);
                    O2den(end+1) = mymean(DATA.DATA(j).RPA_O2den((k-1)*5+1:(k-1)*5+5),1e10);
                    Eden(end+1)  = mymean(DATA.DATA(j).EP_Eden((k-1)*5+1:(k-1)*5+5),1e10);
                    Etemp(end+1) = mymean(DATA.DATA(j).EP_Etemp((k-1)*5+1:(k-1)*5+5),1e10);
                end
            end
        end
    end
    num=give_number(report,param.sattelite_range,R_quake);
    report.data(num).N_sat=param.sattelite_range;
    report.data(num).R=R_quake;
    report.data(num).date{i}=[myformat(param.dates(i).day,2),'.',...
        myformat(param.dates(i).month,2),'.',...
        myformat(param.dates(i).year,4)];
    report.data(num).RPA_HHeden(i)=mymean(HHeden,1e10);
    report.data(num).RPA_O2den(i)=mymean(O2den,1e10);
    report.data(num).EP_Eden(i)=mymean(Eden,1e10);
    report.data(num).EP_Etemp(i)=mymean(Etemp,1e10);
end

report.data(num).RPA_HHeden_lim=mylim(report.data(num).RPA_HHeden,param);
report.data(num).RPA_O2den_lim=mylim(report.data(num).RPA_O2den,param);
report.data(num).EP_Eden_lim=mylim(report.data(num).EP_Eden,param);
report.data(num).EP_Etemp_lim=mylim(report.data(num).EP_Etemp,param);


param.report=report;
save(fname_rep,'report');
end

function out=mylim(y,param)
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);

out=[];
for i=1:length(y)
    if (y(i)<=L1)||(y(i)>=L2)&&(isnan(y(i)))        
        out(i)=(y(i)-Q2)*100/(L1-Q2);% отклонение в процентах
        %out(i)=1;
    else
        out(i)=0;
    end
end
end


function num=give_number(report,f,R)
flag=0; % Проверяет, есть ли в отчете данные со спутника f
[~,b]=size(report.data);

for i=1:b
    if report.data(i).N_sat==f
        if report.data(i).R==R
            flag=1;
            num=i;
        end
    end
end
if flag==0
    num=b+1;
end
end

function DATA=load_file(sattelite, year, month, day, type_sensor)
success=0;
foldername=['DATA/f' num2str(sattelite) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2)];
target=ls(foldername);
[str row]=size(target);
DATA=[];
for i=1:str
    filename=target(i,:);%end-3:end);
    counter=0;
    for j=1:length(filename)
        if filename(j)==' '
            counter=counter+1;
        end
    end
    filename=filename(1:end-counter);
    if length(filename)>7
        keyword=[num2str_new(year,4) num2str_new(month,2) num2str_new(day,2)];
        flag=0;
        for j=1:length(filename)-length(keyword)
            if filename(j:j+length(keyword)-1)==keyword
                flag=1;
            end
            
        end
        if (prod(filename(end-3:end)=='.mat')==1)&&(flag==1)
            disp(['Найден подходящий файл: ',filename])
            DATA = load([foldername,'/',filename], 'DATA'); % открыть для чтения       
            success=1;            
        end        
    end
end
if ~success     
    disp('Нужных файлов не найдено, перехожу к загрузке')
    loaddata_edr(sattelite, year, month, type_sensor);
    DATA=load_file(sattelite, year, month, day, type_sensor);
end
end