function format_geotail_plasma_1(flag)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  function format_geotail_plasma_0(flag)
%  purpose: convert Geotail LEP data CDF file format 
%  into ascii and matlab file out puts. Data is
%  interpolated to one minute and gaps filled with flags.
%  Author: James Weygand
%  flag = yyyymm (year and month)
% Creation: April 1, 2003
% example input: ge_k0_lep_19990911_v01.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.
k_b = 1.380622e-23; 	% boltzman constant
m_p =1.672e-27;	    	% mass of a proton.
re = 6378.0;  	   	% kilometers
pie = 3.141592654;	% pi
ct =1;			% pre-set to one.
degtorad = 0.0174532;	% degrees to radians conversion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.
vmag    = zeros(1,full);	% make maximum size estimate
theta   = zeros(1,full);
phi     = zeros(1,full);
xgse    = zeros(1,full);
ygse    = zeros(1,full);
zgse    = zeros(1,full);
density = zeros(1,full);
therm   = 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_lep_',yyyy,dd,'_v0',ver,'.cdf'];  % of the geotail data and load the
        	outputtxt = [' Input file: ', ssone];a		% data. Version numbers vary
        	disp(outputtxt);				% significantly
		[fid,message] = fopen(ssone,'r');		% No data check
        	if fid == -1
                	disp('Missing Data file!');
                	disp(message);
			version = version + 1;
        	end
	end
		if (fid	>-1),
			[alldata,info] = cdfread(ssone);		% got data well load it and
			len = size(alldata);				% sort it into the arrays I
			for jj = 1:1:len(1)				% understand
				Mtime(ct)   = todatenum(alldata{jj,1});
                		density(ct) = alldata{jj,9}; 		%0.01/cc
                		therm(ct)   = 1.0e+34;     		% No thermal values
                		vmag(ct)    = alldata{jj,10}(1);	%km/s
                		theta(ct)   = alldata{jj,10}(2); 	%0.1 deg
                		phi(ct)     = alldata{jj,10}(3); 	%0.1 deg
                		xgse(ct)    = alldata{jj,6}(1);		%km
                		ygse(ct)    = alldata{jj,6}(2);		%km
                		zgse(ct)    = alldata{jj,6}(3);		%km

				ct = ct+1;
			end
			clear alldata;					% Free Ram
	        end
end	
	newlll = ct-1;				% length of the array with real data	
	Mtime   = Mtime(1:ct-1);		% one tomany data points.
	density = density(1:ct-1)/100.0;
	therm   = therm(1:ct-1);
	vmag    = vmag(1:ct-1);
	theta   = (theta(1:ct-1)/10.0)*degtorad;
	phi     = (phi(1:ct-1)/10.0)*degtorad;
	xgse    = xgse(1:ct-1)/re;
	ygse    = ygse(1:ct-1)/re;
	zgse    = zgse(1:ct-1)/re;

