function format_geotail_mag_1(flag,version)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function format_ace_mag_0(flag)
% purpose: convert CDF file format
% Author: James Weygand
% Creation: Feb 27, 2003
% ex: ac_h0_mfi_20000221_v04.cdf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% givens
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
format long e; % use the long format
months = strvcat('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec');
daysofmonths = [31 28 31 30 31 30 31 31 30 31 30 31]; % number of days per months, normal year
daysofmonths2 = [31 29 31 30 31 30 31 31 30 31 30 31]; % number of days per months, leap year.
re = 6378.0; % Earth's radius kilometers
ct =1 ; % set counter to zero
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load in data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yeardoing = double(int32(flag/100)); % get year we are working on
monthdoing = flag - yeardoing*100; % get month we are working on.
if ( (yeardoing==1976) | (yeardoing==1980) | (yeardoing==1984) | (yeardoing==1988) | (yeardoing==1992) | (yeardoing==1996) | (yeardoing==2000) ),
monthend = daysofmonths2(monthdoing); % find the lasst day of the month.
else % must know if this is a leap yr.
monthend = daysofmonths(monthdoing);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creat arrays (extra long)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
full = 60*60*24*31/3.0; % 60 sec by 60 min by 24 hr by 31 days div by resolution to get max size
Mtime = zeros(1,full); % make arrays now because it saves time later.
bxgse = zeros(1,full); % make maximum size estimate
bygse = zeros(1,full);
bzgse = zeros(1,full);
xgse = zeros(1,full);
ygse = zeros(1,full);
zgse = zeros(1,full);
for ii=1:1:monthend
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine file name and open input file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yyyy=num2str(flag); % change yyyymmm into a string
version = 1; % set version to 1
fid = -1; % set file id ito -1
dd = num2str(ii,'%2.2d'); % convert the day into 2 digit form.
while( (version<=6) & (fid <=-1) ), % The pt here is to step
ver = num2str(version); % through all know possible version
ssone = ['ge_k0_mgf_',yyyy,dd,'_v0',ver,'.cdf']; % of the geotail data and load the
outputtxt = [' Input file: ', ssone]; % data. Version numbers vary
disp(outputtxt); % significantly
[fid,message] = fopen(ssone,'r'); % No data check
if fid == -1
disp('Missing Data file! Try again!');
disp(message);
version = version + 1;
end
end
if (fid >-1), % got data well load it and
[alldata,info] = cdfread(ssone); % sort it into the arrays I
len = size(alldata); % understand
for jj = 1:1:len(1)
Mtime(ct) = todatenum(alldata{jj,1});
bxgse(ct) = double(alldata{jj,5}(1))/10.0; % gse
bygse(ct) = double(alldata{jj,5}(2))/10.0; % gse
bzgse(ct) = double(alldata{jj,5}(3))/10.0; % gse
xgse(ct) = double(alldata{jj,8}(1))/re; % km
ygse(ct) = double(alldata{jj,8}(2))/re; % km
zgse(ct) = double(alldata{jj,8}(3))/re; % km
ct = ct+1;
end
clear allmagdata; % Free RAM
end
end
newlll = ct-1; % length of the array with real data
Mtime = Mtime(1:ct-1); % one tomany data points.
bxgse = bxgse(1:ct-1);
bygse = bygse(1:ct-1);
bzgse = bzgse(1:ct-1);
xgse = xgse(1:ct-1);
ygse = ygse(1:ct-1);
zgse = zgse(1:ct-1);

