program testphaz implicit none save c..exercises the lunar phase routine c..declare character*22 comment character*13 what(4) character*3 amon(12) integer nmax,nfound parameter (nmax = 100) integer iphase(nmax) integer iyear,i,j,iy,im,id,iy1,im1,id1 double precision tjd(nmax),bet(nmax), 1 rhour,xh,xm,xs,rhour1,xh1,xm1,xs1 c..for picking which ephemeris to use integer which_eph common /dop/ which_eph c..initialize some names data what /'new moon ','first quarter', 1 'full moon ','third quarter'/ data amon /'jan' , 'feb' , 'mar' , 'apr' , 'may' , 'jun' , 1 'jul' , 'aug' , 'sep' , 'oct' , 'nov' , 'dec' / c..popular format statements 01 format(1x,t5,a,t21,a,t29,a,t32,a,t36,a,t47,a,t60,a) 02 format(1x,i2.2,' ',a,' ',i2.2,a,i4.4, 1 ' ',i2.2,':',i3.2,':',i2.2,' UT ',0pf5.1,' ',a) c..initialize the name database which_eph = 406 call init_names(0) c..get the input vector 100 write(6,*) 'give the year =>' read(5,*) iyear call moon_phases(iyear,tjd,iphase,bet,nmax,nfound) c..write the results write(6,*) write(6,01) 'event','date','hr','mn','sec','beta','comment' do i=1,nfound call caldat(tjd(i),iy,im,id,rhour) call dh2h(rhour,xh,xm,xs) c..find new moons, potential solar eclipse candidates c if (iphase(i) .eq. 1) then c if (iphase(i) .eq. 3) then c if (iphase(i) .eq. 1 .or. iphase(i) .eq. 3) then c..uncommenting the code below will find all blue moons within the year, c..a blue moon being the second full moon in a calander month c if (iphase(i) .eq. 3) then c j = i + 4 c call caldat(tjd(j),iy1,im1,id1,rhour1) c call dh2h(rhour1,xh1,xm1,xs1) c if (im .eq. im1) then comment = ' ' if (iphase(i) .eq. 1 .and. abs(bet(i)) .le. 1.0) 1 comment = 'total solar eclipse?' if (iphase(i) .eq. 1 .and. 1 abs(bet(i)) .gt. 1.0 .and. abs(bet(i)) .le. 1.5) 2 comment = 'partial solar eclipse?' if (iphase(i) .eq. 3 .and. abs(bet(i)) .le. 1.0) 1 comment = 'lunar eclipse?' write(6,02) i,what(iphase(i)),id,amon(im),iy, 1 int(xh),int(xm),int(xs),bet(i),comment c write(6,02) j,what(iphase(j)),id1,amon(im1),iy1, c 1 int(xh1),int(xm1),xs1 c end if c end if c end if enddo write(6,*) goto 100 end subroutine moon_phases(iyear,tjd,iphase,bet,nmax,nfound) implicit none save c..this routine computes the julian dates of lunar phases c..input: c..iyear = the year of interest c..nmax = dimension of the output arrays c..output: c..tjd(nmax) = ut julian dates of lunar phases. c..iphase(nmax) = phase of the moon for each of the ut julian dates c..bet(nmax) = ecliptic angle for each phase c..nfound = number of dates, phases, and angles found c..declare the pass integer iyear,nmax,nfound integer iphase(nmax) double precision tjd(nmax),bet(nmax) c..for solving the angle equation external msangle logical succes integer niter double precision low,high,msangle,tol,zbrent parameter (tol = 1.0d-6) c..for communicating which lunar phase integer imsangle double precision betcom common /solc11/ betcom,imsangle c..local variables integer nph,iy,im,id,iy1,im1,id1,i,imonth parameter (imonth =1) double precision tjdx,rhour,xh,xm,xs,rhour1 c..initialize nfound = 0 c..get lots of phases do nph=1,55 c..make a darn good guess call guess_moon(iyear,imonth,nph,tjdx) c call caldat(tjdx,iy,im,id,rhour) c..communicate which lunar phase we are doing imsangle = mod(nph,4) if (imsangle .eq. 0) imsangle=4 c..make sure the root is bracketed low = tjdx - 0.05d0 high = tjdx + 0.05d0 call zbrac(msangle,low,high,succes) if (.not. succes) then write(6,*) 'could not bracket in routine moon phases' write(6,*) 'low =',low,' high=',high end if c..solve for the ut julian date of the lunar phase tjdx = zbrent(msangle,low,high,tol,niter) c call caldat(tjdx,iy1,im1,id1,rhour1) c..debug for potential errant root finds c..i've checked it for the limits of the de405 data set c if (iy1 .ne. iy .or. im1 .ne. im1 .or. id1 .ne. id) then c write(6,*) c write(6,*) 'potential mishap' c write(6,222) iy,im,id,rhour,imsangle c write(6,222) iy1,im1,id1,rhour1,nfound c 222 format(1x,i4,' ',i2,' ',i2,' ',1pe12.4,' ',i2) c write(6,*) c read(5,*) c end if c..if its within the year, store the date and phase code call caldat(tjdx,iy,im,id,rhour) if (iy .eq. iyear) then nfound = nfound + 1 tjd(nfound) = tjdx iphase(nfound) = imsangle bet(nfound) = betcom end if enddo return end double precision function msangle(tjd) implicit none save c.for determining exact times of lunar phases c..the ut julian date is iterated on by a root finder. c..declare double precision tjd c..for communicating which lunar phase, and ecliptic angle integer imsangle double precision betcom common /solc11/ betcom,imsangle c..for getting the positional data character*40 aname integer isun,imoon,icoord parameter (isun=1,imoon=2,icoord = 1) double precision glon,glat,height, 1 xequ,yequ,zequ,vxequ,vyequ,vzequ, 2 ra,dec,rah,ram,ras,decd,decm,decs, 3 xecl,yecl,zecl,vxecl,vyecl,vzecl, 4 lam,bet,lamd,lamm,lams,betd,betm,bets, 5 xalt,yalt,zalt,vxalt,vyalt,vzalt, 6 alt,azi,altd,altm,alts,azid,azim,azis, 7 dist,drdt,dsun,dsundt, 8 xmanom,semi_maj,qper,xecen,lan,xinc,aop, 9 period,tjdp parameter (glon = 0.0d0, 2 glat = 0.0d0, 3 height = 0.0d0) c..it doesn't matter if one works in equatorial or ecliptic space, c..so i'll use the ecliptic angle lambda. c.. 0 degrees between sun and moon is a new moon , c.. 90 degrees between sun and moon is a first quarter moon, c..180 degrees between sun and moon is a full moon, c..270 degrees between sun and moon is a trird quarter moon. double precision lammoon,lamsun,diff,angle(4) data angle /0.0d0, 90.0d0, 180.0d0, 270.0d0/ c..get the position of the sun call position( 1 isun,tjd,icoord,glon,glat,height, 2 aname, 3 xequ,yequ,zequ,vxequ,vyequ,vzequ, 4 ra,dec,rah,ram,ras,decd,decm,decs, 5 xecl,yecl,zecl,vxecl,vyecl,vzecl, 6 lam,bet,lamd,lamm,lams,betd,betm,bets, 7 xalt,yalt,zalt,vxalt,vyalt,vzalt, 8 alt,azi,altd,altm,alts,azid,azim,azis, 9 dist,drdt,dsun,dsundt, & semi_maj,qper,xecen,lan,xinc,aop,xmanom,period,tjdp) lamsun = lam c..get the position of the moon call position( 1 imoon,tjd,icoord,glon,glat,height, 2 aname, 3 xequ,yequ,zequ,vxequ,vyequ,vzequ, 4 ra,dec,rah,ram,ras,decd,decm,decs, 5 xecl,yecl,zecl,vxecl,vyecl,vzecl, 6 lam,bet,lamd,lamm,lams,betd,betm,bets, 7 xalt,yalt,zalt,vxalt,vyalt,vzalt, 8 alt,azi,altd,altm,alts,azid,azim,azis, 9 dist,drdt,dsun,dsundt, & semi_maj,qper,xecen,lan,xinc,aop,xmanom,period,tjdp) lammoon = lam betcom = bet c..set the difference between the sun and moon depending on imsangle diff = lammoon - lamsun c..various pathologies to take care of c..because 360 degrees maps back to zero degrees if (imsangle .eq. 1) then if (diff .lt. -180.0) diff = diff + 360.0d0 if (diff .gt. 180.0) diff = diff - 360.0d0 else if (imsangle .eq. 2) then if (diff .lt. 0.0) diff = diff + 360.0d0 if (diff .gt. 180.0) diff = diff - 180.0d0 else if (imsangle .eq. 3) then if (diff .lt. 0.0) diff = diff + 360.0d0 else if (imsangle .eq. 4) then if (diff .lt. 0.0) diff = diff + 360.0d0 if (diff .lt. 180.0) diff = diff + 180.0d0 end if c..and set the equation we seek the root of msangle = diff - angle(imsangle) c write(6,222) lammoon,lamsun,lammoon-lamsun,diff,msangle c 222 format(1x,1p6e14.6) c read(5,*) return end subroutine guess_moon(iyear,imonth,nph,tjde) implicit none save c..this routine makes a darn good guess for the et julian date c..of the moons phases its maximum error is about 30 seconds c..between 1980 and 2020. adapted from meeus page 319. c..input: c..iyear = year c..imonth = month c..nph = type of lunar phase c. 1 = new moon c.. 2 = first quarter c.. 3 = full moon c.. 4 = third quarter c..output c..tjde = et julian date of the desired lunar phase c..declare the pass integer iyear,imonth,nph double precision tjde c..local variables integer i,iwnt double precision rk,t,t2,t3,t4,xe,xm,xmp,xf,xom,a(25), 1 xe2,f1,f2,w,ryear c..initialize ryear = dble(iyear) + dble(30.0d0*(imonth-2))/365.0d0 iwnt = nint( (ryear - 2000.0d0) * 12.3685d0) rk = dble(iwnt) + 0.25d0 * dble(nph-1) t = rk/1236.85d0 t2 = t*t t3 = t2*t t4 = t3*t c..initial guess at the julian date tjde = 2451550.09765d0 + 29.530588853 * rk 1 + 1.337d-4 * t2 2 - 1.50d-7 * t3 3 + 7.3d-10 * t4 c..largest set of corrections are angle dependent, they are xe = 1.0d0 - 2.516d-3 * t - 7.4d-6 * t2 xe2 = xe * xe c..estimate of the sun's mean anomaly at tjde xm = 2.5534 + 29.10535669d0 * rk 1 - 2.18d-5 * t2 2 - 1.1d-7 * t3 xm = mod(xm,360.0d0) if (xm .lt. 0.0) xm = xm + 360.0d0 c..estimate of the moon's mean anomaly at tjde xmp = 201.5643 + 385.81693528d0 * rk 1 + 1.07438d-2 * t2 2 + 1.239d-5 * t3 3 - 5.8d-8 * t4 xmp = mod(xmp,360.0d0) if (xmp .lt. 0.0) xmp = xmp + 360.0d0 c..estimate moon's argument of latitude xf = 160.7108d0 + 390.67050274d0 * rk 1 - 1.6341d-3 * t2 2 - 2.27d-6 * t3 3 + 1.1d-8 * t4 xf = mod(xf,360.0d0) if (xf .lt. 0.0) xf = xf + 360.0d0 c..estimate moon's longitude of the ascending node xom = 124.7746 - 1.56375580d0 * rk 1 + 2.0691d-3 * t2 2 + 2.15d-6 * t3 xom = mod(xom,360.0d0) if (xom .lt. 0.0) xom = xom + 360.0d0 c..new moon and full moon corrections if (mod(nph,4) .eq. 1 .or. mod(nph,4) .eq. 3) then if (mod(nph,4) .eq. 1) then a(1) = -0.40720d0 * sind(xmp) a(2) = 0.17241d0 * xe * sind(xm) a(3) = 0.01608d0 * sind(2.0d0*xmp) a(4) = 0.01039d0 * sind(2.0d0*xf) a(5) = 0.00739d0 * xe * sind(xmp - xm) a(6) = -0.00514d0 * xe * sind(xmp + xm) a(7) = 0.00208d0 * xe2 * sind(2.0d0*xm) else if (mod(nph,4) .eq. 3) then a(1) = -0.40614d0 * sind(xmp) a(2) = 0.17302d0 * xe * sind(xm) a(3) = 0.01614d0 * sind(2.0d0*xmp) a(4) = 0.01043d0 * sind(2.0d0*xf) a(5) = 0.00734d0 * xe * sind(xmp - xm) a(6) = -0.00515d0 * xe * sind(xmp + xm) a(7) = 0.00209d0 * xe2 * sind(2.0d0*xm) end if a(8) = -0.00111d0 * sind(xmp - 2.0d0*xf) a(9) = -0.00057d0 * sind(xmp + 2.0d0*xf) a(10) = 0.00056d0 * xe * sind(2.0d0*xmp + xm) a(11) = -0.00042d0 * sind(3.0d0*xmp) a(12) = 0.00042d0 * xe * sind(xm + 2.0d0*xf) a(13) = 0.00038d0 * xe * sind(xm - 2.0d0*xf) a(14) = -0.00024d0 * xe * sind(2.0d0*xmp - xm) a(15) = -0.00017d0 * sind(xom) a(16) = -0.00007d0 * sind(xmp + 2.0d0*xm) a(17) = 0.00004d0 * sind(2.0d0*xmp - 2.0d0*xf) a(18) = 0.00004d0 * sind(3.0d0*xm) a(19) = 0.00003d0 * sind(xmp + xm - 2.0d0*xf) a(20) = 0.00003d0 * sind(2.0d0*xmp + 2.0d0*xf) a(21) = -0.00003d0 * sind(xmp + xm + 2.0d0*xf) a(22) = 0.00003d0 * sind(xmp - xm + 2.0d0*xf) a(23) = -0.00002d0 * sind(xmp - xm - 2.0d0*xf) a(24) = -0.00002d0 * sind(3.0d0*xmp + xm) a(25) = 0.00002d0 * sind(4.0d0*xmp) w = 0.0d0 c..first and last quarter correction terms else if (mod(nph,4) .eq. 2 .or. mod(nph,4) .eq. 0) then a(1) = -0.62801d0 * sind(xmp) a(2) = 0.17172d0 * xe * sind(xm) a(3) = -0.01183d0 * xe * sind(xmp + xm) a(4) = 0.00862d0 * sind(2.0d0*xmp) a(5) = 0.00804d0 * sind(2.0d0*xf) a(6) = 0.00454d0 * xe * sind(xmp - xm) a(7) = 0.00204d0 * xe2 * sind(2.0d0*xm) a(8) = -0.00180d0 * sind(xmp - 2.0d0*xf) a(9) = -0.00070d0 * sind(xmp + 2.0d0*xf) a(10) = -0.00040d0 * sind(3.0d0*xmp) a(11) = -0.00034d0 * xe * sind(2.0d0*xmp - xm) a(12) = 0.00032d0 * xe * sind(xm + 2.0d0*xf) a(13) = 0.00032d0 * xe * sind(xm - 2.0d0*xf) a(14) = -0.00028d0 * xe2 * sind(xmp + 2.0d0*xm) a(15) = 0.00027d0 * xe * sind(2.0d0*xmp + xm) a(16) = -0.00017d0 * sind(xom) a(17) = -0.00005d0 * sind(xmp - xm - 2.0d0*xf) a(18) = 0.00004d0 * sind(2.0d0*xmp + 2.0d0*xf) a(19) = -0.00004d0 * sind(xmp + xm + 2.0d0*xf) a(20) = 0.00004d0 * sind(xmp - 2.0d0*xm) a(21) = 0.00003d0 * sind(xmp + xm - 2.0d0*xf) a(22) = 0.00003d0 * sind(3.0d0*xm) a(23) = 0.00002d0 * sind(2.0d0*xmp - 2.0d0*xf) a(24) = 0.00002d0 * sind(xmp - xm + 2.0d0*xf) a(25) = -0.00002d0 * sind(3.0d0*xmp + xm) w = 0.00306d0 - 0.00038d0 * xe * cos(xm) + 0.00026d0 * cos(xmp) 1 - 0.00002d0 * cos(xmp - xm) + 0.00002d0 * cos(xmp + xm) 2 + 0.00002d0 * cos(2.0d0*xf) if (mod(nph,4) .eq. 0) w = -w end if c..sum the 25 correction terms f1 = 0.0d0 do i=1,25 f1 = f1 + a(i) enddo c..a smaller set of corrections a(1) = 299.77d0 + 0.107408d0 * rk - 0.009173d0 * t2 a(2) = 251.88d0 + 0.016321d0 * rk a(3) = 251.83d0 + 26.651886d0 * rk a(4) = 349.42d0 + 36.412478d0 * rk a(5) = 84.66d0 + 18.206239d0 * rk a(6) = 141.74d0 + 53.303771d0 * rk a(7) = 207.14d0 + 2.453732d0 * rk a(8) = 154.84d0 + 7.306860d0 * rk a(9) = 34.52d0 + 27.261239d0 * rk a(10) = 207.19d0 + 0.121824d0 * rk a(11) = 291.34d0 + 1.844379d0 * rk a(12) = 161.72d0 + 24.198154d0 * rk a(13) = 239.56d0 + 25.513099d0 * rk a(14) = 331.55d0 + 3.592518d0 * rk do i=1,14 a(i) = mod(a(i),360.0d0) enddo do i=1,14 if (a(i) .lt. 0.0) a(i) = a(i) + 360.0d0 enddo a(1) = 0.000325d0 * sind(a(1)) a(2) = 0.000165d0 * sind(a(2)) a(3) = 0.000164d0 * sind(a(3)) a(4) = 0.000126d0 * sind(a(4)) a(5) = 0.000110d0 * sind(a(5)) a(6) = 0.000062d0 * sind(a(6)) a(7) = 0.000060d0 * sind(a(7)) a(8) = 0.000056d0 * sind(a(8)) a(9) = 0.000047d0 * sind(a(9)) a(10) = 0.000042d0 * sind(a(10)) a(11) = 0.000040d0 * sind(a(11)) a(12) = 0.000037d0 * sind(a(12)) a(13) = 0.000035d0 * sind(a(13)) a(14) = 0.000023d0 * sind(a(14)) f2 = 0.0d0 do i=1,14 f2 = f2 + a(i) enddo c..add in the correction terms tjde = tjde + f1 + f2 + w return end include 'position_routine.f' include 'glue.f'