diff --git a/loaddata_edr_old.m b/loaddata_edr_old.m new file mode 100644 index 0000000..2012f89 --- /dev/null +++ b/loaddata_edr_old.m @@ -0,0 +1,410 @@ +function loaddata_edr_old %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 i10 + 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 \ No newline at end of file