function format_isee3_plasma_4(flag)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function format_imp_mag_0
% author: J.M. Weygand
% purpose: format the ISEE 3 plasma files
% into matlab time format and vx, vy, vz, (GSM)
% X, Y, Z (GSM RE), Density and Vthermal.
% output is a ascii file and a matlab binary file
% FLAG has the format YYYYMM (year/month)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Givens
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
format long e;							% use the long format
k_b   = 1.380622e-23;     					% boltzmann constant
m_p   = 1.672e-27;          					% mass of a proton
eflag = 1.0e+34;						% errror flag
months        = strvcat('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec');	% str of month
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.
doyn = [1,32,60,91,121,152,182,213,244,274,305,335];		% day of year, first each mont, normal year
doyl = [1,32,61,92,122,153,183,214,245,275,306,336];		% day of year, first each mont, leap year
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load in data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yy=flag;							% Get year and month
outputtxt = [' Working on year: ', num2str(yy)]; 		% check inputs
	disp(outputtxt);
if (double(int32(flag/100)) >= 1984 ), 				% check again the inputs (ISEE3 mag has two versions)
	disp('Your using the wrong program! This is only good for year < 1983!');
	keyboard;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine file name and open input file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yyyy=num2str(yy);      						% num to string to get input files
ssone = ['i3_sw_prot_mom_',yyyy(3:6),'.txt'];			% Get file name
outputtxt = [' Input file: ', ssone];				% Check you did it correctly
	disp(outputtxt);
[fidp,message] = fopen(ssone,'r');				% Error test
if fidp == -1
	disp(message);
end
	
ssthree = [yyyy,'_isee3_mag03_1min.asc'];   			% make input file
outputtxt = [' Input file: ', ssthree]; 			% check input 
	disp(outputtxt);
[fidmag,message] = fopen(ssthree,'r');   			% open the file and check errors
if fidmag == -1
	disp(message);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Enter the input plasma file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
allmagdata = fscanf(fidmag,['%d %d %d %d %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n'],[21,inf]);
len = length(allmagdata(1,:));               % Determine array length
% sort data
myrs  = allmagdata(1,:);			% Location determined by getting ascii
doy   = allmagdata(2,:);			% file header
mhrs  = allmagdata(3,:);
mmins = allmagdata(4,:);
bxsc  = allmagdata(5,:);
bysc  = allmagdata(6,:);
bzsc  = allmagdata(7,:);
xgse  = allmagdata(19,:);
ygse  = allmagdata(20,:);
zgse  = allmagdata(21,:);

mdays    = zeros(1,len);     			% saves time 
mons     = zeros(1,len);     			% saves time
lyear = myrs(1);
disp('Getting Days and Months from DOY!'); 	% bob wants YYYY MM DD, but we have yyyy doy

Mmagtime = datenum(myrs,1,doy,mhrs,mmins,0); 	% get matlab time
clear allmagdata;  				% clear some RAM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input plasma data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
allplasmadata = fscanf(fidp,['%d %f %f %f %f %f %f %f\n'],[8,inf]);
yymmdds = allplasmadata(1,:);			% in yymmdd format
totsecs = allplasmadata(2,:);			% sedons of the day
hhmmss  = allplasmadata(3,:);			% stupid hhmmss.ss format
density = allplasmadata(4,:);			% assume #/cm^-3
vmag    = allplasmadata(5,:); 			% bulk speed
vygse   = -1.0*allplasmadata(6,:); % flow azimuth km/s ??
% raw ascii data contains V_azmuth (km/s)
% this is pretty much V_y (gse) to withi fraction
% of a degree.
tempmin = allplasmadata(7,:);			% assume in degrees kelvin
tempmax = allplasmadata(8,:);			% assume in degrees kelvin

