function format_isee1_plasma_4(flag)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function format_imp_mag_0
%
% purpose: format the isee 2 magnetometer files
% into matlab time format and Vx, Vy, Vz, (GSM)
% X, Y, Z (GSM RE) and Density and Temp.
% File produced for ISEE 1 magnetometer data.
% Output columns are Bx, By, Bz, X, Y, Z in GSM
% FLAG requires format YYYYMM (tear/month)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Given 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Jul 14, 1975 62860   64.19 -0.44839 -1.3184 -3.1557 -3.3912 -2.8834 -0.94964 -1.8924 0.087856
format long e;							% use the long format
k_b 	       = 1.380622e-23;					% boltzmans const
m_p            = 1.672e-27;					% mass of a proton.
months         = strvcat('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec');	% String of month
%daysofmonths  = strvcat('31', '28', '31', '30', '31', '30', '31', '31', '30', '31', '30', '31');
%daysofmonths2 = strvcat('31', '29', '31', '30', '31', '30', '31', '31', '30', '31', '30', '31');
daysofmonths   = [31 28 31 30 31 30 31 31 30 31 30 31];		% day of year, normal year
daysofmonths2  = [31 29 31 30 31 30 31 31 30 31 30 31];		% day of year, leap year.
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. (Mag and Plasma)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yyyy=num2str(yy);						% Convert year/month to string	
ssone = ['isee1_plasma_',yyyy,'.txt'];				% Get file name
outputtxt = [' Input file: ', ssone];				% Check you did it correctly
	disp(outputtxt);
[fid,message] = fopen(ssone,'r');				% Error test
if fid == -1
	disp(message);
end
ssthree = ['isee1_mag_',yyyy,'.dat'];				% Get file name
outputtxt = [' Input file: ', ssthree];				% Check you did it correctly
	disp(outputtxt);
[fidmag,message] = fopen(ssthree,'r');
if fidmag == -1
	disp(message);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Enter the input plasma file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alldata = fscanf(fid,['%d %f %f %f %f %f %f %f \n %f %f\n'],[10,inf]);
% get the first line and construct the start Matlab time.
if (length(alldata) == 0),			% if file empty, then make flag file
 	yymmdd(1:2)    = (yy-190000)*100 +1;
	secsoday(1:2)  = 43200;
	UTtime(1:2)    = 120000;
	density(1:2)   = eflag;
	speed(1:2)     = eflag;
	flowazm(1:2)   = eflag;
	flowlat(1:2)   = eflag;
	tempmin(1:2)   = eflag;
	tempmax(1:2)   = eflag;
	alphafrac(1:2) = eflag;
else						% Location determined by getting ascii file header
	yymmdd   = alldata(1,:);		% in yymmdd format
	secsoday = alldata(2,:); 		% secs of the days
	UTtime   = alldata(3,:);		% stupid hhmmss.ss format
	density  = alldata(4,:);		% assume #/cm-3
	speed    = alldata(5,:);		% assume km/s
	flowazm  = alldata(6,:);		% flow azimuth (degs., with 0 deg. corresponding 
						% to flow from the sun [corrected for aberration] 
						% and positive azimuths corresponding to flow 
						% toward dawn),

	flowlat  = -1.0*alldata(7,:);	 	% flow latitude (degs., with positive latitudes 
		    				% corresponding to flow toward the south),
		    				% -1 is to flip data up into northern hemisphere  
		    				% for positive flow lat in input data file.
	tempmin  = alldata(8,:);   		% assume this is deg kelvin
	tempmax  = alldata(9,:);  		% assume this is deg kelvin
	alphafrac = alldata(10,:);
end

yrs  = double(int32(yymmdd/10000));		% get year from str
mons = double(int32(yymmdd/100)) - yrs*100;	% get mnoths from str
days = yymmdd - yrs*10000 - mons*100;		% get days
hrs  = double(int32(UTtime/10000));		% get hours
mins = double(int32(UTtime/100)) - hrs*100;	% get mins
secs = UTtime-hrs*10000-mins*100;		% get secss
yrs  = yrs + 1900; 				% year comes in as YY not YYYY
Mtime   = datenum(yrs,mons,days,hrs,mins,secs);	% get matlab time

