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. Header info c and links are added to convert output to html file. c character*1 junk,keep,aap,aaa,aai,aan,aat,aae,aao, $ acode,pcode,tcode character*8 ref character*10 wds,wdsc(100) character*15 star,starc(100) character*1 decs character*80 header common p,a,xi,xnode,t0,ecc,omega degrad=57.2957795 date0 = 2013 c open(10,file='orb6c.stars',status='UNKNOWN') open(11,file='../Catalog/orb6.master',status='UNKNOWN') open(12,file='orb6cephem.out',status='UNKNOWN') open(13,file='orb6cephem.body',status='UNKNOWN') open(14,file='orb6cephem.top',status='UNKNOWN') open(15,file='orb6cephem.html',status='UNKNOWN') c c write header info for various files c write(13,913) 913 format(''/''/''/'

'/'
')
	do 10 n=1,9
	read(14,914) header
  914	format(a80)
	write(15,914) header
   10	continue
c
c	read in calibration stars
c
	ncal=1
   50	read(10,901,end=80) wdsc(ncal),starc(ncal)
  901	format(a10,1x,a15)
	ncal=ncal+1
	go to 50
c
   80   ncal=ncal-1
	pm=0.0
        equ=2000.
c
c	read header lines from orbit catalog
c
	do 200 n=1,4
  200	read(11,902) junk
  902	format(a1)
	norbits=0
c
c	read star name and orbital elements. Loop back for another 
c	if orbit isn't a "keeper"
c
  500	read(11,903,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
  903	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
	iflag=0
	do 550 nc=1,ncal
	if (wds .ne. wdsc(nc)) go to 550
	if (star .ne. starc(nc)) go to 550
	iflag=1
	go to 570
  550	continue
  570	if (iflag .eq. 0) go to 500
c
	if (aap .eq. ' ') go to 500
	if (aaa .eq. ' ') go to 500
	if (aai .eq. ' ') go to 500
	if (aan .eq. ' ') go to 500
	if (aat .eq. ' ') go to 500
	if (aae .eq. ' ') go to 500
	if (aao .eq. ' ') go to 500
c
        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.
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
	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. 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
	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,904) wds,wds,star,
     $     igrade,ref,date,theta,rho,acode
     	if (n .eq. 1) write(13,904) wds,wds,star,
     $     igrade,ref,date,theta,rho,acode
     	if (n .eq. 1) write(15,904) wds,wds,star,
     $     igrade,ref,date,theta,rho,acode
  904	format('',a10,'',2x,a15,1x,
     $     i2,4x,a8,5x,f7.2,f8.1,f8.3,1x,a1)
	if (n .ne. 1) write(12,905) date,theta,rho,acode
	if (n .ne. 1) write(13,905) date,theta,rho,acode
	if (n .ne. 1) write(15,905) date,theta,rho,acode
  905	format(47x,f7.2,f8.1,f8.3,1x,a1)
  600	continue
c
	go to 500
c
  700	write(13,906)
  906	format('
'/'') do 710 n=1,3 read(14,914) header write(15,914) header 710 continue c 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