lenplasma = length(allplasmadata(6,:));
vzgse = zeros(1,lenplasma);
vxgse = -1.0*sqrt(vmag.*vmag - vygse.*vygse);
% Vmag^2 = Vx^2 + Vy^2 + Vz^2
% have no Vz component so assume zero
% have Vazmuth but assume -vygse
% hence we get the about equation
yrs  = double(int32(yymmdds/10000));
mons = double(int32(yymmdds/100)) - yrs*100;
days = yymmdds - yrs*10000 - mons*100;
hrs  = double(int32(hhmmss/10000));
mins = double(int32(hhmmss/100)) - hrs*100;
secs = hhmmss-hrs*10000-mins*100;
yrs = yrs + 1900; % year comes in as YY not YYYY
Mplasmatime = datenum(yrs,mons,days,hrs,mins,secs);
theirflag = 1.0e+25;
index=find((xgse>=9999.9) | (ygse>=9999.9) | (zgse>=9999.9) | (bxsc>=999.9) | (bysc>=999.9) | (bzsc>=999.9));	%% Change flag value
%% 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) = [];
bxsc(index) = [];
bysc(index) = [];
bzsc(index) = [];
Mmagtime(index) = [];
len = length(Mmagtime);
theirpflag = 1.0e+25;
index=find((vmag>=theirpflag) | (density>=theirpflag) | (tempmin>=theirpflag) | (tempmax>=theirpflag) | (vygse>=theirpflag));	%% Change flag value
vmag(index) = [];
density(index) = [];
tempmin(index) = [];
tempmax(index) = [];
vygse(index) = [];
vxgse(index) = [];
vzgse(index) = [];
Mplasmatime(index) = [];
lenp = length(Mplasmatime);
    
