program temp2 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 2015/02/26: modified extensively by George Kaplan. Converted to c double precision, with IAU expressions for Besselian epochs) c 2020/07/13: minor mod to make ephemeris 2019-2023 rather c than 2016-2020. You must change to the new starting c date on line 44 below, recompile the program, and make c the modification by hand in orb6ephem.txt, c orb6ephem.top, orb6ephem.html and orb6ephem.header. c implicit double precision (a-h,o-z) character*1 junk,keep,aap,aaa,aai,aan,aat,aae,aao, $ acode,pcode,tcode,rflag character*1 dsign character*8 ref character*10 wds character*15 star character*17 note character*138 header dimension theta(25),rho(25),xbess(25) integer ibess(25) common p,a,xi,xnode,t0,ecc,omega c parameter ( pi = 3.14159265358979324D0 ) parameter ( degrad = 180.D0 / pi ) parameter ( b1900 = 15020.31352D0 ) parameter ( tropyr = 365.242198781D0 ) 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') open(16,file='temp2.txt',status='UNKNOWN') c c Give starting date and increment c bess0 = 2015.D0 dbess = 0.5D0 nbess=5 do 100 n=1,nbess xbess(n)=bess0+(n-1)*dbess ibess(n)=nint(xbess(n)) 100 continue 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(''/''/''/'

'/'

')
c
c	read and write header lines for html file
c
  300	read(11,903,end=400) header
  903	format(a138)
  	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            Grade  ',
     $     'Reference',5('   Theta   Rho   '),'  Notes')
	write(13,905) (xbess(n),n=1,5)
  905	format(39x,5f17.1)
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=800) wds,star,grade,keep,ref,p,pcode,a,acode,
     $     xi,xnode,t0,tcode,ecc,omega,aap,aaa,aai,aan,aat,aae,aao,
     $     rah,ram,ras,dsign,decd,decm,decs,eqnx
  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,t1,2f2.0,f5.1,a1,2f2.0,f4.1,t224,f4.0)
  	write(16,906) wds,star,grade,keep,ref,p,pcode,a,acode,
     $     xi,xnode,t0,tcode,ecc,omega,aap,aaa,aai,aan,aat,aae,aao,
     $     rah,ram,ras,dsign,decd,decm,decs,eqnx
	if (wds .eq. '          ') go to 500
	if (keep .ne. 'y') go to 500
c
c        write(6,991) wds,star,grade,p,pcode,a,acode,xi,xnode,t0,tcode,
c     $     ecc,omega,rah,ram,ras,dsign,decd,decm,decs,eqnx
c  991   format(a10,1x,a15,f5.1,f13.6,1x,a1,f10.5,1x,a1,2f9.4,f13.6,1x,
c     $     a1,f9.6,f9.4/2f4.0,f6.1,1x,a1,f3.0,f4.0,f5.1,2x,f5.0)
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 750
	if (aap .eq. ' ') go to 750
	if (aaa .eq. ' ') go to 750
	if (aai .eq. ' ') go to 750
	if (aan .eq. ' ') go to 750
	if (aat .eq. ' ') go to 750
	if (aae .eq. ' ') go to 750
	if (aao .eq. ' ') go to 750
	if (ecc .gt. 1.0) write(6,997) wds
  997	format(a10)
c
	if (acode .eq. 'm') a=a/1000.D0 
	if (pcode .eq. 'c') p=p*100.D0
	if (pcode .eq. 'd') p=p/tropyr
	if (pcode .eq. 'h') p=p/tropyr/24.D0
	if (pcode .eq. 'm') p=p/tropyr/24.D0/60.D0
	if (tcode .eq. 'd') t0=(t0-b1900)/tropyr+1900.D0
	if (tcode .eq. 'm') t0=(t0-b1900+0.5D0)/tropyr+1900.D0 
	if (eqnx .le. 0.D0 ) eqnx=2000.D0
c
	ra = rah + ram/60.D0 + ras/3600.D0
	dec = decd + decm/60.D0 + decs/3600.D0
	if (dsign.eq.'-') dec=-dec 
c
	ra = ra*15.D0/degrad
	dec = 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
c       compute ephemerides for beginning of Besselian years
c
	do 600 n=1,nbess
	bess=bess0+((n-1)*dbess)
	call porbit(bess,pa,rho(n))
	call prec (bess,eqnx,ra,dec,delpa)
	pa=pa+delpa
	if (pa .lt.   0.D0) pa=pa+360.D0
	if (pa .gt. 360.D0) pa=pa-360.D0
	theta(n)=pa
  600   continue
c
c	rescale rho if too large or small
c
	rflag='a'
	rmax=rho(1)
	rmin=rho(1)
	do 610 n=2,nbess
	if (rho(n) .gt. rmax) rmax=rho(n)
	if (rho(n) .lt. rmin) rmin=rho(n)
  610	continue
c
  630	if (rmin .gt. 0.01D0) go to 700
	rflag='m'
