function propagate_all_1(sat,yyyymm, yyyymmb)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unction propagate_all_0(sat,instr,yyyymm)
% Purpose: load and propagate solar wind data to
% bow show at 15 Re.
% Author: J.M. weygand
% SAT = spacecraft
% INSTR = instrument
% YYYYMM = year and mm of data file to process
% YYYYMMB = year and month of data file before. Data from the
% prceeding file is required for the propagation to fill in
% the end effects resulting from the propagation forward
% in time. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Givens
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bow       = 15.0;      % units: Re
re        = 6378.16;   % units: km
secperday = 86400;     % number of seconds per day. Units: sec/day
format long e;         % format for all numbers
Toff      = 1.0;
flag      = 1.0e+25;
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];
daysofmonths2 = [31 29 31 30 31 30 31 31 30 31 30 31];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yyyymmbefor    = num2str(yyyymmb);
yyyymmtxt      = num2str(yyyymm);
inputfileone   = [sat,'plasma',yyyymmbefor,'.mat'];
inputfiletwo   = [sat,'mag',yyyymmbefor,'.mat'];
inputfilethree = [sat,'plasma',yyyymmtxt,'.mat'];
inputfilefour  = [sat,'mag',yyyymmtxt,'.mat'];

inputtxt = [' Input files xpected: ', inputfileone, ' and ', inputfiletwo, ' and ', inputfilethree, ' and ', inputfilefour];
disp(inputtxt);
yrb = double(int32(yyyymmb/100));
monb = yyyymmb-yrb*100;
yr = double(int32(yyyymm/100));
mon = yyyymm-yr*100;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get start and end of the month.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
startyyymmddhhmmss = [num2str(1,'%2.2d'),'-',months(mon,:),'-',num2str(yr),' ','00:00:00.00'];
startend(1) = datenum(startyyymmddhhmmss);
lyear = yr;		% 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(mon)),'-',months(mon,:),'-',num2str(yr),' ','23:59:00.00'];
else
     endyyymmddhhmmss = [num2str(daysofmonths(mon)),'-',months(mon,:),'-',num2str(yr),' ','23:59:00.00'];
end
startend(2) = datenum(endyyymmddhhmmss);   

disp('Propagating Bxgsm values.');
if (yyyymmb == -1),
        disp('No previous file available.');
else
    inplasmafileb = load(inputfileone); 
    inmagfileb = load(inputfiletwo);
    lll = size(inmagfileb.nMtime);
end
%keyboard
inplasmafile = load(inputfilethree); 
inmagfile = load(inputfilefour);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fill arrays
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (yyyymmb > 1900),
    nMtime   = cat(1,inmagfileb.nMtime( lll(1)-180:lll(1) ,1),inmagfile.nMtime);
    nbxgsm   = cat(1,inmagfileb.nbxgsm( lll(1)-180:lll(1) ,1),inmagfile.nbxgsm);
    nbygsm   = cat(1,inmagfileb.nbygsm( lll(1)-180:lll(1) ,1),inmagfile.nbygsm);
    nbzgsm   = cat(1,inmagfileb.nbzgsm( lll(1)-180:lll(1) ,1),inmagfile.nbzgsm);
    nxgsm    = cat(1,inmagfileb.nxgsm( lll(1)-180:lll(1) ,1),inmagfile.nxgsm);
    nygsm    = cat(1,inmagfileb.nygsm( lll(1)-180:lll(1) ,1),inmagfile.nygsm);
    nzgsm    = cat(1,inmagfileb.nzgsm( lll(1)-180:lll(1) ,1),inmagfile.nzgsm);

    ntherm   = cat(1,inplasmafileb.ntherm( lll(1)-180:lll(1) ,1),inplasmafile.ntherm);
    ndensity = cat(1,inplasmafileb.ndensity( lll(1)-180:lll(1) ,1),inplasmafile.ndensity);
    nvxgsm   = cat(1,inplasmafileb.nvxgsm( lll(1)-180:lll(1) ,1),inplasmafile.nvxgsm);
    nvygsm   = cat(1,inplasmafileb.nvygsm( lll(1)-180:lll(1) ,1),inplasmafile.nvygsm);
    nvzgsm   = cat(1,inplasmafileb.nvzgsm( lll(1)-180:lll(1) ,1),inplasmafile.nvzgsm);
else
    nMtime   = inmagfile.nMtime;
    nbxgsm   = inmagfile.nbxgsm;
    nbygsm   = inmagfile.nbygsm;    
    nbzgsm   = inmagfile.nbzgsm;
    nxgsm    = inmagfile.nxgsm;
    nygsm    = inmagfile.nygsm;
    nzgsm    = inmagfile.nzgsm;

    ntherm   = inplasmafile.ntherm;
    ndensity = inplasmafile.ndensity;
    nvxgsm   = inplasmafile.nvxgsm;
    nvygsm   = inplasmafile.nvygsm;
    nvzgsm   = inplasmafile.nvzgsm;