disp('Substituting for error flag with 1.0e+34.');
theirflag=-3200; % the error flag in the geotail CDF files.
flag=1.0e+34; % error flag bob wants
index=find((xgse<=theirflag) | (ygse<=theirflag) | (zgse<=theirflag) | (bxgse<=theirflag/10.0) | (bygse<=theirflag/10.0) | (bzgse<=theirflag/10.0) );
%% eliminate the flags for these. Normally there are very few missing.
%% Removing the flag ones helps prevent corruption of the data on the
%% edges of the flag zones. A 1e34 flag "pulls up" the data on the ends
%% of the flag zones and makes that data worthless. By removing a few
%% flags now it helps to save 3 to 4 data points on each side of the
%% flag. The edges are "pulled up" by the degap program. This program
%% removes bad times, but not bad data.
xgse(index) = [];
ygse(index) = [];
zgse(index) = [];
bxgse(index) = [];
bygse(index) = [];
bzgse(index) = [];
Mtime(index) = [];
newlll = length(Mtime);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check what month we have
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[yrs, monthdoing,dummy, dummy,dummy,dummy] = datevec(Mtime(1)); % Convert time into form I understand
outputtxt = ['Working on month: ', num2str(monthdoing)]; % Check again: what month are we really doing
disp(outputtxt);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine output file name
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ssoutfour = ['geotailmag',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.dat']; % Contruct file names
ssoutfive = ['geotailmag',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.mat']; % Contruct file names
outputtxt = ['Output plasma file: ', ssoutfour];
disp(outputtxt);
% open output file
[fidmagout,message] = fopen(ssoutfour,'w'); % open output file
if fidmagout == -1 % error check
disp(message);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get start and end of the month.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Linearly interpolating to value on the minute (Resampling). ');
startyyymmddhhmmss = [num2str(1,'%2.2d'),'-',months(monthdoing,:),'-',num2str(yrs(1)),' ','00:00:00.00'];
startend(1) = datenum(startyyymmddhhmmss); % Determine start time of month and make string
lyear = yrs(1); % Check if we have a leap year.
if ( (lyear==1976) | (lyear==1980) | (lyear==1984) | (lyear==1988) | (lyear==1992) | (lyear==1996) | (lyear==2000) ),
endyyymmddhhmmss = [num2str(daysofmonths2(monthdoing)),'-',months(monthdoing,:),'-',num2str(yrs(1)),' ','23:59:00.00'];
else
endyyymmddhhmmss = [num2str(daysofmonths(monthdoing)),'-',months(monthdoing,:),'-',num2str(yrs(1)),' ','23:59:00.00'];
end
startend(2) = datenum(endyyymmddhhmmss); % Determine end time of month and make string
% degap requires a specific input
% turn row into column
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert Mag and position from GSE to GSM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Convert Magnetic GSE and position GSE to magnetic field and position GSM');
%RF = [1 2 3; 4 5 6; 7 8 9];
bxgsm = zeros(1,newlll); % these statements save computation time.
bygsm = zeros(1,newlll);
bzgsm = zeros(1,newlll);
xgsm = zeros(1,newlll);
ygsm = zeros(1,newlll);
zgsm = zeros(1,newlll);
for jj=1:1:newlll % convert SC about = GSE to GSM coords
[RF,RI] = gse2gsm(Mtime(jj)); % Builds the matrix
vector = [bxgse(jj),bygse(jj),bzgse(jj)]; % Make a vector
vecout = mat9vec3(RF,vector); % Transfor from gse to gsm
bxgsm(jj) = vecout(1); % Sort our vector components
bygsm(jj) = vecout(2); % Sort our vector components
bzgsm(jj) = vecout(3); % Sort our vector components
vectortwo = [xgse(jj),ygse(jj),zgse(jj)]; % Repeat but for position
vectwoout = mat9vec3(RF,vectortwo);
xgsm(jj) = vectwoout(1);
ygsm(jj) = vectwoout(2);
zgsm(jj) = vectwoout(3);
end
disp('Convertion from GSE to GSM done');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fix the flags
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
flag=1.0e+34;

index=find(bxgsm>=1.0e+20); % When we did GSE to GSM the flag got
bxgsm(index)=repmat(flag,1,length(index)); % all messed up. Fix it.
index=find(bygsm>=1.0e+20);
bygsm(index)=repmat(flag,1,length(index));
index=find(bzgsm>=1.0e+20);
bzgsm(index)=repmat(flag,1,length(index));

index=find(xgsm>=1.0e+20);
xgsm(index)=repmat(flag,1,length(index));
index=find(ygsm>=1.0e+20);
ygsm(index)=repmat(flag,1,length(index));
index=find(zgsm>=1.0e+20);
zgsm(index)=repmat(flag,1,length(index));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% degap requires a specific input construct input array
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nMtime = Mtime'; % rotate arrays. Degape routine wants
nbxgsm = bxgsm'; % a specific format
nbygsm = bygsm';
nbzgsm = bzgsm';
nxgsm = xgsm';
nygsm = ygsm';
nzgsm = zgsm';
inmagarray = [nMtime,nbxgsm,nbygsm,nbzgsm,nxgsm,nygsm,nzgsm]; % make a big matrix
%clear bygsm bzgsm xgsm ygsm zgsm;
outmagarray=degapint(inmagarray,flag,datenum('00:01:00'),startend,datenum('00:10:00')); % Degap it
[maxmagval,lllmag] = max(outmagarray(:,1)); % get output size
nMtime = outmagarray(:,1); % sort back into arrays I
nbxgsm = outmagarray(:,2); % understand
nbygsm = outmagarray(:,3);
nbzgsm = outmagarray(:,4);
nxgsm = outmagarray(:,5);
nygsm = outmagarray(:,6);
nzgsm = outmagarray(:,7);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot out to check
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure; % make fig window
plot(Mtime,bxgsm,'b',nMtime,nbxgsm,'k*'); % plot original and interpolated
datetick('x','dd-mmm-yyyy'); % People friendly date form
if (min(bxgsm) == max(bxgsm)),
ylim([0,1]); % set limits to avoid plotting flags
else
ylim([-50,50]); % set limits to avoid plotting flags
end
clear Mtime vxgsm bxgsm; % free memory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read out to an ascii file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ctout = 1; % set flag
disp('Fixing rounding error in Time.');
newyear = zeros(1,lllmag);
newmon = zeros(1,lllmag);
newday = zeros(1,lllmag);
newhr = zeros(1,lllmag);
newmin = zeros(1,lllmag);
newsec = zeros(1,lllmag);
while( (ctout) <= lllmag),
% fix the rounding error problem.
[newyear(ctout),newmonth(ctout),newday(ctout),newhr(ctout),newmin(ctout),newsec(ctout)]=datevec(outmagarray(ctout,1));
if ( abs(newsec(ctout)-60.0) < 1.0), % here due to rounding error.
newsec(ctout) = 0.000;
end % end of while loop for (ctout) <= lll
if (ctout > 1),
if (newmin(ctout) == newmin(ctout-1)),
newmin(ctout) = newmin(ctout) + 1;
if (newmin(ctout) == 60),
newmin(ctout) = 0;
end
end
end
ctout = ctout + 1; % incriment the counter
atxt = ['Plasma 1: length: ',num2str(lllmag),' ctout: ',num2str(ctout)]; % text out to check
disp(atxt);
end

disp('Making output array');
% contruct output array
outputmagarray = [newday',newmonth',newyear',newhr',newmin',newsec',nbxgsm,nbygsm,nbzgsm,nxgsm,nygsm,nzgsm];
% export output array to file.
ctout = 1; % set flag
while( (ctout) <= lllmag),
count = fprintf(fidmagout,'%2.2d %2.2d %4.4d %2.2d %2.2d %06.3f %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n',outputmagarray(ctout,:));
atxt = ['Mag 2: length: ',num2str(lllmag),' ctout: ',num2str(ctout)];
ctout = ctout + 1; % incriment the counter
end
disp('Closing output file.');
status = fclose (fidmagout);

% read out to matlab binary file.
disp('Making matlab binary file.');
save(ssoutfive,'nMtime','nbxgsm','nbygsm','nbzgsm','nxgsm','nygsm','nzgsm');
% recall we clear the input file varies above
% clear Mtime bxgsm bygsm bzgsm xgsm ygsm zgsm;
% now clear the interpolated data so we can start the next month.
% Recall we will never get the next month.
%clear newyear newday newmonth newhr newmin newsec nbxgsm nbygsm nbzgsm nxgsm nygsm nzgsm;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of Program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%end