c  sunah.f
c  Azimuth and elevation of the Sun for the given location and dates.
c  Azimut a vyska Slunce pro dane stanoviste a datumy.
c  Miroslav Broz (miroslav.broz@email.cz), Jun 28th 2004

      program SUN_AZIMUTH_HEIGHT

c-----------------------------------------------------------------------

      implicit none

c  constants
      real*8 pi, deg, rad
      parameter(pi = 3.1415926535d0, deg = pi/180.d0, rad = 180.d0/pi)

c  input parameters
      character*80 phistr, lambdastr, datumstr, timestr, dtstr
      integer nsteps, precession, nutation, refraction

c  internal variables
      integer i, iu
      real*8 jdate, year, month, day, hour, minute, second
      real*8 phi, lambda, jd0, dt, a, h, ra, de, t, l, b
      character*80 s

c  functions
      real*8 frac, reads3, jd

c-----------------------------------------------------------------------
c
c  read stdin / cti standardni vstup
c
      iu=5      ! stderr - this does not work using astrohk/g77?!
      i=0
      write(*,10) '# zemepisna sirka stanoviste [+DD:MM:SS]: ' 
5     continue
        i=i+1
        read(*,10,end=15,err=15) phistr
10      format(a)
      if (phistr(1:1).eq.'#') goto 5

      i=i+1
      write(*,10) '# zemepisna delka [+DD:MM:SS]: '
      read(*,10,end=15,err=15) lambdastr
      i=i+1
      write(*,10) '# datum [YYYY/MM/DD]: '
      read(*,10,end=15,err=15) datumstr
      i=i+1
      write(*,10) '# pocatecni svetovy cas [HH:MM:SS]: '
      read(*,10,end=15,err=15) timestr
      i=i+1
      write(*,10) '# casovy krok [HH:MM:SS]: '
      read(*,10,end=15,err=15) dtstr
      i=i+1
      write(*,10) '# pocet kroku: '
      read(*,*,end=15,err=15) nsteps
      i=i+1
      write(*,10) '# pocitat s precesi? [0|1] '
      read(*,*,end=15,err=15) precession
      i=i+1
      write(*,10) '# pocitat s nutaci? [0|1] '
      read(*,*,end=15,err=15) nutation
      i=i+1
      write(*,10) '# zapocitat refrakci? [0|1] '
      read(*,*,end=15,err=15) refraction

      goto 17
15    continue
        write(*,16) i
16      format('# SUNAH: Chyba pri cteni vstupnich parametru na radku ',
     :    i2,'.')
        goto 990
17    continue

      phi = reads3(phistr)
      lambda = reads3(lambdastr)
      call reads_split3(datumstr, '/', year, month, day)
      call reads_split3(timestr, ':', hour, minute, second)
      if (month.gt.0) then
        day = day + (hour + minute/60.d0 + second/3600.d0) / 24.d0
        jd0 = jd(year, month, day)
      else
        jd0 = year
      endif
      dt = reads3(dtstr)/pi*180.d0/24.d0

      write(*,20) phi*rad
      write(*,21) lambda*rad
      write(*,22) jd0
      write(*,23) dt, dt*24.d0, dt*1440.d0, dt*86400.d0
      write(*,24) nsteps
      write(*,25) precession
      write(*,26) nutation
      write(*,27) refraction
20    format('# phi = ', f9.5, ' deg')
21    format('# lambda = ', f9.5, ' deg')
22    format('# jd0 = ', f13.5)
23    format('# dt = ', f8.5, ' day = ', f7.3, ' hour = ',
     :  f8.2, ' min = ', f7.0, ' sec')
24    format('# nsteps = ', i4)
25    format('# precession = ', i1)
26    format('# nutation = ', i1)
27    format('# refraction = ', i1)

      write(*,40)
40    format('# JD (UT)     YYYY MM DD HH MM SS',
     :  '   AZ        HT        HeLong    HA        DE        RA')

c-----------------------------------------------------------------------
c
c  calculate ephemeris and write stdout | spocti efemeridu a tiskni standardni vystup
c
      do i = 0, nsteps - 1
        jdate = jd0 + i * dt
        call sunah(jdate, lambda, phi, precession, nutation, refraction,
     :    a, h, ra, de, t, l, b)

c  round to 0.5 sec!
        call od(jdate + 0.5d0/86400.d0, year, month, day)
        call hhms(frac(day)*24.0d0, hour, minute, second)

        write(*,30) jdate, int(year), int(month), int(day),
     :    int(hour), int(minute), int(second),
     :    a*rad, h*rad, l*rad, t*rad, de*rad, ra*rad
30      format(f13.5,1x,i4,1x,i2,1x,i2,1x,
     :    i2,1x,i2,1x,i2,1x,
     :    f10.5,1x,f9.5,4(1x,f9.5))
      enddo

990   continue

      write(*,10) '# stiskni Enter pro ukonceni'
      read(*,10) s

      stop
      end

c-----------------------------------------------------------------------