	subroutine jayday(ny,nm,nd,jnl)
	integer ny,nm,nd,my,mm
        real jnl
	my=ny
	mm=nm
	if (mm .gt. 2 ) go to 12
	mm=mm + 9
	my=my -1
	go to 13
12	mm=mm-3
13      jnl =15078.0 + aint(1461.0*float(my)/4) 
        jnl = jnl + aint((153.0*float(mm)+ 2.0)/5.0) + float(nd)
	return
	end
	subroutine calday(jnl,ny,nm,nd)
        real jnl
	if ( jnl .lt. 15384.0 ) go to 10
	if (jnl .lt. 88067.0) go to 11
10	nm=0
	nd=0
	ny=0
	goto 13
11	nseq=jnl-15078
	ny=(4*nseq-1)/1461
	nd=(4*nseq + 3 -1461*ny)/4
	nm= (5*nd-3)/153
	nd=(5*nd+2-153*nm)/5
	if ( nm .ge. 10) goto 12
	nm=nm +3
	goto 13
12	nm=nm-9
	ny=ny+1
13	return
	end
	function sidtim(jday,fday)
	double precision const0,const1,dthdt,a
        real jday
	dthdt=1.00273790929
	twopi=6.283185308
	const0=.277987616
	const1 = .27379092965e-2
	a= const0 + const1*dble(jday -33282.0)
	a =a + dthdt*dble(fday)
	sidtim= amod(sngl(a),1.0)*twopi
	return
	end
	subroutine eapos(x,y,z,t)
	real w,p,ec
	real *8 mdt,m,twopi,dtu
