	program orb6cephem
c
c	This version of ephem reads several sets of orbital 
c	elements, then determines predicted rho and theta 
c	values over a requested set of dates. 
c
	character*1 pcode,acode,tcode
	character*8 ref,ele1,ele2,ele3
	character*10 star1
	character*12 star2
	common p,a,xi,xnode,t0,ecc,omega
	degrad=57.2957795
	date0 = 2006.
c
	open(10,file='orb6cephem.in',status='UNKNOWN')
	open(12,file='orb6cephem.out',status='UNKNOWN')
c
c	read star name and orbital elements.
c
  500	read(10,906,end=700) star1,star2,igrade,ref,p,pcode,a,acode,xi,
     $     xnode,t0,tcode,ecc,omega
  906	format(t20,a10,1x,a12,t151,i1,t189,a8,t61,f13.6,a1,f11.5,a1,
     $     2f10.4,3x,f12.6,a1,f10.6,f10.4)
c
        if (acode .eq. 'a') acode=' '
        if (pcode .eq. 'c') p=p*100.
        if (pcode .eq. 'd') p=p/365.2421987
	if (pcode .eq. 'h') p=p/365.2421987/24.
	if (pcode .eq. 'm') p=p/365.2421987/24./60.
        if (tcode .eq. 'd') t0=(t0-33282.423)/365.2421987 +1950.
	omega=omega/degrad
	xnode=xnode/degrad
	xi=xi/degrad
c
	ddate=0.05
	if (p .gt.  10.) ddate=0.1
	if (p .gt.  25.) ddate=0.25
	if (p .gt.  50.) ddate=0.5
	if (p .gt. 200.) ddate=1.0
c
	do 600 n=1,20
	date=date0+real(n-1)*ddate
	call porbit(date,theta,rho)
     	if (n .eq. 1) write(12,907) star1,star1,star2,
     $     igrade,ref,date,theta,rho,acode
  907	format('<a name = "',a10,'">',a10,'</a>',2x,a12,1x,
     $     i2,4x,a8,5x,f7.2,f8.1,f8.3,1x,a1)
	if (n .ne. 1) write(12,908) date,theta,rho,acode
  908	format(48x,f7.2,f8.1,f8.3,1x,a1)
  600	continue
c
	go to 500
c
  700	stop
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine porbit(date,theta,rho)
	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
  	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
