program ephem2 c (double precision version, with IAU expressions for c Besselian epochs) 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. c implicit double precision (a-h,o-z) character*23 star,str character*1 pcode,acode,tcode,rflag,dsign character*18 coords character*3 grade character*8 ref dimension theta(25),rho(25),xbess(25) integer ibess(25) common p,a,xi,xnode,t0,ecc,omega parameter ( pi = 3.14159265358979324D0 ) parameter ( degrad = 180.D0 / pi ) parameter ( b1900 = 15020.31352D0 ) parameter ( tropyr = 365.242198781D0 ) c open (10, file='ephem.in', status='OLD') open (11, file='ephem2.out', status='UNKNOWN') c c Give starting date and increment c write(6,991) 991 format('Enter starting date and time increment: ',$) read(5,*) bess0,dbess nbess=8 do 100 n=1,nbess xbess(n)=bess0+(n-1)*dbess ibess(n)=nint(xbess(n)) 100 continue write(11,901) (xbess(n),n=1,nbess) 901 format(20x,8f16.1,' * grade ref') c c read set of elements, derive ephmerides c 200 read(10,902,end=999) star,p,pcode,a,acode,xi,xnode,t0, $ tcode,ecc,omega,rah,ram,ras,dsign,decd,decm,decs, $ eqnx,grade,ref 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, $ t1,2f2.0,f5.1,a1,2f2.0,f4.1,t224,f4.0,t244, $ a3,5x,a8) 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 omega=omega/degrad xnode=xnode/degrad xi=xi/degrad c c compute ephemerides for beginning of Besselian years 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=' ' 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 if (rmax .lt. 999.99D0) go to 630 rflag='M' do 620 n=1,nbess 620 rho(n)=rho(n)/60.D0 go to 700 c 630 if (rmin .gt. 0.01D0) go to 700 rflag='m' do 640 n=1,nbess 640 rho(n)=rho(n)*1000.D0 c c write results to file c 700 write(11,903) star,(theta(n),rho(n),n=1,nbess), $ rflag,grade,ref 903 format(a23,1x,8(f8.1,f8.3),2x,a1,2x,a3,2x,a8) c go to 200 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 ) 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 c 900 write(11,901) date, em, e0, xnu, rho, theta c 901 format(' date, em, e0, xnu, rho, theta = ',6f12.5) c return end c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c 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 PARAMETER ( pi = 3.14159265358979324D0 ) PARAMETER ( degrad = 180.D0 / pi ) dt = date - eqnx c Aitken formula (with constant updated) c delpa = 0.005567D0 * dt * DSIN(ra) / DCOS(dec) 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 900 write(11,901) date, eqnx, ra, dec, delpa c 901 format(' date, eqnx, ra, dec, delpa = ',5f12.5) return end