%
% M-file to evaluate motion equations
clear all;
n = (2*pi)/43200; % Mean motion of satellite 12-hour period
%
% Loop over eccentricity
figure; clf
hold on;
for e = 0.001:0.05:0.99
	% Loop over one orbit
	i = 0;
	for t = 0:600:43199
		i = i + 1;
		M = n*t;   % Mean anomaly
		% Compute eccentric anomaly by iteration
		ek = M; eko = ek;
		for j = 1:100
			ek = M + e*sin(ek);
			dek = ek - eko;
			eko = ek;
		end
		if( abs(dek) > 1.e-6 )
			fprintf(1,'Kepler equation not solved accurately %f\n',dek)
		end

		v(i) = atan2(sqrt(1-e^2)*sin(ek)/(1-e*cos(ek)),(cos(ek)-e)/(1-e*cos(ek)));
		if( v(i) < 0 ) 
			v(i) = v(i) + 2*pi;
		end
%		fprintf(1,' Anomalies Mean %10f Eccentric %10f True %10f \n',M*180/pi,ek*180/pi,v(i)*180/pi)
		Ms(i) = M;
		eks(i) = ek;
		ts(i) = t;
	end
	plot(ts,v-Ms,'r');
	plot(ts,eks-Ms,'g');
end