        program Julian_Day
        implicit none
        external rad
        real*8 rad,pi,JD,T,L,M,Lo,cosu,Dt,HJD
        real*8 a,d,e
        real*8 HoursA,MinutesA,SecondsA,DegreesD,MinutesD,SecondsD
        integer i,N
        parameter (N=45)
        open (1,file="jd")
        open (2,file="HJD")
 7      format(f13.5)

       
        write(*,*) ' Give coordinate RA (hours minutes seconds)'
        read(*,*) HoursA,MinutesA,SecondsA
        write(*,*) 'Give DEC coordinates (degrees minutes seconds)'
        read(*,*) DegreesD,MinutesD,SecondsD

        a= 15.d0*(HoursA + MinutesA/60.d0 + SecondsA/3600.d0)
        d= DegreesD + MinutesD/60.d0 + SecondsD/3600.d0
        e=23.45d0               !Ecliptic: 23o 27'
        pi=acos(-1.d0)


        do i=1,N
          read(1,7) JD         

          T = (JD - 2451545.d0) / 36525.d0

          L = 280.460d0 + (36000.772d0)*T

          M = 357.528d0 + (35999.050d0)*T

          Lo = L + (1.915d0 - 0.0048d0*T)*sin(rad(M)) +
     &         (0.020d0*sin(rad(2.d0*M)))

          cosu = sin(rad(d))*sin(rad(Lo))*sin(rad(e)) +
     &           cos(rad(d))*cos(rad(a))*cos(rad(Lo)) +
     &           cos(rad(d))*sin(rad(a))*sin(rad(Lo))*cos(rad(e))

          Dt = 8.316d0*cosu/1440.d0

          HJD = JD - Dt

          write(2,7) HJD
       enddo
       
       close (1)
       close (2)
       stop
       end
        
        
        real*8 function rad(x)
        implicit none
        real*8 x,pi
        
        pi=acos(-1.d0)
        rad=x*pi/180.d0
        
        return
        end