c
c	write results to various web files
c
  700	note=' '
        if (igrade .eq. 9) note='astrometric orbit'
c
        if (rflag .eq. 'a') write(12,907) wds,wds,star,
     $     igrade,ref,(theta(n),rho(n),n=1,5),note
	if (rflag .eq. 'a') write(13,908) wds,star,igrade,ref,
     $     (theta(n),rho(n),n=1,5),note
     	if (rflag .eq. 'a') write(14,907) wds,wds,star,
     $     igrade,ref,(theta(n),rho(n),n=1,5),note
  907	format('',a10,'',2x,a15,1x,i2,6x,
     $     a8,2x,5(f8.1,f8.3,1x),2x,a17)
  908	format(a10,1x,a15,i4,4x,a8,1x,5(f8.1,f8.3,1x),2x,a17)

     	if (rflag .eq. 'm') write(12,912) wds,wds,star,igrade,ref,
     $     (theta(n),rho(n),n=1,5),note
	if (rflag .eq. 'm') write(13,913) wds,star,igrade,ref,
     $     (theta(n),rho(n),n=1,5),note
     	if (rflag .eq. 'm') write(14,912) wds,wds,star,igrade,ref,
     $     (theta(n),rho(n),n=1,5),note
  912	format('',a10,'',2x,a15,1x,i2,6x,
     $     a8,2x,5(f8.1,f9.4),2x,a17)
  913	format(a10,1x,a15,i4,4x,a8,1x,5(f8.1,f9.4),2x,a17)
	go to 500
c
c       incomplete elements
c
  750	write(12,914) wds,wds,star,igrade,ref
     	write(13,915) wds,star,igrade,ref
  	write(14,914) wds,wds,star,igrade,ref
  914	format('',a10,'',2x,a15,1x,i2,6x,
     $     a8,2x,5('      .     .    '),'  incomplete elements')
  915	format(a10,1x,a15,i4,4x,a8,1x,5('      .     .    '),
     $     '  incomplete elements')
	go to 500
c
  800	write(12,909)
  909	format('

') write(14,910) 910 format(''/'') c 999 stop end c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine porbit(date,theta,rho) implicit double precision (a-h,o-z) common p,a,xi,xnode,t0,ecc,omega parameter ( pi = 3.14159265358979324D0 ) parameter ( twopi = pi + pi ) parameter ( degrad = 180.D0 / pi ) c ee=dsqrt((1.D0+ecc)/(1.D0-ecc)) aee=a*(1.D0-ecc*ecc) cosi=dcos(xi) xmu=twopi/p em=xmu*(date-t0) c e0=em+ecc*dsin(em)+0.5D0*ecc*ecc*dsin(2.D0*em) do 200 l=1,10 em0=e0-ecc*dsin(e0) e0=e0+(em-em0)/(1.D0-ecc*dcos(e0)) 200 continue c xnu=2.D0*datan(ee*dtan(0.5D0*e0)) r=aee/(1.D0+ecc*dcos(xnu)) theta=xnode+datan(cosi*dtan(xnu+omega)) rho=r*dcos(xnu+omega)/dcos(theta-xnode) theta=theta*degrad if(rho .gt. 0.D0) go to 300 rho=-rho theta=theta+180.D0 300 if (theta .lt. 0.D0) theta=theta+360.D0 if (theta .gt. 360.D0) theta=theta-360.D0 c return end c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c subroutine prec (date,eqnx,coords, delpa) subroutine prec (date,eqnx,ra,dec, delpa) c c This subroutine computes the change in the direction of the c celestial pole due to precession, which is a correction to c be applied to the position angle. c IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*1 dsign CHARACTER*18 coords c PARAMETER ( pi = 3.14159265358979324D0 ) PARAMETER ( degrad = 180.D0 / pi ) c dt = date - eqnx c c Extract RA and Dec from character coordinate field c READ ( coords, '(2f2.0,f5.2,a1,2f2.0,f4.1)' ) rah,ram,ras, c $ dsign,decd,decm,decs c ra = rah + ram/60.D0 + ras/3600.D0 c dec = decd + decm/60.D0 + decs/3600.D0 c if (dsign.eq.'-') dec=-dec c c ra = ra*15.D0/degrad c dec = dec/degrad c c Aitken formula (with constant updated) c delpa = 0.005567D0 * dt * DSIN(ra) / DCOS(dec) c c Green formula as reported by Argyle (2004) z = ( 0.00640616D0 * dt + 3.041D-8 * dt*dt ) /degrad theta = ( 0.00556753D0 * dt + 1.185D-8 * dt*dt ) /degrad deltan = DSIN(ra-z)*DSIN(theta)/ $ ( DCOS(dec)*DCOS(theta)+DSIN(dec)*DSIN(theta)*DCOS(ra-z) ) delpa = DATAN ( deltan ) * degrad c c 900 write(11,901) date, eqnx, ra, dec, delpa c 901 format(' date, eqnx, ra, dec, delpa = ',5f12.5) c return end