Source Code for 1-min MAGDAS data

This code can reconstruct the H, D, F, Z magnetogram from the .mgd files

function [MagData,Header,STATDATA,CORRECT_INF]=read_1M(FileName)

%============================================================================

%

% [MagData,Header,STATDATA,CORRECT_INF]=read_1M(FileName)

%

% Read 1-Minute data of MAGDAS staraged data

%

% Input: FileName

% Outout: MagData : Magnetic Data [H,D,Z,F,Inc_EW,Inc_NS,Temp] (1440x7/1day)

% Header : Header Information (Structure array)

% STATDATA : Status Data of the Magnetometer (Structure Array)

% CORRECT_INF : GPS Time Correct Information (Structure Array)

%

% 2003/10/09 Ver. 1.0 by K. Kitamura

% 2004/07/22 Ver. 1.1 (時???Z?ｳ回?狽窿wッダ?[の繰り返し０に対応?j

%============================================================================

%clear all

FileName=('ABU_MIN_201102060000.mgd')

%================= Open file ===================

[Fid,Message]=fopen(FileName,'r','l');

if Fid<0,

disp(Message)

return

end

% ============== Find magnetic data ======================

ReadUnit=1000; FindMag=0; NumRead=0;

Status=fseek(Fid,0,'bof');

while FindMag==0,

[Buff,Count]=fread(Fid,ReadUnit,'uchar');

[Message,Errnum]=ferror(Fid);

if Errnum~=0,

error(Message)

end

NumRead=NumRead+1;

Find1A=find((Buff-26)==0);

if ~isempty(Find1A),

for i=1:length(Find1A),

if Find1A(i)~=ReadUnit,

if Buff(Find1A(i)+1)==0,

MagPos=(NumRead-1)*ReadUnit+Find1A(i)+1;

FindMag=1; break;

end

else

[Buff,Count]=fread(Fid,ReadUnit,'uchar');

[Message,Errnum]=ferror(Fid);

if Errnum~=0,

error(Message)

end

NumRead=NumRead+1;

if Buff(1)==0,

MagPos=(NumRead-1)*ReadUnit+1;

FindMag=1; break;

end

end

end

end

end

fclose(Fid);

% ============== Read header =============================

[Fid,Message]=fopen(FileName,'rt');

if Fid<0,

disp(Message)

return

end

Status=fseek(Fid,0,'bof');

head=fgetl(Fid);

Header.DATA_TYPE=sscanf(head(10:end),'%s');

head=fgetl(Fid);

Header.SERIAL_NUM=str2num(head(11:end));

head=fgetl(Fid);

Header.STATION_NAME=sscanf(head(13:end),'%s');

head=fgetl(Fid);

Header.GEODETIC_LATITUDE=str2num(head(18:end));

head=fgetl(Fid);

Header.GEODETIC_LONGITUDE=str2num(head(19:end));

head=fgetl(Fid);

Header.RECORD_TIME=sscanf(head(12:end),'%s');

head=fgetl(Fid);

Header.REPORTED=sscanf(head(9:end),'%s');

head=fgetl(Fid);

Header.SAMPLING_TIME=str2num(head(14:end));

head=fgetl(Fid);

Header.MAG_RANGE=str2num(head(10:end));

head=fgetl(Fid);

Header.HEADER_DATA_NUM=str2num(head(18:end));

if Header.HEADER_DATA_NUM~=0

for i=1:Header.HEADER_DATA_NUM

head=fgetl(Fid);

STATDATA(i).HEADER_DATA_TIME=sscanf(head(17:end),'%s');

head=fgetl(Fid);

STATDATA(i).BAIAS_H=str2num(head(8:end));

head=fgetl(Fid);

STATDATA(i).BAIAS_D=str2num(head(8:end));

head=fgetl(Fid);

STATDATA(i).BAIAS_Z=str2num(head(8:end));

head=fgetl(Fid);

STATDATA(i).BAIAS_F=str2num(head(8:end));

end

else

STATDATA=[];

end

head=fgetl(Fid);

Header.CORRECT_NUM=str2num(head(12:end));

if Header.CORRECT_NUM~=0

for i=1:Header.CORRECT_NUM

head=fgetl(Fid);

CORRECT_tmp=sscanf(head(12:end),'%s');

CORRECT_INF(i).TIME=CORRECT_tmp(1:8);

CORRECT_INF(i).VALUE=str2num(CORRECT_tmp(10:end));

end

else

CORRECT_INF=[];

end

fclose(Fid);

% HEADER{1}=Header;

% HEADER{2}=STATDATA;

% HEADER{3}=CORRECT_INF;

%================== Read Data ====================

[Fid,Message]=fopen(FileName,'r','l');

% ============== Calculate size of magnetic data =========

Status=fseek(Fid,0,'eof');

EOFid=ftell(Fid);

DataSize=EOFid-MagPos;

%========================================

Status=fseek(Fid,MagPos,'bof');

[MagData,Count]=fread(Fid,[7,DataSize/28],'float32');

%========== convert error data (NaN) to 99999.99=================

[errorrow,errorcol]=find(isnan(MagData)==1);

MagData(errorrow,errorcol)=999999.9;

MagData=MagData';

fclose(Fid);

thr1= [0:60:1440];

thr1= [ 60 120 180 240 300 360 420 480 540 600 660 720 780 840 900 960 1020 1080 1140 1200 1260 1320 1380 1440];

subplot(4,1,1)

plot(MagData(:,1))

set(gca,'XTick',[thr1])

set(gca,'XLim',[0 1440])

set (gca,'FontSize',10)

xlabel('Local time (Hrs)','FontSize',10)

ylabel('H (nT)','FontSize',10)

title('ABU\it{H(nT) 06 Feb 2011}','FontSize',14)

%title('\it{Diurnal variation of Day to day variability of H values 061106/07}','FontSize',14)

%

subplot(4,1,2)

plot(MagData(:,2))

set(gca,'XTick',[thr1])

set(gca,'XLim',[0 1440])

set (gca,'FontSize',10)

xlabel('Local time (Hrs)','FontSize',10)

ylabel('D(nT)','FontSize',10)

title('ABU\it{D(nT) 06 Feb 2011}','FontSize',14)

%title('\it{Diurnal variation of Day to day variability of D values 061106/07}','FontSize',14)

%

subplot(4,1,3)

plot(MagData(:,3))

set(gca,'XTick',[thr1])

set(gca,'XLim',[0 1440])

set (gca,'FontSize',10)

xlabel('Local time (Hrs)','FontSize',10)

ylabel('Z(nT)','FontSize',10)

title('ABU\it{Z(nT) 06 Feb 2011}','FontSize',14)

subplot(4,1,4)

plot(MagData(:,4))

set(gca,'XTick',[thr1])

set(gca,'XLim',[0 1440])

set (gca,'FontSize',10)

xlabel('Local time (Hrs)','FontSize',10)

ylabel('Z(nT)','FontSize',10)

title('ABU\it{Z(nT) 06 Feb 2011}','FontSize',14)

dummy=GenerateMin(MagData);

Hmin=dummy(:,1);

Dmin=dummy(:,2);

Zmin=dummy(:,3);

Fmin=dummy(:,4);

Mag_dat=[Hmin;Dmin;Zmin;Fmin];

% [nrow,ncols]= size(Mag_dat);

Mag_data=(reshape(Mag_dat, length(Mag_dat)/4, 4));

% ddata1=[ddata1;(reshape(data1, ncols, length(data1)/ncols))']

Mag_data=[thr1', Mag_data];

save H06Feb2011ABU.mat;

%title('\it{Diurnal variation of Z values 061106/07}','FontSize',14)

%save dummy.dat Hmin Dmin Zmin Fmin -ascii

%----------------------------------------

% %efforts to evaluate hourly means

% %-----------------------------------

function y=GenerateMin(x)

% Reduce data to every minute average

deltat=60;

Hrpday=24;

for i=1:Hrpday,

for j=1:4,

aux=x((i-1).*deltat+1:i.*deltat,j); % Take every 60 seconds data

I=find(not(isnan(aux)));

if length(I)>=10,

y(i,j)=mean(aux(I));

end

end

end