function format_imp_plasma_5(flag)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function format_imp_plasma_5
% author: J. Weygand
% creation: Mar 1, 2003
% purpose: format the imp 8 plasma files
% into matlab time format thermal speed, 
% density vx, vy, vz, (GSM)
% X, Y, Z (GSM RE)
% NOTE:  this is for file type number 2
% which has no vx, vy, vz in gse or gsm.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Given 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
doyn          = [1,32,60,91,121,152,182,213,244,274,305,335];	% day of year, normal year
doyl          = [1,32,61,92,122,153,183,214,245,275,306,336];	% day of year, leap year.
k_b           = 1.380622e-23;					% boltzzmans const
m_p           = 1.672e-27;					% mass of a proton.
eflag 	      = 1.0e+34;					% error flag 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load in data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yy=flag;							% Get year and month
outputtxt = [' Working on year: ', num2str(yy)];		% Double check
disp(outputtxt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine file name and open input file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yyyy=num2str(yy);						% Convert year/month to string
ssone = ['plasma_',yyyy,'.asc'];				% Get file name
outputtxt = [' Input file: ', ssone];				% Check you did it correctly
disp(outputtxt);
[fid,message] = fopen(ssone,'r');
if fid == -1
	disp(message);
	keyboard;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Enter the first few lines of the input file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alldata = fscanf(fid,['%d %d %d %d %d %d %d %d %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n'],[24,inf]);
% get all the data at once. This saves time.
if (length(alldata)==0), 			% if file empty, then make flag file
	yrstr        = num2str(yyyy);
	yrs(1:2)     = str2num(yrstr(1:4));
	mons(1:2)    = str2num(yrstr(5:6));
	days(1:2)    = 15;
	hrs(1:2)     = 0;
	mins(1:2)    = 0;
	secs(1:2)    = 0;
	xgse(1:2)    = eflag;
	ygsm(1:2)    = eflag;
	zgsm(1:2)    = eflag;
	vxgsm(1:2)   = eflag;
	vygsm(1:2)   = eflag;
	vzgsm(1:2)   = eflag;
	temp(1:2)    = eflag;
	density(1:2) = eflag;
else						% if files has data sort the data
	yrs     = alldata(3,:);			% Location determined by getting ascii
	doy     = alldata(4,:);			% file header
	hrs     = alldata(5,:);
	mins    = alldata(6,:);
	secs    = alldata(7,:);
	msecs   = alldata(8,:);
	xgsm    = alldata(10,:);
	ygsm    = alldata(13,:);
	zgsm    = alldata(14,:);
	speed   = alldata(15,:);
	temp    = alldata(16,:);
	density = alldata(17,:);
	EW_ang  = alldata(18,:);
	NW_ang  = alldata(19,:);
end

len  = length(yrs);
lyear = yrs(1);
disp(' Getting Days and Months from DOY!'); 		% bob wants YYYY MM DD, but we have yyyy doy
Mtime = datenum(yrs,1,doy,hrs,mins,secs);		%construct Matlab time.
[dummy,mmons,dummy,dummy,dummy,dummy]= datevec(Mtime(1));

clear alldata;						% free some RAM
lll = length(Mtime);
diff = Mtime(2:lll) - Mtime(1:lll-1);			% check for time error.
is  = find(diff<0);					% some file catenated together funnny.
if (length(is)>0)					% New files exist be sure your using them
	endtxt = datestr(Mtime(is),0);
	txt = ['Time error at: '];
	disp(strcat(txt,endtxt));
end
flag=1.0e+34;
index=find((xgsm>=9999.0) | (ygsm>=9999.0) | (zgsm>=9999.0) | (speed>=9999.0) | (EW_ang>=9999.0) | (NW_ang>=9999.0) | (density>=9999.0) | (temp>=9999.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.
xgsm(index)    = [];
ygsm(index)    = [];
zgsm(index)    = [];
speed(index)   = [];
EW_ang(index)  = [];
NW_ang(index)  = [];
density(index) = [];
temp(index)    = [];
Mtime(index)   = [];
len            = length(Mtime);
% change the flag. The flag is changed for the temp and not the others
% because there is so many flags for the temp but very little for the
% the others. if vxgse is missing then so usually is vygse and vzgse.
% Denisity is also rarely missing. Temperature has normally 10000's of
% flags so it seemed worth while to sacrifice data points on the ends
% inorder to keep more data in the flaged temp regions.
%index=find(temp>=99999.9);	%% Change flag value
%temp(index)=repmat(flag,1,length(index));        
    
%index=find(xgsm>=9999.);	%% Change flag value
%xgsm(index)=repmat(flag,1,length(index));        
%index=find(ygsm>=9999.);	%% Change flag value
%ygsm(index)=repmat(flag,1,length(index));        
%index=find(zgsm>=9999.);	%% Change flag value
%zgsm(index)=repmat(flag,1,length(index));        
%index=find(density>=9999.);	%% Change flag value
%density(index)=repmat(flag,1,length(index)) ;       
%index=find(temp>=9999.);	%% Change flag value
%temp(index)=repmat(flag,1,length(index));        
%index=find(speed>=9999.);	%% Change flag value
%speed(index)=repmat(flag,1,length(index));        
%index=find(EW_ang>=9999.);	%% Change flag value
%EW_ang(index)=repmat(flag,1,length(index));        
%index=find(NW_ang>=9999.);	%% Change flag value
%NW_ang(index)=repmat(flag,1,length(index));

monthdoing = mmons(1);						% what month are we really doing
outputtxt = [' Working on month: ', num2str(monthdoing)];	
disp(outputtxt);						% what month are we really doing display to screen 
flowazm = 3.141592654*(EW_ang/180.0);				% have degrees, want radians
flowlat = 3.141592654*(NW_ang/180.0);   			% have degrees, want radians
[vxsc,vysc,vzsc]=sph2cart(flowazm,flowlat,speed);
vxsc = -1.0*vxsc;						% points the opposite direction
vysc = -1.0*vysc;						% points the opposite direction
vzsc = -1.0*vzsc;						% points the opposite direction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert Vgse to Vgsm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Convert Velocity SC and position GSE to velocity and position GSM');
%RF = [1 2 3; 4 5 6; 7 8 9];	example of what it looks like
for jj=1:1:len
	[RF,RI] = gse2gsm(Mtime(jj));				% builds the matrix
	vector = [vxsc(jj),vysc(jj),vzsc(jj)];			% make a vector
	vecout = mat9vec3(RF,vector);				% transfor from gse to gsm
	vxgsm(jj) = vecout(1);					% sort our vector components
	vygsm(jj) = vecout(2);					% sort our vector components
	vzgsm(jj) = vecout(3);					% sort our vector components
end
disp('Convertion from GSE to GSM done');			% status report.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine output file name
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ssout    = ['imp8plasma',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.dat']; % Contruct file names
ssouttwo = ['imp8plasma',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.mat'] % Contruct file names;
outputtxt = [' Output file: ', ssout];
disp(outputtxt);

[fidout,message] = fopen(ssout,'w');				% open output file
if fidout == -1							% error check
	disp(message);
	keyboard;
end

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 specifiv input
ntemp    = temp';						% turn row into column
ndensity = density';
nMtime   = Mtime';
nvxgsm   = vxgsm';
nvygsm   = vygsm';
nvzgsm   = vzgsm';
nxgsm    = xgsm';	% xgse = xgsm
nygsm    = ygsm';
nzgsm    = zgsm';
inarray  = [nMtime,ntemp,ndensity,nvxgsm,nvygsm,nvzgsm,nxgsm,nygsm,nzgsm];
clear vygsm vzgsm xgsm ygsm zgsm;				%free up some RAM.
% run the degap program which also interpolates to the minute
outarray=degapint(inarray,1.0e25,datenum('00:01:00'),startend,datenum('00:10:00'));
[maxval,lll] = max(outarray(:,1));
nMtime       = outarray(1:lll,1);
ntemp        = outarray(1:lll,2);
ndensity     = outarray(1:lll,3);
nvxgsm       = outarray(1:lll,4);
nvygsm       = outarray(1:lll,5);
nvzgsm       = outarray(1:lll,6);
nxgsm        = outarray(1:lll,7);
nygsm        = outarray(1:lll,8);
nzgsm        = outarray(1:lll,9);
is = find((nvxgsm > 500.0)|(nvxgsm<-1300.0));
nvxgsm(is) = eflag;
is = find((nvygsm > 600.0)|(nvygsm<-600.0));
nvygsm(is) = eflag;
is = find((nvzgsm > 600.0)|(nvzgsm<-600.0));
nvzgsm(is) = eflag;
is = find((ntemp > 3000.0)|(ntemp < 0.0));
ntemp(is) = eflag;
is = find((ndensity > 200.0)|(ndensity <= 0.0));
ndensity(is) = eflag;
is = find((nxgsm > 300.0)|(nxgsm < -300.0));
nvxgsm(is) = eflag;
is = find((nygsm > 300.0)|(nygsm < -300.0));
nvygsm(is) = eflag;
is = find((nzgsm > 350.0)|(nzgsm < -350.0));
nvzgsm(is) = eflag;

%clear outarray;	%free up some RAM.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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([min(vxgsm),max(vxgsm)]);		% 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.');
while( (ctout) <= lll),
	% 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
	atxt = ['length: ',num2str(lll),' ctout: ',num2str(ctout),' day: ', num2str(newday(ctout)),' hr: ',num2str(newhr(ctout))];
							% text out to check
	disp(atxt);
	ctout = ctout + 1;				% incriment the counter
end
disp('Making output array');
% contruct output array

outputarray = [newday',newmonth',newyear',newhr',newmin',newsec',ntemp(1:lll), ndensity(1:lll), nxgsm(1:lll),nygsm(1:lll),nzgsm(1:lll),nvxgsm(1:lll),nvygsm(1:lll),nvzgsm(1:lll)]; 
							% export output array to file.
ctout = 1;

while( (ctout) <= lll),
	count = fprintf(fidout,'%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,:));
	ctout = ctout + 1;              % incriment the counter
end
disp('Closing output file.');
status = fclose (fidout);	
    
% read out to matlab binary file.
disp('Making matlab binary file.');
save(ssouttwo,'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;

status = fclose (fid);          % close input files
disp('Program complete');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of for loop for years (Yeah, that was a pipe dream).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%end