end
if ~isnan(flag),
% convert the flags to nans
	nMtime(nMtime==flag)    = nan;
	nbxgsm(nbxgsm==flag)    = nan;
	nbygsm(nbygsm==flag)    = nan;
	nbzgsm(nbzgsm==flag)    = nan;
	nxgsm(nxgsm==flag)      = nan;
	nygsm(nygsm==flag)      = nan;
	nzgsm(nzgsm==flag)      = nan;
	ntherm(ntherm==flag)    = nan;
	ndensity(ndensity==flag) = nan;
	nvxgsm(nvxgsm==flag)    = nan;
	nvygsm(nvygsm==flag)    = nan;
	nvzgsm(nvzgsm==flag)    = nan;
end


%%%%% Create the propagation array specific for Bob's program.
disp('Creating propagation array.');
pMtime = inmagfile.nMtime;
datin.dnup = nMtime;
datin.Rup  = [nxgsm,nygsm,nzgsm];
datin.Bup  = [nbxgsm,nbygsm,nbzgsm];
bowshock = 15.0;
pxgsm = bowshock*ones(size(nxgsm));
datin.Rdn  = [pxgsm,nygsm,nzgsm];
datin.Vup  = abs(nvxgsm);

disp('Propagation array completed.');
% Propagate the data. Note we must change the data to be propagated.
disp('Propagate data.');
%%%%%%%%%%%%%Magnetic field and Time%%%%%%%%%%%%%%%
datin.dat  = nbzgsm;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
pMtime = temp(:,1);
pbzgsm = temp(:,2);
datin.dat  = nbxgsm;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
pbxgsm = temp(:,2);
datin.dat  = nbygsm;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
pbygsm = temp(:,2);
%%%%%%%%%%%%%Position%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
datin.dat  = nxgsm;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
pxgsm = temp(:,2);
datin.dat  = nygsm;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
pygsm = temp(:,2);
datin.dat  = nzgsm;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
pzgsm = temp(:,2);
%%%%%%%%%%%%%Therm and Density%%%%%%%%%%%%%%%%%%%%%
datin.dat  = ntherm;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
ptherm = temp(:,2);
datin.dat  = ndensity;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
pdensity = temp(:,2);
%%%%%%%%%%%%%Velocity%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
datin.dat  = nvxgsm;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
pvxgsm = temp(:,2);
datin.dat  = nvygsm;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
pvygsm = temp(:,2);
datin.dat  = nvzgsm;
[temp,fT] = propsw(datin,'Spherical',Toff,startend);
pvzgsm = temp(:,2);

%keyboard
disp('Propagation of data done.');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
%plot(nMtime,nbzgsm,'b',temp(:,1),temp(:,2)-20.0,'r');
plot(nMtime,nbzgsm,'b',pMtime,pbzgsm-20.0,'r');
ylim([-40,20]);
datetick('x','dd/mm/yy');
%keyboard
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read out to an ascii file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ctout = 1;				% set flag
temp = size(pMtime);
lll = temp(1);
disp('Fixing rounding error in Time.');
newyear = zeros(1,lll);
newmonth  = 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(pMtime(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('Rounding error fixed.');
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % make output file names
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 outmagAfile = [sat,'magP',yyyymmtxt,'.dat'];
 outmagMfile = [sat,'magP',yyyymmtxt,'.mat'];
 outputtxt = [' Output plasma file: ', outmagAfile, ' and ', outmagMfile];
 disp(outputtxt);
 outplasmaAfile = [sat,'plasmaP',yyyymmtxt,'.dat'];
 outplasmaMfile = [sat,'plasmaP',yyyymmtxt,'.mat'];
 outputtxt = [' Output plasma file: ', outplasmaAfile, ' and ', outplasmaMfile];
 disp(outputtxt);
 % open output file
 [fidmagout,message] = fopen(outmagAfile,'w');
 if fidmagout == -1
                disp(message);
 end
 [fidplasmaout,message] = fopen(outplasmaAfile,'w');
 if fidplasmaout == -1
                disp(message);
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make new position arrays
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pxgsm = zeros(1,lll) + 15.00;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% output ascii and Mfiles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

outputarray    = [newday',newmonth',newyear',newhr',newmin',newsec',ptherm,pdensity,pvxgsm,pvygsm,pvzgsm,pxgsm',pygsm,pzgsm]; 
outputmagarray = [newday',newmonth',newyear',newhr',newmin',newsec',pbxgsm,pbygsm,pbzgsm,pxgsm',pygsm,pzgsm]; 
ctout = 1;                              % set flag

disp('Writing ASCII propagation files.');
while( (ctout) <= lll),
		count = fprintf(fidplasmaout,'%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 = ['Plasma and Mag 2: length: ',num2str(lll),' ctout: ',num2str(ctout)];
		ctout = ctout + 1;              % incriment the counter
end
disp('Closing output file.');
status = fclose (fidplasmaout);	
status = fclose (fidmagout);	
     
% read out to matlab binary file.
disp('Writing Mfiles.');
save(outplasmaMfile,'pMtime','ptherm','pdensity','pvxgsm','pvygsm','pvzgsm','pxgsm','pygsm','pzgsm');
save(outmagMfile,'pMtime','pbxgsm','pbygsm','pbzgsm','pxgsm','pygsm','pzgsm');
disp('Program complete.');