program testper
      implicit none
      save

c..exercises the perihelion/aphelion (peigee/apogee for the moon) routine

c..declare 
      character*40     aname
      character*15     what(4)
      character*3      amon(12)
      integer          ibody,iyear,i,j,iy,im,id,nmax,nfound,k,nwant
      parameter        (nmax=10000)
      integer          type(nmax)
      double precision tjd(nmax),dist(nmax),rhour,xh,xm,xs,fac


c..for picking which ephemeris to use
      integer          which_eph
      common /dop/     which_eph


c..initializations
      data        what /'perihelion','aphelion',
     1                  'perigee   ','apogee  '/

      data       amon /'jan' , 'feb' , 'mar' , 'apr' , 'may' , 'jun' , 
     1                 'jul' , 'aug' , 'sep' , 'oct' , 'nov' , 'dec' /


c..popular format statements
 01   format(1x,t7,a,t24,a,t34,a,t37,a,t40,a,t52,a)
 02   format(1x,i4.4,' ',a,' ',i2.2,a,i5.4,' ',i2.2,i3.2,0pf7.3,
     1          ' UT ',1pe14.6)



c..use the long ephemeris de406 so a few perihelions of the outer planets
c..as well as  historial moon extrema can be found.

      which_eph = 406

      call init_names(1)



c..get the input 
 100  write(6,*) 'give the body year and number of peri you want=>'
      read(5,*) ibody,iyear,nwant
      nwant = 2*nwant 

      call perihelion(ibody,iyear,tjd,dist,type,aname,nmax,
     1                nwant,nfound)



c..write the results
c..uncommenting the code below will print out any extrema
c..perigees or apogees for the moon. see table 48.c of meuss page 322.

c..for the moon convert au to km
      fac = 1
      if (ibody .eq. 2) fac = 1.49597870691d8

      write(6,*)
      write(6,01) aname
      write(6,01) 'event','date','hr','mn','sec','distance'

      k = 0

      do i=1,nfound
       call caldat2(tjd(i),iy,im,id,rhour)
       call dh2h(rhour,xh,xm,xs)

c       if (type(i) .eq. 1  .and. dist(i)*fac .lt. 356425.0) then
c       if (type(i) .eq. 2  .and. dist(i)*fac .gt. 406710.0) then

       if (ibody .eq. 2) then
        write(6,02) i,what(type(i)+2),id,amon(im),iy,
     1              int(xh),int(xm),xs,dist(i)*fac
       else

c..aphelions/perihelion only
c         if (type(i) .eq. 2) then        
         if (type(i) .eq. 1) then        

          k = k + 1
          write(6,02) k,what(type(i)),id,amon(im),iy,
     1                int(xh),int(xm),xs,dist(i)*fac

         end if

       end if


c       end if

      enddo 
      write(6,*)

      goto 100
      end





      subroutine perihelion(ibody,iyear,tjd,dist,type,aname,nmax,
     1                      nwant,nfound)
      implicit none
      save

c..this routine does a calculation of julian days for perihelion and aphelion
c..or in the case of the moon, perigee and apogee

c..input:
c..ibody  = body of interest
c..iyear  = the year of interest
c..imonth = month of interest
c..nwant  = number of perihelion/aphelions to find

c..output:
c..tjd    = ut julian dates of perihelion/aphelion (perigee/apogee for moon)
c..dist   = distance from the sun in au for perihelion/aphelions
c..         distance from earth in au for perigee/apogee of moon
c..type   = type of event 1=perihelion 2=aphelion 
c..aname  = name of the object
c..nmax   = size of the tjd and dist arrays
c..nfound = number of perihelion and aphelions found


c..declare the pass
      character*40     aname
      integer          ibody,iyear,nmax,nwant,nfound
      integer          type(nmax)
      double precision tjd(nmax),dist(nmax)


c..for communicating which equinox or solstice
      character*40     xname
      integer          ibod,iwhat
      common /perc11/  ibod,iwhat,xname


