program tretro implicit none save c..this program exercises the main position subroutine c..declare c..for calling the retrograde routine character*40 aname integer ibody,nmax,nretro parameter (nmax=100) double precision tjd_beg,tjd_end,tstep, 1 tjd_retro_beg(nmax),tjd_retro_end(nmax) c..for picking which ephemeris to use integer which_eph common /dop/ which_eph c..local variables character*3 amon(12) integer i,iyear_beg,imonth_beg,iday_beg, 1 iyear_end,imonth_end,iday_end, 2 iy,im,id,iy2,im2,id2 double precision rhr,xh,xm,xs,rhr2,xh2,xm2,xs2 c..data statements data amon /'jan' , 'feb' , 'mar' , 'apr' , 'may' , 'jun' , 1 'jul' , 'aug' , 'sep' , 'oct' , 'nov' , 'dec' / c..popular format statements 03 format(1x,i2.2,a,i4,' ',0pf9.4,0pf9.4,' ',a) 04 format(1x,a8,a,i2.2,a,i4,' ',i2.2,':',i2.2,a, 1 i2.2,a,i4,' ',i2.2,':',i2.2) c..initialize the name database which_eph = 405 call init_names(1) c..get the input write(6,*) 'give body number =>' read(5,*) ibody write(6,*) write(6,*) 'give the starting year month and day =>' read(5,*) iyear_beg,imonth_beg,iday_beg write(6,*) write(6,*) 'give the ending year month and day =>' read(5,*) iyear_end,imonth_end,iday_end write(6,*) write(6,*) 'give the tstep to search over' write(6,*) '1.0=1day 2.0d=2 days 0.5 = half days and so on =>' read(5,*) tstep c..get the julian dates of the start and end times call juldat(iyear_beg,imonth_beg,iday_beg,0.0d0,tjd_beg) call juldat(iyear_end,imonth_end,iday_end,0.0d0,tjd_end) call retrograde(ibody,tjd_beg,tjd_end,tstep, 1 aname, 2 tjd_retro_beg,tjd_retro_end,nmax,nretro) c..report the dates do i=1,nretro call caldat(tjd_retro_beg(i),iy,im,id,rhr) call dh2h(rhr,xh,xm,xs) call caldat(tjd_retro_end(i),iy2,im2,id2,rhr2) call dh2h(rhr2,xh2,xm2,xs2) write(6,04) aname,'retro from: ', 1 id,amon(im),iy, 2 int(xh),nint(xm+xs/60.0d0),' to ', 3 id2,amon(im2),iy2, 4 int(xh2),nint(xm2+xs2/60.0d0) enddo stop 'normal termination' end subroutine retrograde(ibody,tjd_beg,tjd_end,tstep, 1 aname, 2 tjd_retro_beg,tjd_retro_end,nmax,nretro) implicit none save c..this routine finds the time periods when a body undergoes retrograde motion c..input: c..ibody = body number c..tjd_beg = ut julian date of the start of the search interval c..tjd_end = ut julian date of the end of the search interval c..tstep = interval in days to conduct the search over c..output: c..aname = name of the body c..tjd_retro_beg(1:nmax) = array od ut julian dates when retrograde starts c..tjd_retro_end(1:nmax) = array od ut julian dates when retrograde ends c..nmax = size of the output arrays c..nretro = number of retrograde periods found c..declare the pass character*40 aname integer ibody,nmax,nretro double precision tjd_beg,tjd_end,tstep, 1 tjd_retro_beg(nmax),tjd_retro_end(nmax) c..for the positions integer iyear,imonth,iday,icoord double precision rhour,glon,glat,height,tjd, 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 semi_maj,qper,xecen,lan,xinc,aop,xmanom, 9 period,tjdp c..assume an geocentric origin parameter (icoord = 1, 1 glon = 0.0d0, 2 glat = 0.0d0, 3 height = 0.0d0) c..for communicating which part of the retrograde integer ibod,iwhat common /retc11/ ibod,iwhat c..local variables external retrofunc logical local_retro,in_retro integer i,j,nstep,niter double precision lamold,ax,bx,cx,fa,fb,fc,retrofunc, 1 mnbrent,xmin,tjdsav,tol parameter (tol = 1.0d-12) c..number of steps to take and other initializations nstep = (tjd_end - tjd_beg)/tstep + 1.0d0 in_retro = .false. nretro = 0 lamold = 0.0d0 tjdsav = 0.0d0 do i=1,nmax tjd_retro_beg(i) = 0.0d0 tjd_retro_end(i) = 0.0d0 enddo c..for every time step in the search window do i=1,nstep tjd = tjd_beg + float(i-1)*tstep c..get the state vectors 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..check if locally this point is retrograde c..and avoid false false retrogrades from periodicity local_retro = .false. if ( (lam-lamold) .le. 0.0 .and. 1 (lamold-lam) .le. 355.0) local_retro = .true. c..mark the begining of a retrograde period if (local_retro) then if (.not.in_retro) tjdsav = tjd in_retro = .true. end if c..note the end of a retrograde period, store it, reset the in_retro flag if (in_retro .and. .not.local_retro) then nretro = nretro + 1 if (nretro .gt. nmax) stop 'nretro > nmax' tjd_retro_beg(nretro) = tjdsav tjd_retro_end(nretro) = tjd in_retro = .false. end if c..save the current ecliptic nagle and go bak for another date lamold = lam enddo c..set the common block body number ibod = ibody c..find each extrema do i=1,nretro do j=1,2 iwhat = mod(j,2) if (iwhat .eq. 0) iwhat = 2 c..bracket the extrema if (iwhat .eq. 1) then ax = tjd_retro_beg(i) - 2.0d0*tstep bx = tjd_retro_beg(i) + 2.0d0*tstep else ax = tjd_retro_end(i) - 2.0d0*tstep bx = tjd_retro_end(i) + 2.0d0*tstep end if call mnbrak(ax,bx,cx,fa,fb,fc,retrofunc) c..get the extrema xmin = mnbrent(ax,bx,cx,retrofunc,tol,tjdsav,niter) if (iwhat .eq. 1) then tjd_retro_beg(i) = tjdsav else tjd_retro_end(i) = tjdsav end if c..back for another beg or end enddo enddo return end double precision function retrofunc(tjd) implicit none save c.for determining exact times of retrograde c..the ut julian date is iterated on by a minima finder. c..declare double precision tjd c..for communicating which part of the retrograde integer ibod,iwhat common /retc11/ ibod,iwhat c..for getting the positional data character*40 aname integer icoord 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 c..geocentric parameter (icoord = 1, 1 glon = 0.0d0, 2 glat = 0.0d0, 3 height = 0.0d0) c..get the objects coordinates call position( 1 ibod,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..for the beginning of retrograde if (iwhat .eq. 1) then retrofunc = 1.0d0/lam c..for the end of a retrograde else if (iwhat .eq. 2) then retrofunc = lam c..otherwise somethings not set right else stop 'blown iwhat in function retrofunc' end if return end include 'position_routine.f' include 'glue.f'