program lunar_time

      integer yr,mon,day,sec
      double precision jd,t
	real lun_lt(61),lun,nu,lt

C Test to see if full moon on April 6 2012  (nu=180deg is full moon!)
        yr = 2012
        mon = 4
        day = 6
        sec = 0
	lt = 0.
        call ymd2jd(yr,mon,day,sec,jd)
        t = jd - 2415020.
        t = t / 36525. ! convert to julian centuries
        nu = -9.26009 + 445267.12165*t+0.00168*t*t
        nu = mod(nu,360.)
        lun = lt - nu*24./360.
        lun = mod(lun,24.)+24.
        write (6,"(1x,'yr mon day sec jd t nu lun lt =',4i6,d17.8,
     |    d17.8,2e13.4,f6.1)") yr,mon,day,sec,jd,t,nu,lun,lt

C 07325=11/21, 07335=12/01
C  calc for noon at Jicamarca (-11.98,-76.87)
        lt = 12.
        sec = (lt + 76.87/15.)*3600.
        yr = 2007.
        mon = 11.
        nd = 0
        do 1000 day=21,30
        nd = nd + 1
        call ymd2jd(yr,mon,day,sec,jd)
        t = jd - 2415020.
        t = t / 36525. ! convert to julian centuries
        nu = -9.26009 + 445267.12165*t+0.00168*t*t
        nu = mod(nu,360.)
        lun = lt - nu*24./360.
        lun_lt(nd) = mod(lun,24.)+24.
 1000   continue
        write (6,"(1x,'nd yr mon day sec jd t nu lun lt =',5i6,d17.8,
     |    d17.8,2e13.4,f6.1)") nd,yr,mon,day,sec,jd,t,nu,lun_lt(nd),lt
        mon = 12.	
        do 1001 day=1,31
        nd = nd + 1
        call ymd2jd(yr,mon,day,sec,jd)
        t = jd - 2415020.
        t = t / 36525. ! convert to julian centuries
        nu = -9.26009 + 445267.12165*t+0.00168*t*t
        nu = mod(nu,360.)
        lun = lt - nu*24./360.
        lun_lt(nd) = mod(lun,24.)+24.
 1001   continue
        write (6,"(1x,'nd yr mon day sec jd t nu lun lt =',5i6,d17.8,
     |    d17.8,2e13.4,f6.1)") nd,yr,mon,day,sec,jd,t,nu,lun_lt(nd),lt
        yr = 2008
        mon = 1.	
        do 1002 day=1,20
        nd = nd + 1
        call ymd2jd(yr,mon,day,sec,jd)
        t = jd - 2415020.
        t = t / 36525. ! convert to julian centuries
        nu = -9.26009 + 445267.12165*t+0.00168*t*t
        nu = mod(nu,360.)
        lun = lt - nu*24./360.
        lun_lt(nd) = mod(lun,24.)+24.
        write (6,"(1x,'nd yr mon day sec jd t nu lun lt =',5i6,d17.8,
     |    d17.8,2e13.4,f6.1)") nd,yr,mon,day,sec,jd,t,nu,lun_lt(nd),lt
 1002   continue
	write (6,"(1x,'lun =',10f6.1)") lun_lt
        stop
        end

      subroutine ymd2jd(yr,mon,day,sec,jd)
      implicit none
      integer itimes(3),j,yr,mon,day,sec
      double precision jd,a,y,m

      if (yr .le. 1000) then
        if (yr .le. 40) then
          yr = yr + 2000
        else
          yr = yr + 1900
        endif
      endif

      itimes(1) = int(mon)
      itimes(2) = int(day)
      itimes(3) = int(yr)

      a = INT((14.-dble(itimes(1)))/12.)
      y = itimes(3)+4800-a
      m = itimes(1)+12*a-3
      j = itimes(2) + INT((153.*dble(m)+2.) / 5.) + 365*y +
     .   INT(dble(y)/4.) - INT(dble(y)/100.) + INT(dble(y)/400.)- 32045
      jd = j + dble(sec)/86400.

      end