program testequ implicit none save c..exercises the equinox and solstice routine c..declare character*15 what(4) character*3 amon(12) integer iyear,i,iy,im,id double precision tjd(4),rhour,xh,xm,xs c..for picking which ephemeris to use integer which_eph common /dop/ which_eph data what /'spring equinox ','summer solstice', 1 'autum equinox ','winter solstice'/ data amon /'jan' , 'feb' , 'mar' , 'apr' , 'may' , 'jun' , 1 'jul' , 'aug' , 'sep' , 'oct' , 'nov' , 'dec' / c..popular format statements 01 format(1x,t2,a,t18,a,t28,a,t31,a,t34,a) 02 format(1x,a,' ',i2.2,a,i4.4,' ',i2.2,i3.2,0pf7.3,' UT') c..initialize the name database which_eph = 406 call init_names(0) 100 write(6,*) 'give the year =>' read(5,*) iyear call equinox(iyear,tjd) c..write the results write(6,*) write(6,01) 'event','date','hr','mn','sec' do i=1,4 call caldat(tjd(i),iy,im,id,rhour) call dh2h(rhour,xh,xm,xs) write(6,02) what(i),id,amon(im),iy,int(xh),int(xm),xs enddo write(6,*) goto 100 end subroutine equinox(iyear,tjd) implicit none save c..this routine computes the julian dates of equinoxes and solstices c..input: c..iyear = the year of interest c..output: c..tjd = ut julian dates of spring equinox, summer solstice, c.. fall equinox, and winter solstice, in that order. c..declare the pass integer iyear double precision tjd(4) c..for solving angle equation external sunangle logical succes integer niter double precision low,high,sunangle,tol,zbrent parameter (tol = 1.0d-6) c..for communicating which equinox or solstice integer isunangle common /solc11/ isunangle c..local variables integer i double precision tjdx c..get each equinox or solstice of the year do i=1,4 c..make a darn good guess call guess_equinox(i,iyear,tjdx) c..communicate which equinox or solstice we are doing isunangle = i c..make sure the root is bracketed low = tjdx - 0.5d0 high = tjdx + 0.5d0 call zbrac(sunangle,low,high,succes) if (.not. succes) then write(6,*) 'could not bracket in routine equinox' write(6,*) 'low =',low,' high=',high end if c..solve for the ut julian date of the equinox or solstice tjd(i) = zbrent(sunangle,low,high,tol,niter) enddo return end double precision function sunangle(tjd) implicit none save c.for determining exact times of equinox and solstice c..the ut julian date is iterated on by a root finder. c..declare double precision tjd c..for communicating which equinox or solstice integer isunangle common /solc11/ isunangle c..for getting the positional data character*40 aname integer ibody,icoord parameter (ibody=1,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 is vernal equinox, 90 degrees is summer solstice c..180 degrees is fall equinox, 270 degrees is winter solstice double precision angle(4) data angle /0.0d0, 90.0d0, 180.0d0, 270.0d0/ c..get the position of the sun call position( 1 ibody,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) c..the spring equinox at zero degrees is tough on circular functions c..in this case move the branch cut to pi if (isunangle .eq. 1 .and. lam .gt. 180.0) lam = lam - 360.0d0 c..set the root depending on the value of isunangle sunangle = angle(isunangle) - lam return end subroutine guess_equinox(icode,iyear,tjde) implicit none save c..this routine makes a darn good guess for the et julian date c..of solstices and equinoxes. its maximum error is about 1 minute c..between 1951 and 2050. adapted from meeus page 165. c..input: c..icode = type of equinox or solstice desired c.. 1 = vernal equinox, c.. 2 = summer solstice, c.. 3 = fall equinox c. 4 = winter solstice c..iyear = year of interest c..output c..tjde = et julian date of the solstice or equinox c..declare the pass integer icode,iyear double precision tjde c..math constants double precision pi,a2rad,rad2a parameter (pi = 3.1415926535897932384d0, 1 a2rad = pi/180.0d0, 2 rad2a = 180.0d0/pi) c..local variables integer i,nmax parameter (nmax=24) double precision y,tjde0,t,w,dl,s,a(nmax),b(nmax),cx(nmax) c..fitting coefficients data a/ 1 485.0d0, 203.0d0, 199.0d0, 182.0d0, 156.0d0, 136.0d0, 2 77.0d0, 74.0d0, 70.0d0, 58.0d0, 52.0d0, 50.0d0, 3 45.0d0, 44.0d0, 29.0d0, 18.0d0, 17.0d0, 16.0d0, 4 14.0d0, 12.0d0, 12.0d0, 12.0d0, 9.0d0, 8.0d0/ data b/ 1 324.96d0, 337.23d0, 342.08d0, 27.85d0, 73.14d0, 171.52d0, 2 222.54d0, 296.72d0, 243.58d0, 119.81d0, 297.17d0, 21.02d0, 3 247.54d0, 325.15d0, 60.93d0, 155.12d0, 288.79d0, 198.04d0, 4 199.76d0, 95.39d0, 287.11d0, 320.81d0, 227.73d0, 15.45d0/ data cx/ 1 1934.136d0, 32964.467d0, 20.186d0, 445267.112d0, 2 45036.886d0, 22518.443d0, 65928.934d0, 3034.906d0, 3 9037.513d0, 33718.147d0, 150.678d0, 2281.226d0, 4 29929.562d0, 31555.956d0, 4443.417d0, 67555.328d0, 5 4562.452d0, 62894.029d0, 31436.921d0, 14577.848d0, 6 31931.756d0, 34777.259d0, 1222.114d0, 16859.074d0/ c..spring equinoxes are in march if (icode .eq. 1) then if (iyear .lt. 1000) then y = dble(iyear) * 1.0d-3 tjde0 = 1721139.29189d0 + y*(365242.13740d0 + y*(0.06134d0 + 1 y*(0.00111d0 - y*0.00071d0))) else y = dble(iyear-2000) * 1.0d-3 tjde0 = 2451623.80984d0 + y*(365242.37404d0 + y*(0.05169d0 + 1 y*(-0.00411d0 - y*0.00057d0))) end if c..summer solsticees are in june else if (icode .eq. 2) then if (iyear .lt. 1000) then y = dble(iyear) * 1.0d-3 tjde0 = 1721233.25401d0 + y*(365241.72562d0 + y*(-0.05323d0 + 1 y*(0.00907d0 + y*0.00025d0))) else y = dble(iyear-2000) * 1.0d-3 tjde0 = 2451716.56767d0 + y*(365241.62603d0 + y*(0.00325d0 + 1 y*(0.00888d0 - y*0.00030d0))) end if c..fall equinoxes are in september else if (icode .eq. 3) then if (iyear .lt. 1000) then y = dble(iyear) * 1.0d-3 tjde0 = 1721325.70455d0 + y*(365242.49558d0 + y*(-0.11677d0 + 1 y*(-0.00297d0 + y*0.00074d0))) else y = dble(iyear-2000) * 1.0d-3 tjde0 = 2451810.21715d0 + y*(365242.01767d0 + y*(-0.11575d0 + 1 y*(0.00337d0 + y*0.00078d0))) end if c..winter solstices are in december else if (icode .eq. 4) then if (iyear .lt. 1000) then y = dble(iyear) * 1.0d-3 tjde0 = 1721414.39987d0 + y*(365242.88257d0 + y*(-0.00769d0 + 1 y*(-0.00933d0 - y*(0.00006d0)))) else y = dble(iyear-2000) * 1.0d-3 tjde0 = 2451900.05952d0 + y*(365242.74049d0 + y*(-0.06223d0 + 1 y*(-0.00823d0 + y*0.00032d0))) end if c..a blown input else stop 'unknown icode in routine guess_equinox' end if c..sum terms t = (tjde0 - 2451545.0d0)/36525.0d0 w = 35999.373d0 * t - 2.47d0 dl = 1.0d0 + 3.34d-2 * cos(a2rad*w) + 7.0d-3 * cos(a2rad*2.0d0*w) s = 0.0d0 do i=1,24 s = s + a(i) * cos(a2rad*(b(i) + cx(i)*t)) enddo s = int(s) tjde = tjde0 + s * 1.0d-5/dl return end include 'position_routine.f'