function loaddata_edr %DATA=loaddata_edr(filename)
%% By L_DelOff
% Загрузка данных из файла *.EDR
% Описание файла приведено в:
% !docs/!AFRL ASCII and Binary File Format Descriptions.pdf
%% Поехали
warning off
sattelite_range=15:18; % Выбор с какого спутника требуются данные
year_range=2015; %
month_range=9; %
%% Скачать данные в формате архива
if 1
    for sattelite=sattelite_range %6:20
        for year=year_range % года
            for month=month_range % месяца
                download_fun(sattelite, year, month); %скачивание и обработка нужных данных
            end
        end
    end
end
%% Распаковать архив
if 1
    for sattelite=sattelite_range %6:20
        for year=year_range % года
            for month=month_range % месяца
                extract_data_fun(sattelite, year, month); %скачивание и обработка нужных данных
            end
        end
    end
end
%% Преобраховать EDR в mat
if 1
    for sattelite=sattelite_range %6:20
        for year=year_range % года
            for month=month_range % месяца
                convert_data_fun(sattelite, year, month); %скачивание и обработка нужных данных
            end
        end
    end
end


end

function convert_data_fun(sattelite, year, month)
type_sensor='ssies';
foldername=['DATA/f' num2str(sattelite) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2)];
target=ls(foldername);
[str row]=size(target);
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)>4
        if filename(end-3:end)=='.EDR'            
            message=['Расшифровка файла:',filename,'\n'];
            fprintf(message);
            
            fid = fopen([foldername,'/',filename], 'rt'); % открыть для чтения
            counter=1;
            while ~feof(fid)    % пока не конец файла
                frame=[];       % здесь будет 1 кадр (1 минута) в текстовом формате
                for i=1:114
                    if ~feof(fid)
                        frame(i).str = fgets(fid);  % Последовательно читаем из файла строки
                    end
                end
                
                if counter==1
                    DATA=getframe_edr(frame);          % Считываем кадр
                else
                    DATA(counter)=getframe_edr(frame);          % Считываем кадр
                end
                counter=counter+1;
                if mod(counter,100)==0
                    disp(num2str(counter));
                end
            end
            
            save([foldername,'/',filename(1:end-4),'.mat'],'DATA');
            
            message=['Расшифровка файла закончена\n'];
            fprintf(message);
        end
    end
end
end

function extract_data_fun(sattelite, year, month)
type_sensor='ssies';
foldername=['DATA/f' num2str(sattelite) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2)];
target=ls(foldername);
[str row]=size(target);
for i=1:str
    if target(i,end-2:end)=='.gz'
        message=['Распаковка:',target(i,:),'\n'];
        fprintf(message);
        gunzip([foldername,'/',target(i,:)],foldername);
        message=['Распакован\n'];
        fprintf(message);
    end
end
end

function download_fun(sattelite, year, month)
% https://satdat.ngdc.noaa.gov/dmsp/data/f18/ssies/2020/05/
% https://satdat.ngdc.noaa.gov/dmsp/data/f18/ssies/2020/05/PS.CKGWC_SC.U_DI.A_GP.SIES3-F18-R99990-B9999090-APGA_AR.GLOBAL_DD.20200501_TP.000001-235959_DF.EDR.gz
type_sensor='ssies';
message=['Ищем данные спутника F' num2str(sattelite) ' по дате:' num2str_new(year,4) '/' num2str_new(month,2) '...'];
fprintf(message)
url_site=['https://satdat.ngdc.noaa.gov/dmsp/data/f' num2str_new(sattelite,2) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2) '/'];
fullList = webread(url_site);
message='готово\n';
fprintf(message)
url=extract_url(fullList);

for i=1:length(url)
    message=['число ' num2str(i) '...'];
    fprintf(message);
    url_file = [url_site url(i).file];
    filename = ['DATA/f' num2str(sattelite) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2) '/' url(i).file];
    mkdir(['DATA/f' num2str(sattelite) '/' type_sensor '/' num2str_new(year,4) '/' num2str_new(month,2)])
    outfilename = websave(filename,url_file);
    message='скачано\n';
    fprintf(message);
    
end

a=1;
end

