      program HW02_03

      implicit none

*     Program to solve the N-body gravitational problem using Runge-Kutta
*     numerical integration.
*     Arguments for the program are passed through runsting:
*     % HW02_03 <Body file> <duration (days)> <error tolerance>
*

      include 'NBody.h'

* MAIN Program variables
      
      real*8 t     ! Current time (days)
     .,      step  ! Current step size (days)
     .,      new_step  ! Updated step size that will yield desired accuracy

      real*8 poss(2,max_body), vels(2,max_body)  ! Position and velocity at start
                           ! of integration step
      real*8 posc(2,max_body), velc(2,max_body)  ! Position and velocity at end
                           ! of integration step when done as one step

      real*8 posa(2,max_body), vela(2,max_body)  ! Position and velocity at end
                           ! of first half/step integration
      real*8 posb(2,max_body), velb(2,max_body)  ! Position and velocity at end
                           ! of second half/step integration (posb should = posc 
                           ! within tolerance.
      real*8 minstep       ! Minum step size needed.

      integer*4 sc    ! Counter for step so that we don't get stuck



      logical step_not_OK  ! Logical set true while step size is still being
                           ! evaluated
     .,        no_step_increase  ! Set true to stop step from inceasing.


****  The runstring arguments passed by user
      call read_runstring( infile, dur, tol)

****  Get the information on the initial conditions
      call read_nbody

****  OK Start the integration
      step = 1.0d0
      minstep = step 
*     Initializa time
      t = 0.d0
      closest = 1.d20
*     Initialize starting position and velocity
      poss(:,1:num_body) = posi(:,1:num_body)
      vels(:,1:num_body) = veli(:,1:num_body)

*     Report initial conditions
      call report_ICS(t, poss, vels, 'Initial Conditions')
      call animate(t, step,  poss)

      do while ( t.lt. dur )

*        Nake sure that the next step will nor overshoot the end
*        of the time (this can happen because are are varying the
*        step size during the run.
         no_step_increase = .false.
         if ( t + step .gt.dur ) then
             step = dur - t
             no_step_increase = .true.
         endif

*        Iterate on step size to get the value needed to maintain the
*        requested tolerance
         step_not_OK = .true.
         sc  = 0
         do while ( step_not_OK .and. sc.lt.15  )
            sc = sc + 1

*           Step with curret step size
            call int_step( posc, velc, poss, vels, step, t )
*           Now take two steps with half size
            call int_step( posa, vela, poss, vels, step/2, t )
            call int_step( posb, velb, posa, vela, step/2, t+step/2 )
*           Compare the two results and see if the meet tolerances or if
*           results too good can we increase step
            call eval_step( t, posc, posb, step, new_step, sc)

*           Make sure we have not increased step at wrong time
            if( no_step_increase .and. new_step.gt.step ) then
                new_step = step
            endif

*           If the new_step is the same as step then we are done.
            if( new_step.eq.step ) then
                step_not_OK = .false.
            else
                step = new_step
            endif
         end do
*        save smallest step size
         if( step.lt. minstep ) minstep = step 

****     OK step hase been correctly taken, update state and time
         t = t + step
         poss(:,1:num_body) = posb(:,1:num_body)
         vels(:,1:num_body) = velb(:,1:num_body)

*****    Use more man's graphics with vt100 sequences
         call animate(t, step, poss)
      end do

****  Output the final values
      call report_ICS(t, poss, vels,'Final conditions')
      write(*,320) minstep, closest, close_body
 320  format('Smallest step size needed ',F9.6,' days',/,
     .       'Closest approach distance was ',ES14.5,
     .       ' km Bodies ',I2,' and ',i2)


      end


      subroutine read_runstring( infile, dur, tol )

      implicit none

*     Routine to read the runsting for hw02_03
*     % hw02_03 <in file> <dur> <tol>

      real*8 dur ! Duration of integration days
     .,      tol  ! Relative tolerence of integration

      character*(*) infile  ! Name of input file

* LOCAL VARIABLES

      integer*4 ierr ! IOSTAT errors
     .,         lenline  ! Length of string
      
      character*64 runstring

****  Get the input file name
      call getarg(1,infile)
      if( lenline(infile).eq.0 ) then
           write(*,120)
 120       format('HW02_03:N-body problem.  Correct runstring:',/,
     .            '% HW02_03 <infile> <duration (days)> <tolerace>',/,
     .            'Default duration is 550 days, ',
     .            'default tolerance 1.e-6')
           stop 'HW02_03: No file name passed'
      endif

****  See if other arguments are passed
      call getarg(2,runstring)
      if( lenline(runstring).gt.0 ) then
          read(runstring,*,iostat=ierr) dur
          if( ierr.ne.0 ) then
              write(*,220) ierr, runstring(1:lenline(runstring))
 220          format('IOSTAT error ',i5,' reading duration from ',a)
              stop 'HW02_3: Error reading duraton'
          endif
      else
          dur = 515.d0
      endif

*     Check tolerance
      call getarg(3,runstring)
      if( lenline(runstring).gt.0 ) then
          read(runstring,*,iostat=ierr) tol
          if( ierr.ne.0 ) then
              write(*,240) ierr, runstring(1:lenline(runstring))
 240          format('IOSTAT error ',i5,' reading tolerance from ',a)
              stop 'HW02_3: Error reading tolerance'
          endif
      else
          tol = 1.d-6
      endif

***** Thats all
      return
      end

      subroutine read_nbody

      implicit none

*     Routine to read the input data file.
*     File format lines are 
*     Mass (kg)  Xpos (km)  Ypos (km)  Xvel (km/s)  Yvel (km/s)
*
      include 'NBody.h'

      integer*4 n  ! Counter for body number
     .,         ierr ! IOSTAT error
     .,         jerr ! IOSTAT error reading data line
     .,         lenline  ! Fuction to return length of line

      character*128 line


***** Open the file and readlines
      open(20,file=infile,iostat=ierr,status='old')
      if( ierr.ne. 0 ) then
          call report_error(ierr,
     .        'Opening N-body file ' // infile(1:lenline(infile)))
          stop 'HW02_03: Error opening N-body file'
      endif

***** Read file
      n = 0
      do while ( ierr.eq.0 .and. n.lt.max_body ) 
         read(20,'(a)',iostat=ierr) line
         if( ierr.eq.0 .and. line(1:1).eq.' ' ) then
*            read the line
             n = n + 1
             read(line,*,iostat=jerr) mass(n), posi(:,n), veli(:,n)
             if( jerr.ne.0 ) then
                call report_error(jerr,
     .             'Reading N-body line ' // line(1:lenline(line)))
                stop 'HW02_03: Error reading N-body line'
             endif
         endif
      enddo

      if( n.eq. max_body ) then
          write(*,220) max_body
 220      format('HW02_03: Maximum number of bodies ',i5,' reached')
      endif

      num_body = n

****  Thats all 
      return
      end 
 
             
      subroutine report_ICS(t, posr, velr, type)

      implicit none

*     Routine to report status of bodies at time t ond of type type

      include 'NBody.h'

* PASSED VARIABLES
      real*8 t  ! Current time days
      real*8 posr(2,num_body)   ! Positions of bodies (km)
     .,      velr(2,num_body)   ! Velocities (km/s)
       
      character*(*) type   ! Type of output.: String passed by user

* LOCAL VARIABLES
      integer*4 i


      write(*,120) type, t
 120  format('+ 12.010 HW 02 Q 03: ',a,' At time ',F14.5,' days',/,
     .       'Body     Mass (kg)   PosX (km)     PosY (km)    ',
     .       'VelX  (km/s)  VelY   (km/s)')

      do i = 1, num_body
         write(*,140) i, mass(i), posr(:,i),velr(:,i)
 140     format(I4,2x,ES18.6,2x,2(ES18.7,1x),2x,2(ES18.7,1x))
      end do

      return
      end



      subroutine report_error(ierr,mess)

      implicit none

*     Routine to report IOSTAT errors.  Initially developed
*     for 12.010 HW2 Problem 2.

* PASSED VARIABLES

      integer*4 ierr
      character*(*) mess

* LOCAL VARIABLES
* None
      
***** See if the IOSTAT error was non-zero
      if( ierr.ne.0 ) then
         write(*,120) ierr, mess
 120     format(/,'IOSTAT Error ',i4,' occurred ',a)
      end if

      return
      end
 
      integer*4 function lenline(in)

      implicit none

*     Returns length of used portion of string.  Originally
*     developed for HW2 Problem 2.  f90 has an intrisic called
*     len_trim(<string>) that performs the same function.

* PASSED VARIABLES
      
      character*(*) in   ! The string passed that we need to find
                          ! the last character of

* LOCAL VARIABLES
      integer*4 len_in    ! Declared length of in string
      integer*4 i         ! Counter used to work backwards through string

***** Get the declared length of string
      len_in = LEN(in)
      i = len_in
      do while ( i.gt.0 .and. in(i:i).eq.' ')
          i = i - 1
      end do

      lenline = i

      return
      end


      subroutine animate(t, step, pos)

      implicit none

*     Routine to animate the changes in position

      include 'NBody.h' 

* PASSED Variables
      real*8 t                ! Current time
     .,      step             ! Current step size
     .,      pos(2,num_body)  ! Positions of bodies (also use mass to
                              ! get size
      real*8 maxx, maxy, minx,miny  ! Min and max coordintes in initial state
     .,      maxa, mina  ! Max and min of both x and y

      real*8 scalex, scaley  ! Scales to convert x and y to screen coordinates
     .,      offx, offy      ! Offsets for origin in X and Y

      integer*4 SWX, SHY    ! Screen width and height
     .,         IX, IY      ! Coorinates in screen coordinatea
     .,      i   ! Loop counter

       character*1 sym(max_body)  ! Symbols for each body


       character*85 screen(45)   ! 45-lines of 85-characters to output
       character*6 mv45   ! Escape sequence to move up 51 lines

       save sym, scalex, scaley, offx, offy, mv45

****   Initialize when t = 0;
       if( t .eq. 0 ) then
*          First call; so don't move cursor back up and compute size we need
*          Get o or . based on mass
           do i = 1, num_body
              if( mass(i).gt. 1.e30 ) then
                  sym(i) = 'O'
              elseif( mass(i).gt. 1.e28 ) then
                  sym(i) = 'o'
              else
                  sym(i) = '.'
              endif
           end do
*          Now get the scale
           SWX = 85  ! 85 characeter accross (5/3 ratio)
           SHY = 45  ! 45 lines down 
          
           maxx = pos(1,1)
           minx = pos(1,1)
           maxy = pos(2,1)
           miny = pos(2,1)
           do i = 2,num_body
              maxx = max(maxx, pos(1,i))
              minx = min(minx, pos(1,i))
              maxy = max(maxy, pos(2,i))
              miny = min(miny, pos(2,i))
           enddo

           maxa = max(maxx, maxy)
           mina = min(minx, miny)

*          Set Scale in X and Y
           Scalex = SWX/(maxa-mina) * 0.95d0
           Scaley = SHY/(maxa-mina) * 0.95d0 
           offx = minx* 1.05d0
           offy = (miny - 1.5e8 ) * 1.05d0
           write(mv45,100) char(27)
 100       format(a1,'[47A')

      else
*          Move the cursor back 52-lines
           write(*,'(a)') mv45
      endif
c      print *,'Scale ',scalex, scaley, offx, offy

      write(*,120) t, step
 120  format('TIME ',F13.6,' days, Step ',F10.6, ' days')

*     OK form up the output lines
      do i = 1,45
         screen(i) = ' '
      enddo
      do i = 1,num_body
         IX = nint((pos(1,i)-offx)*scalex)
         IY = 45-nint((pos(2,i)-offy)*scaley)
c        print *,'Pos ', IX,IY, ' ',sym(i)
         if( IX.ge.1 .and. IX.le.85 .and.
     .       IY.ge.1 .and. IY.le.45 ) then
             if( screen(IY)(IX:IX).eq.' ' ) then
                screen(IY)(IX:IX) = sym(i)
             endif
         endif
      enddo

*     write lines
      do i = 1,45
        write(*,'(a)') screen(i)
      enddo

****  Thats all
      return
      end

      subroutine int_step( pose, vele, poss, vels, step, t )

      implicit none

*     Routine to integrate 1 time step of duration step in days

      include 'NBody.h'

*     Passed Variables
      real*8 pose(2,num_body)  ! Position at end of step
     .,      vele(2,num_body)  ! velocity at end of step
     .,      poss(2,num_body)  ! Position at start of step
     .,      vels(2,num_body)  ! Velocity at emd of step
     .,      step   ! Step size in days
     .,      t      ! Current time

****  Local variables
      real*8 SS ! Step size in seconds
     .,    P(2,num_body), V(2,num_body), A(2,num_body)  ! Positon,
                ! velcity and acceleration at intermediate points
     .,    k1(2,num_body), k2(2,num_body), k3(2, num_body) 
     .,    k4(2,num_body)   ! Intermediate integration points in
                ! Runge-Kutta integration.

****  Step in seconds
      SS = step*86400.d0

****  Get the accelerations of all bodies at start
      call get_accel(t, poss, vels, A)
      k1 = SS*A

*     Step system to mid-point
      P = poss + SS*vels/2.d0 + SS*k1/8.d0
      V = vels + k1/2.d0

*     Get acceleration at this new point
      call get_accel(t+step/2, P, V, A)
      k2 = SS*A
      P = poss + SS*vels/2.d0 + SS*k1/8.d0
      V = vels + k2/2.d0

*     Nexy value at mid point
      call get_accel(t+step/2, P, V, A)
      k3 = SS*A
      P = poss + SS*Vels + SS*k3/2.d0
      V = vels + k3

*     Now do final end point
      call get_accel(t+step, P, V, A)
      k4 = SS*A
      pose = poss + SS*(vels + (k1+k2+k3)/6.d0)
      vele = vels + (k1+2*k2+2*k3+k4)/6.d0

****  That all
      return 
      end

      subroutine get_accel(t, P, V, A)

      implicit none

*     Routine to compute accelerations

      include 'NBody.h'

* PASSED VARIABLES
      real*8 t  ! Time in days (not used)
     .,   P(2,num_body), V(2,num_body) ! Position and Velocity (not used)
                ! to compute accelerations
     .,   A(2,num_body)   ! Accerations at each body

* LOCAL VARIABLES

      integer*4 i, j   ! Counters over the bodies.
      real*8 dP(2,num_body)  ! Position difference for each bodu
     .,   R(num_body)   ! Distances bewteen current body and others



****  Loop over bodies
      do i = 1, num_body
         dP(1,:) = P(1,:) - P(1,i)
         dP(2,:) = P(2,:) - P(2,i)

         R = sqrt(dP(1,:)**2+dP(2,:)**2)

***      Now get vector sum of accelerations
         A(:,i) = 0.d0
         do j = 1,num_body
            if ( i.ne.j ) then
               A(:,i) = A(:,i) + G_univ*mass(j)*dP(:,j)/R(j)**3
            endif
        enddo
      enddo

****  Thats all
      return
      end

      subroutine eval_step( t, posc, posb, step, new_step, sc)

      implicit none

*     Routine to see if step size meetds tolerances
  
      include 'NBody.h'

* PASSED VARIABLES
      real*8 t   ! current time in days
     .,      posc(2,num_body)  ! Position with 1-step
     .,      posb(2,num_body)  ! Position with 2-step
     .,      step, new_step    ! Current and new estimate of step size
                               ! new is factor of two different or the same

      integer*4 sc    ! Step size evalution counter.  If two any iterations then
                      ! we stop (so debug out as well)

* LOCAL VARIABLE
      real*8 dPe(2,num_body)    ! Difference in position estimates
     .,      R     ! distance betweeb bodies
     .,      Rmin  ! Smallest difference   
     .,      err   ! Position error in km  
     .,      minerr, maxerr ! Min and max relative error   

      integer*4 i, j  ! Loop counters                  


****  Get the difference and see how low
      dPe = posc - posb
*     See if meets tolerance
      maxerr = 0
      minerr = 1.d20
      do i = 1,num_body
         err = sqrt(dPe(1,i)**2+dPe(2,i)**2)
         Rmin = 1.d20
         do j = 1,num_body
             if( j.ne.i ) then 
                R = sqrt((posb(1,j)-posb(1,i))**2 + 
     .                   (posb(2,j)-posb(2,i))**2)
                if( R .lt.Rmin ) Rmin = R
                if( Rmin.lt.closest ) then
                    closest = Rmin
                    close_body(1) = i
                    close_body(2) = j
                endif
             endif
         enddo 
*        Now see if this smallest or largest error
         if( err/Rmin.lt.minerr) minerr = err/Rmin
         if( err/Rmin.gt.maxerr) maxerr = err/Rmin

      enddo

****  See if we have meet tolerance
      if( sc.gt.11 ) then
         write(*,120) t, sc, step, maxerr, minerr, tol
 120     format('UpateStep Time ',F12.6,1x, i3,' Size ',F9.6,' Errs ',
     .           2E13.3,' Tol ',E13.2)
      endif

      new_step = step
      if( maxerr.lt.tol/40 .and. step.lt. 1.d0 ) then
          new_step = step*2
      endif
      if( maxerr.gt.tol .and. step.gt. 1.d0/2.d0**16 ) then
          new_step = step/2
      endif

****  Thats all
      return
      end



              

         
     
