!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
! Copyright D.J. Jeffery, 2015jan01.
!
! A fortran-95 program to compute sunrise times
! for any latitude in the northern hemisphere.
! The compile command is f95 report.f.
! The run command is a.out. 
!
! To get accurate values go to 
!    http://aa.usno.navy.mil/data/docs/RS_OneYear.php
!
! For L=0:  errors of  =< 60 min.
! For L=10:  errors of  =< 40 min.
! For Las Vegas, L=36:  errors of  =< 30 min.
! For L=50:  errors of  =< 60 min.
! For L=70:  errors of  =< 120 min.
! For L=80:  errors of  =< 80 min.
!
!-----------------------------------------------------------------------
!
      program sunrise 
      implicit none
!
      integer :: i,j,k,l,m,n
      integer :: i_hour 
!
      real (kind=kind(0.d0)) :: coef 
      real (kind=kind(0.d0)) :: coef_f  ! The fitted coef.
      real (kind=kind(0.d0)) :: pi2 
      real (kind=kind(0.d0)) :: t_hour
      real (kind=kind(0.d0)) :: t_month 
      real (kind=kind(0.d0)) :: t_sunrise
      real (kind=kind(0.d0)) :: t_vernal
      real (kind=kind(0.d0)) :: xl_arctic
!      real (kind=kind(0.d0)) :: xlatitude=0.0d0
!      real (kind=kind(0.d0)) :: xlatitude=10.0d0
      real (kind=kind(0.d0)) :: xlatitude=36.0d0
!      real (kind=kind(0.d0)) :: xlatitude=50.0d0
!      real (kind=kind(0.d0)) :: xlatitude=70.0d0
!      real (kind=kind(0.d0)) :: xlatitude=80.0d0
!      real (kind=kind(0.d0)) :: xlatitude=90.0d0
      real (kind=kind(0.d0)) :: xminute
!
      t_vernal=78.75d0/(365.25d0/12.d0)
      coef=(24.d0/360.d0)*23.4d0
      pi2=2.*acos(-1.0d0)
      xl_arctic=(360.d0/pi2)*atan(coef/6.d0)
      coef_f=6.d0*tan((pi2/360.d0)*(23.4d0))
      print*,'t_vernal,coef,xl_arctic,coef_f'
      print*,t_vernal,coef,xl_arctic,coef_f
!   2.5872689938398357        1.5599999999999998        14.574216198038735        2.5964318534845550
      print*
!
      do i=1,12
        t_month=real(i,kind(0.d0))
        t_hour=t_sunrise(t_month,xlatitude) 
        i_hour=floor(t_hour)
        xminute=modulo(t_hour,1.d0)*60.d0
        write(*,'(i5,f6.2,f10.3,i5,f6.2)')                               &
     &    i,t_month,t_hour,i_hour,xminute
      end do
!
      end program sunrise 
!
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
!
      function t_sunrise(t_month,xlatitude)
!
      real (kind=kind(0.d0)) :: coef 
      real (kind=kind(0.d0)) :: coef_f 
      real (kind=kind(0.d0)) :: pi2 
      real (kind=kind(0.d0)) :: t_hour
      real (kind=kind(0.d0)) :: t_month
      real (kind=kind(0.d0)) :: t_sunrise
      real (kind=kind(0.d0)) :: xlatitude
!
      t_vernal=78.75d0/(365.25d0/12.d0)
      pi2=2.*acos(-1.0d0)
      coef_f=6.d0*tan((pi2/360.d0)*(23.4d0))
      coef=(24.d0/360.d0)*23.4d0
!      coef=coef_f
!
      t_sunrise=6.d0-coef*sin(pi2*(t_month-t_vernal)/12.d0)                   &
     &                   /tan(pi2*(90.d0-xlatitude)/360.d0)
!
!      print*,sin(pi2*(t_month-2.7d0)/12.d0)
!      print*,tan(pi2*(90.d0-xlatitude)/360.d0)
!      print*,t_month,t_hour
!      print*,coef  ! 1.5599999999999998
!
      return
      end function t_sunrise
!
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12