flowazm = 3.141592654*(flowazm/180.0);		% get flow angle as deg we have rad
flowlat = 3.141592654*(flowlat/180.0);		% get flow angle as deg we have rad
[vxsctemp,vysctemp,vzsc]=sph2cart(flowazm,flowlat,speed);
vxsc    = -1.0*vxsctemp;			% convert to m/s
vysc    = -1.0*vysctemp;			% convert to m/s
clear alldata vxsctemp vysctemp;		% Free RAm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Enter the input plasma file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
allmagdata = fscanf(fidmag,['%d %d %d %d %d %d %d %f %f %f %f %f %f %f %f %f %f %f %f %f \n'],[20,inf]);
myrs     = allmagdata(1,:);			% Location determined by getting ascii
mmons    = allmagdata(2,:);			% file header
mdays    = allmagdata(3,:);
mhrs     = allmagdata(4,:);
mmins    = allmagdata(5,:);
msec     = allmagdata(6,:);
mmsec    = allmagdata(7,:);
bxgsm    = allmagdata(12,:);
bygsm    = allmagdata(13,:);
bzgsm    = allmagdata(14,:);
xgsm     = allmagdata(18,:);
ygsm     = allmagdata(19,:);
zgsm     = allmagdata(20,:);
Mmagtime = datenum(myrs,mmons,mdays,mhrs,mmins,msec);
clear allmagdata;
theirflag=1.0e+25;
index=find((xgsm>=theirflag) | (ygsm>=theirflag) | (zgsm>=theirflag) | (bxgsm>=theirflag) | (bygsm>=theirflag) | (bzgsm>=theirflag) );
%keyboard
%% 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)  = [];
bxgsm(index) = [];
bygsm(index) = [];
bzgsm(index) = [];
Mmagtime(index) = [];
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check what month we have
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
monthdoing = mons(1);
outputtxt = [' Working on month: ', num2str(monthdoing)];
	disp(outputtxt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert Vgse to Vgsm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Convert Velocity SC and position GSE to velocity and position GSM');
lenplasma = length(yrs);
vxgsm = zeros(1,lenplasma);
vygsm = zeros(1,lenplasma);
vzgsm = zeros(1,lenplasma);
%RF = [1 2 3; 4 5 6; 7 8 9];
for jj=1:1:lenplasma
	[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
	tempmax(jj) = (sqrt( (3*k_b*0.5*(tempmin(jj)+tempmax(jj)))/m_p ))/1000; % convert from Kelvin to km/s
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine output file name
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ssout    = ['isee1plasma',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.dat'];	% Contruct file names
ssouttwo = ['isee1plasma',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.mat'];	% Contruct file names
outputtxt = [' Output plasma file: ', ssout];
	disp(outputtxt);
% open output file
[fidout,message] = fopen(ssout,'w');		% open output file
if fidout == -1					% error check
	disp(message);
end
ssoutfour = ['isee1mag',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.dat'];	% Contruct file names
ssoutfive = ['isee1mag',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 strin
							% degap requires a specifiv input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% degap requires a specific input construct input array
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nMtime   = Mtime';					% turn row into column
nvxgsm   = vxgsm';
nvygsm   = vygsm';
nvzgsm   = vzgsm';
nden     = density';
ntempmax = tempmax';
inarray  = [nMtime,nvxgsm,nvygsm,nvzgsm,nden,ntempmax];
%lenplasma = length(yrs);
clear vygsm vzgsm tempmin tempmax density speed flowazm flowlat;	%free up some RAM.

% run the degap program which also interpolates to the minute
outarray=degapint(inarray,eflag,datenum('00:01:00'),startend,datenum('00:10:00'));
[maxval,lll] = max(outarray(:,1));			% Degap done get new size.
nMtime = outarray(:,1);					% sort out data.
nvxgsm = outarray(:,2);					% I understand it better this way.
nvygsm = outarray(:,3);
nvzgsm = outarray(:,4);
nden   = outarray(:,5);
ntempmax  = outarray(:,6);
% 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) = 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((ntempmax > 3000.0)|(ntempmax < 0.0));
ntempmax(is) = eflag;
is = find((nden > 200.0)|(nden <= 0.0));
nden(is) = eflag;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now magnetic field data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nMmagtime  = Mmagtime';				% turn row into column
nbxgsm     = bxgsm';
nbygsm     = bygsm';
nbzgsm     = bzgsm';
nxgsm      = xgsm';
nygsm      = ygsm';
nzgsm      = zgsm';
inmagarray = [nMmagtime,nbxgsm,nbygsm,nbzgsm,nxgsm,nygsm,nzgsm];
clear bygsm bzgsm xgsm ygsm zgsm;

outmagarray=degapint(inmagarray,eflag,datenum('00:01:00'),startend,datenum('00:10:00'));
[maxmagval,lllmag] = max(outmagarray(:,1));		% Degap done get new size.
nMmagtime = outmagarray(:,1);				% sort out data.
nbxgsm    = outmagarray(:,2);				% I understand it better this way.
nbygsm    = outmagarray(:,3);
nbzgsm    = outmagarray(:,4);
nxgsm     = outmagarray(:,5);
nygsm     = outmagarray(:,6);
nzgsm     = outmagarray(:,7);
is = find((nbxgsm > 2000.0)|(nbxgsm<-2000.0));		
nvxgsm(is) = eflag;
is = find((nbygsm > 2000.0)|(nbygsm<-2000.0));
nvygsm(is) = eflag;
is = find((nbzgsm > 2000.0)|(nbzgsm<-2000.0));
nvzgsm(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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% another check
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (length(nMmagtime) ~= length(nMtime)),
	disp('Mtime length of two different interpolated arrays does not match!!');
	keyboard;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot out to check
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;						% make fig window
subplot(2,1,1);
plot(Mtime,vxgsm,'b',nMtime,nvxgsm,'k*');	% plot original and interpolated
datetick('x','dd/mm/yyyy');			% 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
subplot(2,1,2);
plot(Mmagtime,bxgsm,'b',nMmagtime,nbxgsm,'k*'); % plot original and interpolated
datetick('x','dd/mm/yyyy');			% show date in dd/mm/yy form
if (min(bxgsm) == max(bxgsm)),
	ylim([0,1]);            		% set limits to avoid plotting flags
else
	ylim([min(bxgsm),min(bxgsm)+1000]);     % 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,lll);
newmon  = zeros(1,lll);
newday  = zeros(1,lll);
newhr   = zeros(1,lll);
newmin  = zeros(1,lll);
newsec  = zeros(1,lll);
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
	ctout = ctout + 1;				% incriment the counter
	atxt = ['Plasma 1: length: ',num2str(lll),' ctout: ',num2str(ctout)];	% text out to check
	disp(atxt);
end

disp('Making output array');
% contruct output array
outputarray = [newday',newmonth',newyear',newhr',newmin',newsec',ntempmax,nden,nvxgsm,nvygsm,nvzgsm,nxgsm,nygsm,nzgsm]; 
outputmagarray = [newday',newmonth',newyear',newhr',newmin',newsec',nbxgsm,nbygsm,nbzgsm,nxgsm,nygsm,nzgsm]; 
					% export output array to file.
ctout = 1;                              % set flag
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,:));
	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 = ['Plasm and Maga 2: length: ',num2str(lll),' ctout: ',num2str(ctout)];
	ctout = ctout + 1;              % incriment the counter
end
%keyboard
disp('Closing output file.');
status = fclose (fidout);	
status = fclose (fidmagout);	
    
% read out to matlab binary file.
disp('Making matlab binary files.');
save(ssouttwo,'nMtime','ntempmax','nden','nvxgsm','nvygsm','nvzgsm','nxgsm','nygsm','nzgsm');
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 nvxgsm nvygsm nvzgsm;

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