C     
C ------Co-ordinates of earth in Ecliptic reference frame of date
C ------ only approximate ( i.e. longitude accuracy +/- 1 minute of Arc
C ------ Reference Almanac for computers -1980
	data mdt,twopi/6.283019465D+2,6.28318531D+0/
	dtu = dble(t-15019.5)/3.6525D+4
	m= dble(-.026601727) + dtu*mdt
	m = dmod(m,twopi)
	w = sngl(m)
	w = w +.0334409*sin(w) + 3.4907E-4*sin(2.0*w)
	ec =.01675104 -.418E-4*sngl(dtu)
	p= 1.0 -ec*ec
	om = 1.766636813 + .030005*sngl(dtu)
	xw =p*cos(w)/(1.0 + ec*cos(w))
	yw = p*sin(w)/(1.0 + ec*cos(w))
	co = cos(om)
	so = sin(om)
	x = xw*co -yw*so
	y = xw*so + yw*co
	z = 0
	return
	end
	function dot(v1,v2)
	dimension v1(3),v2(3)
C -----  Vector dot product
	dot = 0.0
	do 1 i =1,3
1	dot = dot + v1(i)*v2(i)
	return
	end
	subroutine cross(v1,v2,v3)
C  
C ----- Cross product of two vectors
	dimension v1(3),v2(3),v3(3)
	v3(1) = v1(2)*v2(3) -v1(3)*v2(2)
	v3(2) =v1(3)*v2(1) - v1(1)*v2(3)
	v3(3) = v1(1)*v2(2) -v1(2)*v2(1)
	return
	end
 	function decml(ch,n)
C	DECML .. function for BCD to floating point binary conversion
C	Nicole Simon
C 
C	Revised for FORTRAN-80 by Anthony BERESFORD 30 DEC 1979
	logical*1 ch(1),bcd(13)
	data bcd/1h1,1h2,1h3,1h4,1h5,1h6,1h7,1h8,1h9,1h-,1h.,1he/
	data bcd(13)/1hE/
C   
 	L=lstrng(ch,n)
6	in=n
7	kexpon = 0
8	kpoint = 0
C	convert bcd to fixed point number
10	decml= 0
20	fmlplr = 10.**l
25	last = n+l -1
30	do 110 k= in,last
50	do 90 intger= 1,9
60	IF (ch(k) .ne. bcd(intger))goto 90
	fmlplr = fmlplr/10.
	decml= decml + float(intger)*fmlplr
	goto 110
90	continue
	if (ch(k) .eq. bcd(10))goto 100
	if (ch(k) .eq. bcd(11))goto 105
	if(ch(k) .eq. bcd(12) .or. ch(k) .eq. bcd(13))goto 115
	goto 101
100	fmlplr = -fmlplr
101	fmlplr = fmlplr/10.
	goto 110
105	kpoint = k
110	continue
C   
111	goto 120
C        
115	kexpon = k
	in = k + 1
	il = last -kexpon
120	if ( kpoint .ne. 0 ) goto 140
	if (kpoint .eq. 0 .and. kexpon .eq. 0)goto 210
	kpoint =kexpon
140	dvisor = 10.**(last-kpoint+1)
	decml =decml/dvisor
C       
	if ( kexpon .eq. 0)goto 210
C	convert to floating point number,shifting for exponent
	kpoint = kdecml(ch,in,il)
	decml =decml*10.**kpoint
C       
210	n = n + l
220	return
	end
	function idecml(ch,n)
	logical*1 ch(1)
C        IDECML  .... function for bcd to fixed point binary conversion
C        Nicole Simon
	data limit/6/
C	limit = maximum no. of digits allowed for integers
5	l =lstrng(ch,n)
	if (l .lt. limit )goto 10
	idecml =kdecml(ch,n,limit)
	goto 220
10	idecml = kdecml(ch,n,l)
220	n= n + l
   	return
	end
	function kdecml(ch,n,l)
	logical*1 ch(1),bcd(10)
	real mulplr
	data bcd/1h1,1h2,1h3,1h4,1h5,1h6,1h7,1h8,1h9,1h-/
C        
	kdecml = 0
	mulplr = 10.**l
	last = n+ l -1
	do 110 k =n,last
 	mulplr =mulplr/10.
	do 90 intger = 1,9
	if (ch(k) .ne. bcd(intger))goto 90
	kdecml = kdecml + ifix(float(intger)*mulplr)
	goto 110
90	continue
	if (ch(k) .ne. bcd(10) ) goto 110
	mulplr = - mulplr
110	continue
	return
	end
	function lstrng(ch,n)
C       LSTRNG  ... function to search for length of input string
C       starting in column  N
	logical *1 ch(1),blank
	data last,label,blank/80,81,1h /
	do 50 k = n,last
	if ( ch(k) .ne. blank )goto 50
	lstrng = k -n
	if( lstrng .ne. 0) goto 70
	n = n +1
50	continue
	lstrng = label -n
70	return
	end
	subroutine eaxyz(x,y,z,t,day)
C * * * to give Earths's co-ordinates in Cartesian form
C * * * for the reference frame 1950.0
	real*4 w,ec,lam,l1950,kep,ld1950
	real*8 mdt,m,twopi,dtu
	data mdt,twopi/6.28301945717D+2,6.2831853071795864D+0/
	data t1950,prec/33282.0,2.435637E-4/
	data obli0,obli1,ecrot/.4093197,-2.27111E-4,2.28105E-6/
	dtu =(dble(t-15019.) +dble(day-0.5))/3.6525D+4
	ec = .01675104 -.418E-4*sngl(dtu)
	m= dble(-.026601727) + dtu*mdt
	m=dmod(m,twopi)
	w = sngl(m)
	ecan = kep(w,ec)
	se = sin(ecan)
	ce = cos(ecan)
	r= 1.0 -ec*ce
	upsil = atan2(sqrt(1.-ec*ec)*se,ce-ec)
	om = 1.766636813 + .030005*sngl(dtu)
	dy = (t-t1950)/365.25
	lam = om + upsil
	l1950 = lam - prec*dy
	beta =-ecrot*dy*sin(lam +.091367)
	xe = r*cos(l1950)*cos(beta)
	ye =r*sin(l1950)*cos(beta)
	ze =r*sin(beta)
	obli =obli0 +obli1*sngl(dtu)
	x= xe
	y =ye*cos(obli) -ze*sin(obli)
	z= ye*sin(obli) +ze*cos(obli)
	return
	end
	real function Kep(am,ecc)
	data eps/1.0E-8/
	eni = am + ecc*sin(am)
10	de =(eni -ecc*sin(eni)-am)/(1.0 -ecc*cos(eni))
	eni = eni -de
	if (abs(de) .gt. eps)goto 10
  	kep = eni
	return
	end