c..for picking which ephemeris to use
      integer          which_eph
      common /dop/     which_eph

c..for finding the minima
      external         perifunc
      integer          niter
      double precision ax,bx,cx,fa,fb,fc,xmin,perifunc,mnbrent,tol
      parameter        (tol = 1.0d-14)


c..local variables
      integer          i,imonth,iy,im,id
      parameter        (imonth = 1)
      double precision tjdx,rhour


c..initialize
      nfound = 0

c..get lots of perihelions and aphelions
c      if (ibody .eq. 2) then
c       iend = 40
c      else if (ibody .le. 6) then
c       iend= 20
c       iend= 162
c      else
c       iend = 9
c      end if



      do i = 1,nwant


c..make a fair guess
       call guess_peri(i,ibody,iyear,imonth,tjdx)


c..bail if the implied year is greater than the limits of the ephemris       
       call caldat2(tjdx,iy,im,id,rhour)
        if (which_eph .eq. 405) then
         if (iy .lt. 1600 .or. iy .ge. 2200) return
        else if (which_eph .eq. 406) then
         if (iy .lt. -3000 .or. iy .ge. 3000) return
        end if 


c..communicate what body and what perihelion or aphelion we are doing
c..map the earth body code into the sun body code
       ibod   = ibody
       if (ibody .eq. 5) ibod = 1
       iwhat  = mod(i,2)
       if (iwhat .eq. 0) iwhat = 2
      

c..bracket the minima
       ax = tjdx - 1.0d0
       bx = tjdx + 1.0d0
       call mnbrak(ax,bx,cx,fa,fb,fc,perifunc)


c..get the minima
       xmin = mnbrent(ax,bx,cx,perifunc,tol,tjdx,niter) 


c..for the moon, mercury, venus, and earth 
c..only store it if its within the specified year
       if (ibod.eq.2 .or. ibod.eq.3 .or. ibod.eq.4 .or. ibod.eq.1) then
        call caldat2(tjdx,iy,im,id,rhour)

c        if (iy .eq. iyear) then
         nfound       = nfound + 1
         if (nfound .gt. nmax) stop 'nfound > nmax in equinox'
         tjd(nfound)  = tjdx
         type(nfound) = iwhat
c         write(6,*) nfound

         if (iwhat .eq. 1) then
          dist(nfound) = xmin
         else if (iwhat .eq. 2) then
          dist(nfound) = 1.0d0/xmin 
         endif
c        end if



c..for any other body, just store it
       else

        nfound        = nfound + 1
        if (nfound .gt. nmax) stop 'nfound > nmax in equinox'
        tjd(nfound)   = tjdx
        type(nfound)  = iwhat
c        write(6,*) nfound
 
        if (iwhat .eq. 1) then
         dist(nfound) = xmin
        else if (iwhat .eq. 2) then
         dist(nfound) = 1.0d0/xmin 
        endif
       end if

c..save the name of the object
       aname   = xname
       if (ibody .eq. 5) aname = 'earth'


