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+34;
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 expected: ', 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.');
%keyboard
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
%%keyboard
if ~isnan(flag),
% convert the flags to nans
	nMtime(find(nMtime==flag))    = nan;
	nbxgsm(find(nbxgsm==flag))    = nan;
	nbygsm(find(nbygsm==flag))    = nan;
	nbzgsm(find(nbzgsm==flag))    = nan;
	nxgsm(find(nxgsm==flag))      = nan;
	nygsm(find(nygsm==flag))      = nan;
	nzgsm(find(nzgsm==flag))      = nan;
	ntherm(find(ntherm==flag))    = nan;
	ndensity(find(ndensity==flag)) = nan;
	nvxgsm(find(nvxgsm==flag))    = nan;
	nvygsm(find(nvygsm==flag))    = nan;
	nvzgsm(find(nvzgsm==flag))    = nan;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Elimnate points in the magnetosphere
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the next three constants are taken from Slavin and Holzer, JGR v86
% pp 11401--11418 table 4. They are part of a model to determine the
% location of the bowshock. we do not want to include points inside the bow
% shock in the triggering study
eccent = 1.19; 		% eccentricity
LL     = 23.1; 		% stand off distance of magnetopause
xnaught = 3.0;		% constant? in Re
%%%% determine the bowshock location
ii = 1;
for deg = -145:1:145
	theta = 3.141592*(deg/180.0);
	rad(ii) = LL/(1 + eccent*cos(theta));
	yy(ii) = rad(ii)*sin(theta);
	xx(ii) = rad(ii)*cos(theta) + xnaught;
	ii = ii+1;
end
for ii = 1:1:length(nxgsm)
    xbow(ii) = interp1(yy,xx,nygsm(ii),'nearest');  
end
xbow(isnan(nygsm))      = nan;
index = find(nxgsm <= xbow');
    nMtime(index)    = nan;
	nbxgsm(index)    = nan;
	nbygsm(index)    = nan;
	nbzgsm(index)    = nan;
	nxgsm(index)      = nan;
	nygsm(index)      = nan;
	nzgsm(index)      = nan;
	ntherm(index)    = nan;
	ndensity(index) = nan;
	nvxgsm(index)    = nan;
	nvygsm(index)    = nan;
	nvzgsm(index)    = nan;
figure;
plot(nMtime,nxgsm,'b',nMtime,xbow,'k',nMtime(index),nxgsm(index),'r*','LineWidth',3);
ylim([-50,300]);


%keyboard

%%%%% 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.');