	program dephem
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. Header info 
c	and links are added to convert output to html file.
c
	character*23 star
	character*1 pcode,acode,tcode
	dimension theta1(1000),rho1(1000),bess(1000)
	common p,a,xi,xnode,t0,ecc,omega
	degrad=57.2957795
c
	open(10,file='dephem.in',status='UNKNOWN')
	open(11,file='dephem.out',status='UNKNOWN')
c
c	Give starting date and increment
c
	bess0 = 2015.06
	dbess = 0.0001
	nbess=100
	do 100 n=1,nbess
  	bess(n)=bess0+(n-1)*dbess
  100   continue
c
c	read first set of elements, derive ephemerides
c
   	read(10,902,end=999) star,p,pcode,a,acode,xi,xnode,t0,
     $     tcode,ecc,omega
  902	format(t20,a23,t81,f12.6,a1,t106,f9.5,a1,t126,f8.4,
     $     t144,f8.4,t163,f12.6,a1,t188,f8.6,t206,f8.4)
        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 (acode .eq. 'm') a=a/1000.
        if (tcode .eq. 'd') t0=(t0-33282.423)/365.2421987+1950.
        if (tcode .eq. 'm') t0=(t0+0.5-33282.423)/365.2421987+1950.
c
	omega=omega/degrad
	xnode=xnode/degrad
	xi=xi/degrad
c
	do 200 n=1,nbess
	call porbit(bess(n),theta1(n),rho1(n))
  200	continue
c
	do 400 n=1,nbess
	write(11,903) bess(n),theta1(n),rho1(n)
  903	format(1x,f10.4,f10.2,f11.5)
  400	continue
c
  999	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