disp('Substituting for error flag with 1.0e+34.');
theirflag=-3200;				% the error flag in the geotail CDF files.
flag=1.0e+34;a					% error flag bob wants
index=find((xgse<=theirflag) | (ygse<=theirflag) | (zgse<=theirflag) | (vmag<=theirflag) | (theta<=degtorad*theirflag/10.0) | (phi<=degtorad*theirflag/10.0) | (density<=theirflag/100.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) 	= [];
vmag(index) 	= [];
theta(index) 	= [];
phi(index) 	= [];
density(index)  = [];
therm(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 = ['geotailplasma',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.dat'];	% Contruct file names
ssoutfive = ['geotailplasma',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.mat'];	% Contruct file names
outputtxt = ['Output plasma file: ', ssoutfour];
disp(outputtxt);
[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);a		% 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];
vxgsm = zeros(1,newlll);
vygsm = zeros(1,newlll);
vzgsm = 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
	vectortwo = [xgse(jj),ygse(jj),zgse(jj)];	% Make a vector
	vectwoout = mat9vec3(RF,vectortwo);		% Transfor from gse to gsm
	xgsm(jj) = vectwoout(1);			% Sort our vector components
	ygsm(jj) = vectwoout(2);			% Sort our vector components
	zgsm(jj) = vectwoout(3);			% Sort our vector components
	if (vmag(jj) ~= 1.0e+34),			% repeat for velocity and convert to XYZ
		vxgsm(jj) = vmag(jj)*sin(pie/2 - theta(jj))*cos(phi(jj)); 	% km/s
		vygsm(jj) = vmag(jj)*sin(pie/2 - theta(jj))*sin(phi(jj));  	% km/s
		vzgsm(jj) = vmag(jj)*cos(pie/2 - theta(jj)); 			% km/s
	else
		vxgsm(jj) = 1.0e+34;		% error flag
		vygsm(jj) = 1.0e+34;		% error flag
		vzgsm(jj) = 1.0e+34;		% error flag
	end
end
disp('Convertion from GSE to GSM done');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% degap requires a specific input construct input array
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nMtime   = Mtime';
ndensity = density';
ntemp   = therm';
nvxgsm   = vxgsm';
nvygsm   = vygsm';
nvzgsm   = vzgsm';
nxgsm    = xgsm';
nygsm    = ygsm';
nzgsm    = zgsm';
inarray  = [nMtime,ntemp,ndensity,nvxgsm,nvygsm,nvzgsm,nxgsm,nygsm,nzgsm];
clear vygsm vzgsm xgsm ygsm zgsm density therm;
        
outarray=degapint(inarray,flag,datenum('00:01:00'),startend,datenum('00:10:00'));
[maxval,lllmag] = max(outarray(:,1));		% Degap done get new size.
nMtime   = outarray(:,1);			% sort out data.
ntemp   = outarray(:,2);			% I understand it better this way.
ndensity = outarray(:,3);
nvxgsm   = outarray(:,4);
nvygsm   = outarray(:,5);
nvzgsm   = outarray(:,6);
nxgsm    = outarray(:,7);
nygsm    = outarray(:,8);
nzgsm    = outarray(:,9);

% Here we are removing spikes in the data. Originally this was
% associated with the "pull" up edges arounded the originally flagged
% regions, but now it just removes spikes.
is = find((nvxgsm > 500.0)|(nvxgsm<-1300.0));
nvxgsm(is) = 1e+34;
is = find((nvygsm > 600.0)|(nvygsm<-600.0));
nvygsm(is) = 1e+34;
is = find((nvzgsm > 600.0)|(nvzgsm<-600.0));
nvzgsm(is) = 1e+34;
is = find((ntemp > 3000.0)|(ntemp < 0.0));
ntemp(is) = 1e+34;
is = find((ndensity > 200.0)|(ndensity <= 0.0));
ndensity(is) = 1e+34;
is = find((nxgsm > 300.0)|(nxgsm < -300.0));
nvxgsm(is) = 1e+34;
is = find((nygsm > 300.0)|(nygsm < -300.0));
nvygsm(is) = 1e+34;
is = find((nzgsm > 350.0)|(nzgsm < -350.0));
nvzgsm(is) = 1e+34;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot out to check
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;                                     	% make fig window
plot(Mtime,vxgsm,'b',nMtime,nvxgsm,'k*');   	% plot original and interpolated
datetick('x','dd/mm/yy');			% show date in dd/mm/yy form
if (min(vxgsm) == max(vxgsm)),
        ylim([0,1]);            		% set limits to avoid plotting flags
else
	ylim([-1000,200]);          		% set limits to avoid plotting flags
end
clear Mtime vxgsm;                            	% 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(outarray(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
outputarray = [newday',newmonth',newyear',newhr',newmin',newsec',ntemp,ndensity,nvxgsm,nvygsm,nvzgsm,nxgsm,nygsm,nzgsm];
					% export output array to file.
ctout = 1;                              % set counter
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 %12.6e %12.6e\n',outputarray(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','ntemp','ndensity','nvxgsm','nvygsm','nvzgsm','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 nvxgsm nvygsm nvzgsm nxgsm nygsm nzgsm ndensity ntemp;
disp('Program complete');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%end