	function [clksv, a, e, E, f, i, omega, l, r] = bcast(t, eph, toe)
%	bcast	 [clksv, a, e, E, f, i, omega, l, r] = bcast(t, eph, toe)
%
%	Computes Kepelerian elements from broadcast message
%	Note: Takes a matrix of ephemeris input.
%
%	INPUTS		
%		t 		 = (tx1) vector of epochs
%		eph		 = (31xs) matrix of ephemeris values
%		toe 		 = (sx1) vector of ephemeris times
%
%	OUTPUTS
%		clksv	= (txs) Satellite clock offset
%		a	= (1xs) Semi-major axis
%		e	= (1xs) eccentricity
%		E	= (txs) Eccentric anomaly
%		f	= (txs) True anomaly
%		i	= (txs) inclination
%		omega	= (txs) argument of perigee
%		l	= (txs) Longitude of ascending node  
%		r 	= (txs) range including  cosine and sine terms
%
%	s: number of satellites
%	t: number of epochs

%Set constants, etc.
	
	week0 = 44244.0;
	GM2=1.996498184321739e+007;
	omegaE=7292115.1467e-11;
	j=size(t,1);
	vec2mat=ones(j,1);

%Calculate time difference between current time and time of ephemeris

	dt = (t*ones(1,length(toe)) - vec2mat*toe')*86400.0;
		
%Calculate time difference between current time and beginning of GPS week

	gpsweek = eph(22);
	t0 = week0 + (gpsweek*7);
	dt0 = (t - t0)*86400;
		
%Calculate semi-major axes

	a = vec2mat*(eph(11,:).*eph(11,:));

%Get eccentricities

	e = vec2mat*eph(9,:);

%Get orbital inclinations, rates of change

	i0 = vec2mat*eph(16,:);
	idot = vec2mat*eph(20,:);

%Get right ascensions of the ascending nodes, rates of change

	OMEGA0 = vec2mat*eph(14,:);
	OMEGAdot = vec2mat*eph(19,:);

%Get arguments of perigee

	omega0 = vec2mat*eph(18,:);

%Get mean motion differences;

	deltan = vec2mat*eph(6,:);

%Get mean anomalies at reference epoch

	M0 = vec2mat*eph(7,:);

%Get correction coefficients (argument of perigee)
	
	Cuc = vec2mat*eph(8,:);
	Cus = vec2mat*eph(10,:);

%Get correction coefficients (geocentric distance)
	
	Crc = vec2mat*eph(17,:);	
	Crs = vec2mat*eph(5,:);

%Get correction coefficients (inclination)
	
	Cic = vec2mat*eph(13,:);
	Cis = vec2mat*eph(15,:);

%Calculate anomalies

	M = M0 + ((GM2.*a.^(-1.5)+deltan).*dt);
	E = eccentric(M,e);
	f = 2*atan(sqrt((1+e)./(1-e)).*tan(E/2));

%Calculate the argument of latitude

	u = omega0 + f;

%Calculate the remaining Keplerian elements

	i = i0 + Cic.*cos(2*u) + Cis.*sin(2*u) + idot.*dt;
	omega = omega0 + Cuc.*cos(2*u) + Cus.*sin(2*u);
	l = OMEGA0 + OMEGAdot.*dt - omegaE.*dt0*ones(1,length(toe));


	
%Calculate range

	r = a.*(1-e.*cos(E))+ Crc.*cos(2*u) + Crs.*sin(2*u);

%Calculate clock polynomial coefficient

	clksv = vec2mat*eph(1,:) + vec2mat*eph(2,:).*dt + vec2mat*eph(3,:).*(dt.^2) ;

%Revert time-independent variables to vectors

	a=a(1,:);
	e=e(1,:);