c..and back for another perihelion or aphelion
      enddo
      return
      end





      double precision function perifunc(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
      character*40     xname
      integer          ibod,iwhat
      common /perc11/  ibod,iwhat,xname



c..for getting the positional data
      character*40     aname
      integer          icoord
      parameter        (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..get the position of the body
      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)

      xname = aname


c..for perihelion
      if (iwhat .eq. 1) then

c..for the earth, which was mapped into the sun body number, or the moon
       if (ibod .eq. 1 .or. ibod .eq. 2) then
        perifunc = dist

c..for any other planet
       else 
        perifunc = dsun
       end if

c..for aphelion
      else if (iwhat .eq. 2) then
       if (ibod .eq. 1 .or. ibod .eq. 2) then
        perifunc = 1.0d0/dist
       else
        perifunc = 1.0d0/dsun
       endif
      else
       stop 'blown iwhat in function perifunc'
      end if

      return
      end





      subroutine guess_peri(icode,ibody,iyear,imonth,tjde)
      implicit none
      save

c..this routine makes a fair guess for the et julian date of
c..perihelion and aphelion of the planets, or perigee and apogee 
c..in the case of the moon. the dates can be in error by hours
c..for the planets and within minutes for the moon.
c..adapted rom meuss page 253 and page 325.



c..input:
c..icode  = type of event
c..         1 = perihelion
c..         2 = aphelion
c..ibody  = planet of interest
c..iyear  = year of interest


c..output
c..tjde = et julian date of the perihelion or aphelion


c..declare the pass
      integer          icode,ibody,iyear,imonth
      double precision tjde


c..local variables
      integer          k
      double precision ryear,rk,t,t2,t3,t4,d,m,f,f1


c..initialize
       ryear = dble(iyear) + dble(30.0d0*(imonth-3))/365.0d0


c..moon
      if (ibody .eq. 2) then
       k    = 13.2555d0 * (ryear - 1999.97d0)
       rk   = dble(k) + 0.5d0 * dble(icode - 1)
       t    = rk/1325.55d0
       t2   = t**2
       t3   = t * t2
       t4   = t2 * t2

c..this is not a bad guess as it is, but we'll make a few refinements
       tjde = 2451534.6698d0 + 27.55454988d0 * rk
     1                       - 0.0006886d0 * t2
     2                       - 0.000001098d0 * t3
     3                       + 0.0000000052d0 * t4

c..mean elongation, mean anomaly, and argument of latitude
       d = 171.9179d0 + 335.9106046 * rk
     1                - 0.0100250d0 * t2
     2                - 0.00001156d0 * t3
     3                + 0.000000055d0 * t4
       m = 347.3477d0 + 27.1577721d0 * rk
     1                - 0.0008323d0 * t2
     2                - 0.0000010d0 * t3
       f = 316.6109d0 + 364.5287911d0 * rk
     1                - 0.0125131d0 * t2
     2                - 0.0000148d0 * t3

c..for perigees
      if (mod(icode,2) .eq. 1) then
       f1  = -1.6769d0 * sind(2.0d0*d)
     1      + 0.4589d0 * sind(4.0d0*d)
     2      - 0.1896d0 * sind(6.0d0*d)
     3      + 0.0883d0 * sind(8.0d0*d)
     4      + (-0.0773d0 + 0.00019d0*t) * sind(2.0d0*d - m)
     5      + ( 0.0502d0 - 0.00013d0*t) * sind(m)
     6      + 0.0460d0 * sind(10.0d0*d)
     7      + ( 0.0422d0 - 0.00011d0*t) * sind(4.0d0*d - m)
     8      - 0.0256d0 * sind(6.0d0*d - m)
     9      + 0.0253d0 * sind(12.0d0*d)
     &      + 0.0237d0 * sind(d)
     1      + 0.0162d0 * sind(8.0d0*d - m)
     2      - 0.0145d0 * sind(14.0d0*d)
     3      + 0.0129d0 * sind(2.0d0*f)
     4      - 0.0112d0 * sind(3.0d0*d)
     5      - 0.0104d0 * sind(10.0d0*d - m)
     6      + 0.0086d0 * sind(16.0d0*d)
       f1   =  f1
     1      + 0.0069d0 * sind(12.0d0*d - m)
     2      + 0.0066d0 * sind(5.0d0*d)
     3      - 0.0053d0 * sind(2.0d0*d + 2.0d0*f)
     4      - 0.0052d0 * sind(18.0d0*d)
     5      - 0.0046d0 * sind(14.0d0*d - m)
     6      - 0.0041d0 * sind(7.0d0*d)
     7      + 0.0040d0 * sind(2.0d0*d + m)
     8      + 0.0032d0 * sind(20.0d0*d)
     9      - 0.0032d0 * sind(d + m)
     &      + 0.0031d0 * sind(16.0d0*d - m)
     1      - 0.0029d0 * sind(4.0d0*d + m)
     2      + 0.0027d0 * sind(9.0d0*d)
     3      + 0.0027d0 * sind(4.0d0*d + 2.0d0*f)
     4      - 0.0027d0 * sind(2.0d0*d - 2.0d0*m)
     5      + 0.0024d0 * sind(4.0d0*d - 2.0d0*m)
     6      - 0.0021d0 * sind(6.0d0*d - 2.0d0*m)
       f1   =  f1
     1      - 0.0021d0 * sind(22.0d0*d)
     2      - 0.0021d0 * sind(18.0d0*d - m)
     3      + 0.0019d0 * sind(6.0d0*d + m)
     4      - 0.0018d0 * sind(11.0d0*d)
     5      - 0.0014d0 * sind(8.0d0*d + m)
     6      - 0.0014d0 * sind(4.0d0*d - 2.0d0*f)
     7      - 0.0014d0 * sind(6.0d0*d + 2.0d0*f)
     8      + 0.0014d0 * sind(3.0d0*d + m)
     9      - 0.0014d0 * sind(5.0d0*d + m)
     &      + 0.0013d0 * sind(13.0d0*d)
     1      + 0.0013d0 * sind(20.0d0*d - m)
     2      + 0.0011d0 * sind(3.0d0*d + 2.0d0*m)
     3      - 0.0011d0 * sind(4.0d0*d + 2.0d0*f - 2.0d0*m)
     4      - 0.0010d0 * sind(d + 2.0d0*m)
     5      - 0.0009d0 * sind(22.0d0*d - m)
     6      - 0.0008d0 * sind(4.0d0*f)
       f1   =  f1
     1      + 0.0008d0 * sind(6.0d0*d - 2.0d0*f)
     2      + 0.0008d0 * sind(2.0d0*d - 2.0d0*f + m)
     3      + 0.0007d0 * sind(2.0d0*m)
     4      + 0.0007d0 * sind(2.0d0*f - m)
     5      + 0.0007d0 * sind(2.0d0*d + 4.0d0*f)
     6      - 0.0006d0 * sind(2.0d0*f - 2.0d0*m)
     7      - 0.0006d0 * sind(2.0d0*d - 2.0d0*f + 2.0d0*m)
     8      + 0.0006d0 * sind(24.0d0*d)
     9      + 0.0005d0 * sind(4.0d0*d - 4.0d0*f)
     &      + 0.0005d0 * sind(2.0d0*d + 2.0d0*m)
     1      - 0.0004d0 * sind(d - m)

c..for apogees
      else if (mod(icode,2) .eq. 0) then
       f1  =  0.4392d0 * sind(2.0d0*d)
     1      + 0.0684d0 * sind(4.0d0*d)
     2      + (0.0456d0 - 0.00011d0*t) * sind(m)
     3      + (0.0426d0 - 0.00011d0*t) * sind(2.0d0*d - m)
     4      + 0.0212d0 * sind(2.0d0*f)
     5      - 0.0189d0 * sind(d)
     6      + 0.0144d0 * sind(6.0d0*d)
     7      + 0.0113d0 * sind(4.0d0*d - m)
     8      + 0.0047d0 * sind(2.0d0*d + 2.0d0*f)
     9      + 0.0036d0 * sind(d + m)
     &      + 0.0035d0 * sind(8.0d0*d)
     1      + 0.0034d0 * sind(6.0d0*d - m)
     2      - 0.0034d0 * sind(2.0d0*d - 2.0d0*f)
     3      + 0.0022d0 * sind(2.0d0*d - 2.0d0*m)
     4      - 0.0017d0 * sind(3.0d0*d)
     5      + 0.0013d0 * sind(4.0d0*d + 2.0d0*f)
     6      + 0.0011d0 * sind(8.0d0*d - m)
       f1  = f1
     1      + 0.0010d0 * sind(4.0d0*d - 2.0d0*m)
     2      + 0.0009d0 * sind(10.0d0*d)
     3      + 0.0007d0 * sind(3.0d0*d + m) 
     4      + 0.0006d0 * sind(2.0d0*m)
     5      + 0.0005d0 * sind(2.0d0*d + m)
     6      + 0.0005d0 * sind(2.0d0*d + 2.0d0*m)
     7      + 0.0004d0 * sind(6.0d0*d + 2.0d0*f)
     8      + 0.0004d0 * sind(6.0d0*d - 2.0d0*m)
     9      + 0.0004d0 * sind(10.0d0*d - m)
     &      - 0.0004d0 * sind(5.0d0*d)
     1      - 0.0004d0 * sind(4.0d0*d - 2.0d0*f)
     2      + 0.0003d0 * sind(2.0d0*f + m)
     3      + 0.0003d0 * sind(12.0d0*d)
     4      + 0.0003d0 * sind(2.0d0*d + 2.0d0*f - m)
     5      - 0.0003d0 * sind(d - m)
       end if



c..add in the correction term
       tjde = tjde + f1


c..mercury
      else if (ibody .eq. 3) then
       k     = 4.15201d0 * (ryear - 2000.12d0) 
       rk    = dble(k) + 0.5d0 * dble(icode - 1)
       tjde  = 2451590.257d0  + 87.96934963d0 * rk

c..venus
      else if (ibody .eq. 4) then
       k     = 1.62549d0 * (ryear - 2000.53d0) 
       rk    = dble(k) + 0.5d0 * dble(icode - 1)
       tjde  = 2451738.233d0 
     1         + 224.700818d0 * rk
     2         - 0.0000000327d0 * rk**2    

c..earth
      else if (ibody .eq. 5) then
       k     = 0.99997d0 * (ryear - 2000.01d0) 
       rk    = dble(k) + 0.5d0 * dble(icode - 1)
       tjde  = 2451547.507d0
     1         + 365.2596358d0 * rk
     2         - 0.0000000158d0 * rk**2

c..mars
      else if (ibody .eq. 6) then
       k     = 0.53166d0 * (ryear - 2000.78d0) 
       rk    = dble(k) + 0.5d0 * dble(icode - 1)
       tjde  = 2452195.026d0
     1         + 686.995784d0 * rk
     2         - 0.0000001187d0 * rk**2

c..jupiter
      else if (ibody .eq. 7) then
       k     = 0.08430d0 * (ryear - 2011.20d0) 
       rk    = dble(k) + 0.5d0 * dble(icode - 1)
       tjde  = 2455636.938d0
     1         + 4332.897080d0 * rk
     2         + 0.0001368d0 * rk**2

c..saturn
      else if (ibody .eq. 8) then
       k     = 0.03393d0 * (ryear - 2003.52d0) 
       rk    = dble(k) + 0.5d0 * dble(icode - 1)
       tjde  = 2452830.11d0
     1         + 10764.21731d0 * rk
     2         + 0.000826d0 * rk**2

c..uranus
      else if (ibody .eq. 9) then
       k     = 0.01190d0 * (ryear - 2051.1d0) 
       rk    = dble(k) + 0.5d0 * dble(icode - 1)
       tjde  = 2470213.5d0
     1         + 30694.8767d0 * rk
     2         - 0.00541d0 * rk**2


c..neptune
      else if (ibody .eq. 10) then
       k     = 0.00607d0 * (ryear - 2047.5d0) 
       rk    = dble(k) + 0.5d0 * dble(icode - 1)
       tjde  = 2468895.7d0
     1         + 60190.32d0 * rk
     2         + 0.3175d0 * rk**2

c..pluto
      else if (ibody .eq. 11) then
       k     = 0.00607d0 * (ryear - 2047.5d0) 
       rk    = dble(k) + 0.5d0 * dble(icode - 1)
       tjde  = 2468895.7d0
     1         + 60190.32d0 * rk
     2         + 0.3175d0 * rk**2


c..a blown body number
      else
       stop 'unknown icode in routine guess_peri'
      end if

      return
      end



      include 'position_routine.f'
      include 'glue.f'