	program orb6ephem
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
c       2014/09/25: modified to deal with grade 7 correctly
c
c       2022/01/26: use orb6ephemx.f. This version of program
c       retained for historical purposes.
c 
	character*1 junk,keep,aap,aaa,aai,aan,aat,aae,aao,
     $              acode,pcode,tcode
	character*8 ref
	character*10 wds
	character*15 star
	character*131 header
        character*1 decs
	dimension theta(10),rho(10)
	integer ibess(10)
	common p,a,xi,xnode,t0,ecc,omega
	degrad=57.2957795
c
	open(10,file='../Catalog/orb6.master',
     $     status='UNKNOWN')
	open(11,file='orb6ephem.header',status='UNKNOWN')
	open(12,file='orb6ephem.html',status='UNKNOWN')
	open(13,file='orb6ephem.txt',status='UNKNOWN')
	open(14,file='orb6ephem.body',status='UNKNOWN')
	open(15,file='orb6ephem.check',status='UNKNOWN')
c
c	Give starting date and increment
c
        pm=0.0
        equ=2000.
c
	bess0 = 2021.
	dbess = 1.0
	do 100 n=1,5
  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 orb6ephem2.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(a131)
  	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 ',
     $     5('   Theta  Rho  '))
	write(13,905) (ibess(n),n=1,5)
  905	format(32x,5i15)
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,pcode,a,acode,
     $     xi,xnode,t0,tcode,ecc,omega,aap,aaa,aai,aan,aat,aae,aao,
     $     rah,ram,decs,decd,decm
  906	format(t20,a10,1x,a15,t244,f3.1,3x,a1,1x,a8,t81,f12.6,a1,
     $     t106,f9.5,a1,t126,f8.4,t144,f8.4,t163,f12.6,a1,t188,f8.6,
     $     t206,f8.4,t85,a1,t108,a1,t128,a1,t146,a1,t167,a1,t188,a1,
     $     t208,a1,t20,f2.0,f3.1,a1,2f2.0)
	if (wds .eq. '          ') go to 500
	if (keep .ne. 'y') go to 500
c
	igrade=nint(grade)
	if (grade .lt. 1.49) igrade=1
	if (grade .gt. 4.51) igrade=5
	if (grade .gt. 5.51) igrade=6
	if (grade .gt. 6.51) igrade=7
	if (grade .gt. 7.51) igrade=8
	if (grade .gt. 8.51) igrade=9
c
	if (igrade .eq. 7) go to 650
	if (aap .eq. ' ') go to 650
	if (aaa .eq. ' ') go to 650
	if (aai .eq. ' ') go to 650
	if (aan .eq. ' ') go to 650
	if (aat .eq. ' ') go to 650
	if (aae .eq. ' ') go to 650
	if (aao .eq. ' ') go to 650
	if (ecc .gt. 1.0) write(6,997) wds
  997	format(a10)
        if (acode .eq. 'm') a=a/1000.
        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.
        if (tcode .eq. 'm') t0=(t0+0.5-33282.423)/365.2421987 +1950.
c
	ra=rah*15. + ram/4.
	dec=decd + decm/60.
	if (decs .eq. '-') dec=-dec
	corr = (0.0056*sin(ra/degrad)/cos(dec/degrad))
     $     +0.00417*pm*sin(dec/degrad)
c
        norbits=norbits+1
 	write(15,911) wds,star,p,a,xi,xnode,t0,ecc,omega
  911	format(1x,a10,a15,7f15.7)
	omega=omega/degrad
	xnode=xnode/degrad
	xi=xi/degrad
c
	do 600 n=1,5
	bess=bess0+real(n-1)*dbess
	call porbit(bess,corr,theta(n),rho(n))
  600	continue
c
     	if (igrade .ne. 9) write(12,907) wds,wds,star,
     $     igrade,ref,(theta(n),rho(n),n=1,5)
     	if (igrade .ne. 9) write(14,907) wds,wds,star,
     $     igrade,ref,(theta(n),rho(n),n=1,5)
  907	format('<a name = "',a10,'">',a10,'</a>',2x,a15,1x,i2,6x,
     $     a8,2x,5(f8.1,f7.3))
c
	if (igrade .ne. 9) write(13,908) wds,star,ref,
     $     (theta(n),rho(n),n=1,5)
  908	format(a10,1x,a15,1x,a8,1x,5(f8.1,f7.3))
     	if (igrade .eq. 9) write(12,912) wds,wds,star,igrade,ref,
     $     (theta(n),rho(n),n=1,5)
c
     	if (igrade .eq. 9) write(14,912) wds,wds,star,igrade,ref,
     $     (theta(n),rho(n),n=1,5)
  912	format('<a name = "',a10,'">',a10,'</a>',2x,a15,1x,i2,6x,
     $     a8,2x,5(f8.1,f7.3),'  astrometric orbit')
	if (igrade .eq. 9) write(13,913) wds,star,ref,
     $     (theta(n),rho(n),n=1,5)
  913	format(a10,1x,a15,1x,a8,1x,5(f8.1,f7.3),'  astrometric orbit')
	go to 500
c
  650	write(12,914) wds,wds,star,igrade,ref
  914	format('<a name = "',a10,'">',a10,'</a>',2x,a15,1x,i2,6x,
     $     a8,2x,5('      .    .   '),'  incomplete elements')
     	write(13,915) wds,star,ref
  915	format(a10,1x,a15,1x,a8,4x,'incomplete elements')
  	write(14,914) wds,wds,star,igrade,ref
	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,corr,theta,rho)
	common p,a,xi,xnode,t0,ecc,omega
	degrad=57.2957795
        equ=2000.
	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 + corr*(date-equ)
	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