%%%keyboard
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert Vgse to Vgsm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Convert Velocity SC and position GSE to velocity and position GSM');
vxgsm = zeros(1,lenp);						% save time
vygsm = zeros(1,lenp);						% save time
vzgsm = zeros(1,lenp);						% save time
%RF = [1 2 3; 4 5 6; 7 8 9];
for jj=1:1:lenp
	[RF,RI] = gse2gsm(Mplasmatime(jj));			% builds the matrix
	vector = [vxgse(jj),vygse(jj),vzgse(jj)];		
	vecout = mat9vec3(RF,vector);				% make a vector
	vxgsm(jj) = vecout(1);					% sort our vector components
	vygsm(jj) = vecout(2);					% sort our vector components
	vzgsm(jj) = eflag;          				% vecout(3); not real, flag it.
	temp(jj) = (sqrt( (3*k_b*0.5*(tempmin(jj)+tempmax(jj)))/m_p ))/1000; % convert from Kelvin to km/s
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check what month we have
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
monthdoing = mons(1);  						% check, check and recheck
outputtxt = [' Working on month: ', num2str(monthdoing)];	% examine output
	disp(outputtxt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine output file name
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ssout    = ['isee3plasma',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.dat'];	% Contruct file names
ssouttwo = ['isee3plasma',num2str(yrs(1)),num2str(monthdoing,'%2.2d'),'.mat'];	% Contruct file names
outputtxt = [' Output plasma file: ', ssout];
	disp(outputtxt);

[fidpout,message] = fopen(ssout,'w');				% open output file
if fidpout == -1						% error check
	disp(message);
end
ssoutfour = ['isee3mag',num2str(myrs(1)),num2str(monthdoing,'%2.2d'),'.dat'];	% Contruct file names
ssoutfive = ['isee3mag',num2str(myrs(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    					% check if they really open.
	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(myrs(1)),' ','00:00:00.00'];
startend(1) = datenum(startyyymmddhhmmss);			% Determine start time of month and make string
lyear = myrs(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(myrs(1)),' ','23:59:00.00'];
else
	endyyymmddhhmmss = [num2str(daysofmonths(monthdoing)),'-',months(monthdoing,:),'-',num2str(myrs(1)),' ','23:59:00.00'];
end
startend(2) = datenum(endyyymmddhhmmss);			% Determine end time of month and make string
	                                                        % degap requires a specifiv input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert Mag and position from GSE to GSM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Convert Magnetic SC and position GSE to magnetic field and position GSM');
%RF = [1 2 3; 4 5 6; 7 8 9];
bxgsm = zeros(1,len); 						% these statements save computation time.
bygsm = zeros(1,len);
bzgsm = zeros(1,len);
xgsm = zeros(1,len);
ygsm = zeros(1,len);
zgsm = zeros(1,len);
for jj=1:1:len   						% convert SC about = GSE to GSM coords
	[RF,RI] = gse2gsm(Mmagtime(jj));			% build matrix
	vector = [bxsc(jj),bysc(jj),bzsc(jj)];			% build vector
	vecout = mat9vec3(RF,vector);
	bxgsm(jj) = vecout(1);					% sort our vector components
	bygsm(jj) = vecout(2);
	bzgsm(jj) = vecout(3);
	vectortwo = [xgse(jj),ygse(jj),zgse(jj)];		% build vector
	vectwoout = mat9vec3(RF,vectortwo);
	xgsm(jj) = vectwoout(1);				% sort our vector components
	ygsm(jj) = vectwoout(2);
	zgsm(jj) = vectwoout(3);
end
disp('Convertion from GSE to GSM done');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% degap requires a specific input construct input array
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nMmagtime = Mmagtime';   					% Bobs interpolation program needs a specific format
nbxgsm    = bxgsm';						% turn row into column
nbygsm    = bygsm';
nbzgsm    = bzgsm';
nxgsm     = xgsm';
nygsm     = ygsm';
nzgsm     = zgsm';
inmagarray = [nMmagtime,nbxgsm,nbygsm,nbzgsm,nxgsm,nygsm,nzgsm];
clear bygsm bzgsm xgsm ygsm zgsm;  				% free RAM

outmagarray=degapint(inmagarray,eflag,datenum('00:01:00'),startend,datenum('00:10:00'));
[maxmagval,lllmag] = max(outmagarray(:,1));
nMmagtime = outmagarray(:,1);     			% put in to something more readable for this stupid programmer
nbxgsm    = outmagarray(:,2);
nbygsm    = outmagarray(:,3);
nbzgsm    = outmagarray(:,4);
nxgsm     = outmagarray(:,5);
nygsm     = outmagarray(:,6);
nzgsm     = outmagarray(:,7);
% 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((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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% degap requires a specific input construct input array
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nMplasmatime = Mplasmatime';						% turn row into cloumn
nvxgsm = vxgsm';
nvygsm = vygsm';
nvzgsm = vzgsm';
nden   = density';
ntemp  = temp';
inarray = [nMplasmatime,nvxgsm,nvygsm,nvzgsm,nden,ntemp];
%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));
%keyboard
nMtime = outarray(:,1);
nvxgsm = outarray(:,2);
nvygsm = outarray(:,3);
nvzgsm = outarray(:,4);
nden   = outarray(:,5);
ntemp  = outarray(:,6);
is = find((nvxgsm > 2000.0)|(nvxgsm<-2000.0));
nvxgsm(is) = eflag;
is = find((nvygsm > 2000.0)|(nvygsm<-2000.0));
nvygsm(is) = eflag;
is = find((nvzgsm > 2000.0)|(nvzgsm<-2000.0));
nvzgsm(is) = eflag;

is = find((nden > 200.0)|(nden < 0.0));
nden(is) = eflag;
is = find((ntemp > 3000.0)|(ntemp < 0.0));
ntemp(is) = eflag;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot out to check
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;							% make fig window
subplot(2,1,1);
plot(Mplasmatime,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([-500.0,500.0]);          			% 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');
if (min(bxgsm) == max(bxgsm)),
	ylim([0,1]);            			% set limits to avoid plotting flags
else
	ylim([-200,200]);          			% 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]; 
outputplasmaarray = [newday',newmonth',newyear',newhr',newmin',newsec',ntemp,nden,nvxgsm,nvygsm,nvzgsm,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,:));
	count = fprintf(fidpout,'%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',outputplasmaarray(ctout,:));
	atxt = ['Plasma and 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,'nMmagtime','nbxgsm','nbygsm','nbzgsm','nxgsm','nygsm','nzgsm');
save(ssouttwo,'nMmagtime','ntemp','nden','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 nden ntemp nbxgsm nbygsm nbzgsm;

status = fclose (fidmag);          % close input files
status = fclose (fidp);          % close input files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of for loop for years (Yeah, that was a pipe dream).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%end