function url=extract_url(text)
counter=0;
i=1;
while i<length(text)
    if text(i)=='"'
        a=i;
        i=i+1;
        while text(i)~='"'
            i=i+1;
        end
        b=i;
        if b-a>10
            c=text(b-6:b-1);
            if c=='EDR.gz'
                counter=counter+1;
                url(counter).file=text(a+1:b-1);
            end
        end
    end
    i=i+1;
end
end



% %% Инициализация
% %load DATA
% %% Открытие файла если уже нет DATA
% if exist('DATA')~=1
%     fid = fopen('DATA/PS.CKGWC_SC.U_DI.A_GP.SIES3-F18-R99990-B9999090-APGA_AR.GLOBAL_DD.20140831_TP.000001-235959_DF.EDR', 'rt'); % открыть для чтения
%     counter=1;
%     while ~feof(fid)    % пока не конец файла
%         frame=[];       % здесь будет 1 кадр (1 минута) в текстовом формате
%         for i=1:114
%             if ~feof(fid)
%                 frame(i).str = fgets(fid);  % Последовательно читаем из файла строки
%             end
%         end
%
%         if counter==1
%             DATA=getframe_edr(frame);          % Считываем кадр
%         else
%             DATA(counter)=getframe_edr(frame);          % Считываем кадр
%         end
%         counter=counter+1;
%         if mod(counter,100)==0
%             disp(num2str(counter));
%         end
%     end
% end
%
% for i=1:length(DATA)
%     x(i)=(DATA(i).GeoLong(1));
%     y(i)=(DATA(i).GeoLat(1));
% end
% plot(x,y)
% xlabel('Geographic longitude (degrees, east)')
% ylabel('Geographic latitude (degrees, north)')
%end

function FRAME=getframe_edr(frame)
HEALTH=0; % "здоровье кадра"
%% строка 1: Пустая строка (Blank line)
if isspace(frame(1).str) HEALTH=HEALTH+1;
end
%% строка 2: Шапка + данные о ПО(если это первый кадр в файле)
if frame(2).str~="" HEALTH=HEALTH+1; end
%% строка 3: 5 значений
z = textscan(frame(3).str,'%f','EmptyValue',-Inf);
FRAME.N_record_bin = z{1}(1);   %1. Номер записи в двоичном файле, откуда был взят этот EDR.(Number of the record in the binary file from where this EDR was taken.)
FRAME.N_EDR_bin = z{1}(2);      %2. Номер EDR в записи в двоичном файле (1-3).(Number of the EDR within the record in the binary file (1-3).)
FRAME.SatID = z{1}(3);          %3. ID спутника.(Satellite Flight ID (two digit integer))
FRAME.year = fix(z{1}(4)/1e4);  %4. Год, месяц, день, час, минута
FRAME.month = mod(fix(z{1}(4)/1e2),1e2);
FRAME.day = mod(z{1}(4),1e2);
FRAME.hour = fix(z{1}(5)/1e2);
FRAME.minute = mod(z{1}(5),1e2);
if length(z{1})==5 HEALTH=HEALTH+1; end
%% строка 4:  ЭФЕМЕРИДА (EPHEMERIS)
if frame(4).str(1:end-1)=="EPHEMERIS" HEALTH=HEALTH+1; end
%% строки 5-7: по 6 значений
for i=0:2
    a = textscan(frame(5+i).str,'%f','EmptyValue',-Inf);
    if length(a{1})==6 HEALTH=HEALTH+1; end
    % данные на секунде 00 20 40 (i=0 1 2)
    FRAME.GeoLat(i+1)=a{1}(1);     % 1. Географическая широта (градусы, север)(Geographic latitude (degrees, north))
    FRAME.GeoLong(i+1)=a{1}(2);    % 2. Географическая долгота (градусы, восток)(Geographic longitude (degrees, east))
    FRAME.ApLat(i+1)=a{1}(3);      % 3. Вершинная широта (градусы, север)(Apex latitude (degrees, north))
    FRAME.ApLong(i+1)=a{1}(4);     % 4. Вершинная долгота (градусы, восток)(Apex longitude (degrees, east))
    FRAME.ApLocalTime(i+1)=a{1}(5);% 5. Местное время Apex (часы)(градусы, восток)(Apex local time (hours))
    FRAME.SatAlt(i+1)=a{1}(6);     % 6. Высота (км)(Satellite altitude (km))
