	program ephem
c
c	This program reads a set of orbital elements, then determines
c	predicted rho and theta values over a requested set of dates.
c
	common p,a,xi,xnode,t0,ecc,omega
	degrad=57.2957795
c
c	read star name and orb. if orb=yes then read orbital elements
c
	write(6,901) 
  901	format(' Enter orbital elements (P, a, i, W, T e, w): ')
	read(5,*) p,a,xi,xnode,t0,ecc,omega
	write(6,902) p,a,xi,xnode,t0,ecc,omega
  902	format(/' Orbital elements are as follows:'/
     $     '  P = ',f12.4,'  a = ',f12.4,
     $     '  i = ',f12.4,'  W = ',f12.4/
     $     '  T = ',f12.4,'  e = ',f12.4,
     $     '  w = ',f12.4,/)
	omega=omega/degrad
	xnode=xnode/degrad
	xi=xi/degrad
c
c	Get starting date and increment
c
	write(6,903)
  903	format(' Enter starting date and increment in years: ',$)
	read(5,*) bess0,dbess
c
	do 200 n=1,25
	bess=bess0+real(n-1)*dbess
	call porbit(bess)
  200	continue
c
  999	stop
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine porbit(date)
	common p,a,xi,xnode,t0,ecc,omega
	degrad=57.2957795
	ee=sqrt((1.+ecc)/(1.-ecc))
	aee=a*(1.-ecc*ecc)
	cosi=cos(xi)
	xmu=6.2831853/p
	em=xmu*(date-t0)
c
	e0=em+ecc*sin(em)+0.5*ecc*ecc*sin(2.*em)
	do 200 l=1,10
	em0=e0-ecc*sin(e0)
	e0=e0+(em-em0)/(1.-ecc*cos(e0))
  200	continue
c
	xnu=2.*atan(ee*tan(0.5*e0))
	r=aee/(1.+ecc*cos(xnu))
	theta=xnode+atan(cosi*tan(xnu+omega))
	rho=r*cos(xnu+omega)/cos(theta-xnode)
	theta=theta*degrad
	if(rho .gt. 0.) go to 300
	rho=-rho
	theta=theta+180.
  300	if(theta .lt.   0.) theta=theta+360.
	if(theta .gt. 360.) theta=theta-360.
c
c
	write(6,901) date,theta,rho
  901	format(1x,f9.4,2x,f10.2,f12.4)
  	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
