	program orb7ephem
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*1 junk,keep
	character*8 ref
	character*10 wds
	character*12 star
	character*17 ele
	character*130 header
	dimension theta(10),rho(10)
	integer ibess(10)
	common p,a,xi,xnode,t0,ecc,omega
	degrad=57.2957795
c
	open(10,file='orb6.master',status='UNKNOWN')
	open(11,file='/data/Orbits/Orb6/orb6ephem.header',status='UNKNOWN')
	open(12,file='/home/bdm/orb7ephem.html',status='UNKNOWN')
	open(13,file='/home/bdm/orb7.ephem',status='UNKNOWN')
	open(14,file='/home/bdm/orb7ephem.body',status='UNKNOWN')
	open(15,file='/home/bdm/orb7ephem.check',status='UNKNOWN')
c
c	Give starting date and increment
c
	bess0 = 2001.
	dbess = 1.0
	do 100 n=1,25
  100	ibess(n)=nint(bess0+(n-1)*dbess)
c
c	read header lines from orbit catalog
c
	do 200 n=1,4
  200	read(10,901) junk
  901	format(a1)
	norbits=0
c
c	write header lines for hmw5ephem.body
c
	write(14,902)
  902	format('<html>'/'<target="_top">'/'<body bgcolor="white"',
     $     ' text="black" link="blue" vlink="#0066CC">'/'<h4>'/'<pre>')
c
c	read and write header lines for html file
c
  300	read(11,903,end=400) header
  903	format(a130)
  	write(12,903) header
	go to 300
c
c	write header lines for ascii file
c
  400	write(13,904)
  904	format('Sixth Catalog of Orbits of Visual Binary Stars:',
     $     ' Ephemerides'//'WDS       Name         Reference',
     $     25('   Theta  Rho  '))
	write(13,905) (ibess(n),n=1,25)
  905	format(28x,25i15)
c
c	read star name and orbital elements. Loop back for another 
c	if orbit isn't a "keeper"
c
  500	read(10,906,end=700) wds,star,grade,keep,ref,p,a,xi,
     $     xnode,t0,ecc,omega,ele
  906	format(t23,a10,a12,t169,f3.1,t175,a1,1x,a8,t81,f12.6,
     $     f10.7,f9.5,f7.3,2x,f10.5,f9.6,f8.4,t149,a17)
	if (wds .eq. '          ') go to 500
	if (keep .ne. 'y') go to 500
	if (p .lt. 0.) go to 500
	if (a .lt. 0.) go to 500
	if (xi .lt. 0.) go to 500
	if (xnode .lt. 0.) go to 500
	if (t0 .lt. 0.) go to 500
	if (ecc .lt. 0.) go to 500
	if (omega .lt. 0.) go to 500
	if (ecc .gt. 1.0) write(6,997) wds
  997	format(a10)
	norbits=norbits+1
 	write(15,911) wds,star,p,a,xi,xnode,t0,ecc,omega
  911	format(1x,a10,a12,7f15.7)
	omega=omega/degrad
	xnode=xnode/degrad
	xi=xi/degrad
	igrade=nint(grade)
	if (grade .lt. 1.49) igrade=1
	if (grade .gt. 4.51) igrade=5
	if (grade .gt. 7.51) igrade=8
	if (grade .gt. 8.51) igrade=9
c
	do 600 n=1,25
	bess=bess0+real(n-1)*dbess
	call porbit(bess,theta(n),rho(n))
  600	continue
c
     	if (igrade .ne. 9) write(12,907) wds,wds,star,
     $     igrade,ref,ele,(theta(n),rho(n),n=1,25)
     	if (igrade .ne. 9) write(14,907) wds,wds,star,
     $     igrade,ref,ele,(theta(n),rho(n),n=1,25)
  907	format('<a name = "',a10,'">',a10,'</a>',2x,a12,1x,i2,2x,
     $     a8,2x,a17,1x,5(f8.1,f7.3))
	if (igrade .ne. 9) write(13,908) wds,star,ref,
     $     (theta(n),rho(n),n=1,25)
  908	format(a10,a12,1x,a8,1x,25(f8.1,f7.3))
     	if (igrade .eq. 9) write(12,912) wds,wds,star,
     $     igrade,ref,ele
     	if (igrade .eq. 9) write(14,912) wds,wds,star,
     $     igrade,ref,ele
  912	format('<a name = "',a10,'">',a10,'</a>',2x,a12,1x,i2,2x,
     $     a8,2x,a17,4x,'astrometric orbit')
	if (igrade .eq. 9) write(13,913) wds,star,ref
  913	format(a10,a12,1x,a8,4x,'astrometric orbit')

	go to 500
c
  700	write(12,909)
  909	format('</pre>')
     	write(14,910)
  910	format('</pre>'/'</html>')
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