end
%% строка 8:  Потенциал спутника (“SATELLITE POTENTIAL, LAST = SOURCE”)
if frame(8).str(1:end-1)=="SATELLITE POTENTIAL, LAST = SOURCE" HEALTH=HEALTH+1; end
%% строки 9-10:  16 значений (2х8)
a = textscan(frame(9).str,'%f','EmptyValue',-Inf);
b = textscan(frame(10).str,'%f','EmptyValue',-Inf);
if length(a{1})==8 HEALTH=HEALTH+1; end
if length(b{1})==8 HEALTH=HEALTH+1; end
FRAME.SatPot(1:8)=a{1}(1:8);    % значения для потенциала спутника (Vbias + VIP) в вольтах каждые 4 секунды
FRAME.SatPot(9:15)=b{1}(1:7);   % values for the satellite potential (Vbias+VIP) in volts every 4 seconds
FRAME.SatPotSrc=b{1}(8);        % Измеритель потенциала (an integer (1-2) indicating the satellite potential source where:
% 1 - as set by the on-board microprocessor
% 2 - as set by the SENPOT sensor)
%% строка 11:  Потенциал спутника (“PRIMARY PLASMA DENSITY, THEN SOURCE”)
if frame(11).str(1:end-1)=="PRIMARY PLASMA DENSITY, THEN SOURCE" HEALTH=HEALTH+1; end
%% строки 12-21:  60 значений (6х10)
% Односекундные средние значения плотности первичной плазмы (# / см3)
% One-second averages of the primary plasma density (#/cm3)
z=[];
for i=12:21
    a = textscan(frame(i).str,'%f','EmptyValue',-Inf);
    if length(a{1})==6 HEALTH=HEALTH+1; end
    z=[z a{1}.'];
end
FRAME.PrimPlasmaDen=z;
%% строка 22:  Источник данных по плазме
% An integer (1-3) indicating the plasma density source where:
% 1 - Ion density from SM sensor
% 2 - Ion density from DM sensor
% 3 - Electron density from EP sensor (DC Mode)
a = textscan(frame(22).str,'%f','EmptyValue',-Inf);
if length(a{1})==1 HEALTH=HEALTH+1; end
FRAME.PrimPlasmaDenSrc=a{1};
%% строка 23:  Горизонтальный ионный дрейф (“HORIZONTAL ION DRIFT VELOCS”)
if frame(23).str(1:end-1)=="HORIZONTAL ION DRIFT VELOCS" HEALTH=HEALTH+1; end
%% строки 24-33:  60 значений (6х10)
% Односекундные значения скорости горизонтального дрейфа (м/с) для времен от HHMM00 до HHMM59 с 6 значениями в строке и всего 60 значениями.
% One-second values for the horizontal drift speed (m/s) for times HHMM00 through HHMM59 with 6 values per line and 60 values total.
z=[];
for i=24:33
    a = textscan(frame(i).str,'%f','EmptyValue',-Inf);
    if length(a{1})==6 HEALTH=HEALTH+1; end
    z=[z a{1}.'];
end
FRAME.HorIonDrift=z;
%% строка 34:  Вертикальный ионный дрейф (“VERTICAL ION DRIFT VELOCS”)
if frame(34).str(1:end-1)=="VERTICAL ION DRIFT VELOCS" HEALTH=HEALTH+1; end
%% строки 35-44:  60 значений (6х10)
% Односекундные значения скорости вертикального дрейфа (м/с) для времен от HHMM00 до HHMM59 с 6 значениями в строке и всего 60 значениями.
% One-second values for the vertical drift speed (m/s) for times HHMM00 through HHMM59 with 6 values per line and 60 values total.
z=[];
for i=35:44
    a = textscan(frame(i).str,'%f','EmptyValue',-Inf);
    if length(a{1})==6 HEALTH=HEALTH+1; end
    z=[z a{1}.'];
end
FRAME.VerIonDrift=z;
%% строка 45:  (“CKL ANALYSES, THEN SOURCE”)
if frame(45).str(1:end-1)=="CKL ANALYSES, THEN SOURCE" HEALTH=HEALTH+1; end
%% строка 46-63:
for i=1:6 % для секунд: 5 15 25 35 45 55
    z=[];
    a = textscan(frame(46+(3*(i-1))).str,'%f','EmptyValue',-Inf);
    if length(a{1})==4 HEALTH=HEALTH+1; end
    FRAME.CKL_RMS(i)=a{1}(1);   % CKL Analysis: (RMS ΔN)/N (%)
    FRAME.CKL_T1(i)=a{1}(2);    % CKL Analysis: T1
    FRAME.CKL_p1(i)=a{1}(3);    % CKL Analysis: p1
    FRAME.CKL_CKL(i)=a{1}(4);   % CKL Analysis: CKL
    a = textscan(frame(47+(3*(i-1))).str,'%f','EmptyValue',-Inf);
    if length(a{1})==8 HEALTH=HEALTH+1; end
    z=[z a{1}.'];
    a = textscan(frame(48+(3*(i-1))).str,'%f','EmptyValue',-Inf);
    if length(a{1})==8 HEALTH=HEALTH+1; end
    z=[z a{1}.'];
    FRAME.CKL_PDS(i,:)=z(1:end-1);       %CKL Analysis: Decimated power density spectrum (PDS) for time period.
end
%% строка 64:  Источник данных по CLK
% An integer (1-3) indicating data source used for CKL calculation where:
% 1 - SM density data only
% 2 - SM density and filter data
% 3 - EP DC mode density data
a = textscan(frame(64).str,'%f','EmptyValue',-Inf);
if length(a{1})==1 HEALTH=HEALTH+1; end
FRAME.CLK_Src=a{1};
%% строка 112 (потом понадобится)
a = textscan(frame(112).str,'%f','EmptyValue',-Inf);
if length(a{1})==7 HEALTH=HEALTH+1; end
FRAME.ADC_temp=a{1}(2);     %2: ADC temperature (degrees C)
FRAME.SEP_temp=a{1}(3);     %3: SEP temperature
FRAME.DM_offsetV=a{1}(4);   %4: DM offset voltage (volts) if IES-2 or RPA plasma plate potential (volts) if SSIES-3
FRAME.DM_mode=a{1}(5);      %5: DM mode (integer)
% SSIES-2:
% -1: Undefined
% 0: Normal
% 1-8: H+
% 9: FIBA
% SSIES-3:
% 0: Slow
% 1: Normal
FRAME.EP_mode=a{1}(6);      % 6: EP mode (0-6 : A,B,BS,C,D,DS,E) (integer)
FRAME.EP_VIP=a{1}(7);       % 7: VIP at EDR start (volts)
%% строки 65-80
if (FRAME.EP_mode==0)||(FRAME.EP_mode==1)||(FRAME.EP_mode==2)||(FRAME.EP_mode==6)
    % строка 65 (“EP SWEEP ANALYSES SETS”)
    if frame(65).str(1:end-1)=="EP SWEEP ANALYSES SETS" HEALTH=HEALTH+1; end
    % EP Sweep analyses (every 4 seconds). There are 15 EP sweep analysis
    % sets. Each is valid for either 4 (modes A, B and BS) or 2 (mode E)
    % seconds centered on the time specified in the set. Each analysis set
    % contains the following parameters:
    for i=1:15
        a = textscan(frame(66+(i-1)).str,'%f','EmptyValue',-Inf);
        if length(a{1})==6 HEALTH=HEALTH+1; end
        FRAME.EP_SCT(i)=a{1}(1);    % 1: Sweep center time (UT, seconds) (integer)
        FRAME.EP_Eden(i)=a{1}(2);   % 2: Electron density (el/cm3)
        FRAME.EP_Etemp(i)=a{1}(3);  % 3: Electron temperature (degrees K)
        FRAME.EP_SatPot(i)=a{1}(4); % 4: Satellite potential (volts)
        FRAME.EP_Qualif(i)=a{1}(5); % 5: Analysis qualifier (integer). Set to zero if the on-board microprocessor did not perform the analysis, per the flag in element 415. If the on-board microprocessor was in use, then it is set to the MP EP flags word from Word 11 of Cycle 1 in the telemetry. This word can also be zero.
        FRAME.EP_PES(i)=a{1}(6);    % 6: EP photo-electron surrogate value.
    end
end
%% строка 81 (“EP ANALYSES SOURCE”)
if frame(81).str(1:end-1)=="EP ANALYSES SOURCE" HEALTH=HEALTH+1; end
%% строка 82
% An Integer (1-2) indicating EP analysis source where
% 1 – Ground processing analysis
% 2 – On-board microprocessor analysis
a = textscan(frame(82).str,'%f','EmptyValue',-Inf);
if length(a{1})==1 HEALTH=HEALTH+1; end
FRAME.EP_src=a{1};
%% строка 83 (“RPA SWEEP ANALYSES SETS, THEN SOURCE”)
if frame(83).str(1:end-1)=="RPA SWEEP ANALYSES SETS, THEN SOURCE" HEALTH=HEALTH+1; end
%% строка 84-98 (“RPA SWEEP ANALYSES SETS, THEN SOURCE”)
% RPA Sweep analyses (every 4 seconds). There are 15 RPA sweep analysis
% sets. Each is valid for the 4 seconds centered on the time specified in
% the set. Each analysis set contains the following parameters:
for i=1:15
    a = textscan(frame(84+(i-1)).str,'%f','EmptyValue',-Inf);
    if length(a{1})==8 HEALTH=HEALTH+1; end
    FRAME.RPA_SCT(i)=a{1}(1);       % 1: Sweep center time (UT, seconds) (integer)
    FRAME.RPA_O2den(i)=a{1}(2);     % 2: O+ density (ion/cm3)
    FRAME.RPA_HHeden(i)=a{1}(3);    % 3: Total (H+ + He+) density (ion/cm3)
    FRAME.RPA_LIF(i)=a{1}(4);       % 4: Light ion flag (integer)
    % 0 - No light ion
    % 1 - Light ion is H+
    % 2 - Light ion is He+
    % 3± = 3 + 10000 x (H+ fraction)
    FRAME.RPA_IonNemp(i)=a{1}(5);   % 5: Ion temperature (degrees K)
    FRAME.RPA_RamIonDrift(i)=a{1}(6);% 6: Ram ion drift velocity (m/s)
    FRAME.RPA_qualif(i)=a{1}(7);    % 7: Analysis qualifier (integer)
    % 0 - Analysis terminated unsuccessfully
    % 1 - Successful analysis
    FRAME.RPA_totIonDens(i)=a{1}(8);% 8: RPA-derived total ion density
end
%% строка 99
% An integer (1-2) indicating RPA analysis source where:
% 1 - Ground processing analysis
% 2 - On-board microprocessor analysis
a = textscan(frame(99).str,'%f','EmptyValue',-Inf);
if length(a{1})==1 HEALTH=HEALTH+1; end
FRAME.RPA_src=a{1};
%% строка 100 (“DM ION DENSITY”)
if frame(100).str(1:end-1)=="DM ION DENSITY" HEALTH=HEALTH+1; end
%% строки 101-110:  60 значений (6х10)
%DM one-second average ion density (ion/cm3) with 6 values per line and 60
% values total
z=[];
for i=101:110
    a = textscan(frame(i).str,'%f','EmptyValue',-Inf);
    if length(a{1})==6 HEALTH=HEALTH+1; end
    z=[z a{1}.'];
end
FRAME.DM_IonDen=z;
%% строка 111 (“ENGINEERING DATA”)
if frame(111).str(1:end-1)=="ENGINEERING DATA" HEALTH=HEALTH+1; end
%% строка 113 (“FILLER”)
if frame(113).str(1:end-1)=="FILLER" HEALTH=HEALTH+1; end
%% строка 114 (“FILLER”)
a = textscan(frame(114).str,'%f','EmptyValue',-Inf);
if length(a{1})>=1 HEALTH=HEALTH+1; end
FRAME.Filler=a{1};

FRAME.HEALTH=HEALTH;
end