!        Galactic Travel via Hohmann transfers
!     Given the starting location and destinations
!    this program will calculate the transfer time,
!     orbital transfer speeds, and plot the orbit 
!       using the concepts of Hohmann transfers
!  assuming all planets have planar, circular orbits

!                 Created by 
!                Derek Schutz
!                2 - 28 - 07

!jc  Good header comment.

PROGRAM project
IMPLICIT NONE

!jc  More comments in the variable declarations are needed.
!jc  Explain what you variables are, as well as their units.

  CHARACTER (len=10)::destination, departure
  INTEGER :: i, j, k, m, n, p, h, IER, PGBEG, numpts, ypt
  REAL (KIND=8) :: km, a, b, t, d, v1, v2, delv1, delv2, vt1, vt2, time(1:9,1:9), e, alpha, semiminor_axis,&
 ypts(-100:100), xpts(-100:100), az, bz, negypts(-100:100)

  TYPE :: planet
    REAL (KIND=8) :: distance_to_sun
    CHARACTER(len=7) :: name
    INTEGER :: name_length
  END TYPE

!jc  Nice use of structures.

  TYPE ( planet ) :: planets(9) !mercury, venus, earth, mars, jupiter, saturn, uranus, neptune, pluto
    planets(1) = planet(5.74099E10, "mercury", 7)
    planets(2) = planet(1.09209E11, 'venus', 5)
    planets(3) = planet( 1.49598E11 , 'earth', 5)
    planets(4) = planet(2.27937E11, 'mars', 4)
    planets(5) = planet(7.7857E11, 'jupiter', 7)
    planets(6) = planet(1.42673E12, 'saturn', 6)
    planets(7) = planet(2.87097E12, 'uranus', 6)
    planets(8) = planet(4.498E12, 'neptune', 7)
    planets(9) = planet(5.906E12, 'pluto', 6)

 PRINT * , ''
 PRINT * , 'Days for transfers'					!Generates a table for the dayseach transfer takes
 PRINT * ,'      |Mercury|   Venus|   Earth|    Mars| Jupiter|  Saturn|  Uranus| Neptune|   Pluto|'
 PRINT * ,'------|--------------------------------------------------------------------------------'

  DO m = 1, 9
    DO n=1, 9

      IF (m==n) THEN
        time(n,m)=0.
      ELSE IF(n/=m) THEN
        km = 7.53E-21 !m/k in s^2/m^3
        time(n,m) = 3.14159265*SQRT((( planets(n)%distance_to_sun + planets(m)%distance_to_sun )**3.)*(km)/8.)/86400.
      END IF

    END DO

  END DO

  DO p = 1, 9
    WRITE (*,150) planets(p)%name, time(1:9,p)
    150 FORMAT ('',(A), '|',9(F7.1'| '))
  END DO

  PRINT * , ''
  PRINT * , 'Destination and Departure choices'

  DO i=1,9
    PRINT * , i, planets(i)%name
  END DO

  PRINT * , 'Where would you like to go?'
  READ * , j


!jc  This should really be in a loop, so if they put in a bad choice, they
!jc  have to continue entering it till they get it correct.

  IF ( j < 1 ) THEN			!If statement takes out possibilities that are not choices for planets
    PRINT * , 'This is not an option, please select a different one' 
    READ * , j
  ELSE IF ( j > 9) THEN
    PRINT * , 'This is not an option, please select a different one'
    READ * , j
  END IF

  WRITE (*,100) planets(j)%name
  100 FORMAT('So you want to go to ', (A))

  PRINT * , 'Where are you going to be leaving from?'
  READ * , k
  PRINT * , ''

  IF ( k < 1 ) THEN   			!If statement takes out possibilities that are not choices for planets
    PRINT * , 'This is not an option, please select a different one'
    READ * , k
  ELSE IF ( k > 9 ) THEN
    PRINT * , 'This is not an option, please select a different one'
    READ * , k
  ELSE IF ( k == j ) THEN
    PRINT * , 'You picked the same planets, please select a different one'
    READ * , k
  END IF

  WRITE (*,110) planets(k)%name, planets(j)%name

  110 FORMAT('Going from ',(A), 'to ',(A),' will take:')


  a = .5*( planets(k)%distance_to_sun + planets(j)%distance_to_sun ) !semimajor axis of the transfer in meters
  t = 3.14159265*( km )**.5 *( a )**1.5                      	    !gives the time the transfer will take in seconds
  d=t/60./60./24.							  !converts seconds into days

  WRITE(*,120) d
  120 FORMAT( (F6.0), ' days to get there')

  PRINT * , ''

								!calculates the orbital transfer speeds for
  v1 = SQRT(1./(km*planets(k)%distance_to_sun))			!departure and return
  vt1 = SQRT(2.*((planets(j)%distance_to_sun)/( planets(k)%distance_to_sun + planets(j)%distance_to_sun )))*v1

  delv1 = (vt1-v1)/1000. 		!delv1 gives the difference between the planet's orbit speed and the ships transfer speed to get into the orbit to the other planet


  v2 = SQRT(1./(km*planets(j)%distance_to_sun))
  vt2 = SQRT(2.*((planets(k)%distance_to_sun)/( planets(k)%distance_to_sun + planets(j)%distance_to_sun )))*v2

  delv2 = (vt2 - v2)/1000.	!delv1 gives the difference between the planet's orbit speed and the ships transfer speed to get into the orbit to the other planet


!jc  Calling your delta v, "speed transfer" is a bit misleading.

  WRITE (*, 130) planets(k)%name, delv1
  130 FORMAT('The speed transfer from ', (A), '  is ', (F6.1), ' km/s')	

  WRITE (*, 140) planets(k)%name, delv2
  140 FORMAT('The speed transfer to ', (A), '  is ', (F6.1), ' km/s')

  IF (planets(j)%distance_to_sun < planets(k)%distance_to_sun) THEN

	e = 1. - (planets(j)%distance_to_sun) / a
  ELSE IF (planets(k)%distance_to_sun < planets(j)%distance_to_sun) THEN

	e = 1. - (planets(k)%distance_to_sun) / a
  END IF
PRINT *, 'The eccentricity of the transfer orbit is ', e


  IF (planets(j)%distance_to_sun < planets(k)%distance_to_sun) THEN	!Calculate alpha for finding the semiminor axis

	alpha = planets(j)%distance_to_sun*(1. + e)
  ELSE IF (planets(k)%distance_to_sun < planets(j)%distance_to_sun) THEN

	alpha = planets(k)%distance_to_sun*(1. + e)
  END IF

  PRINT *, 'The semimajor axis for this transfer orbit is ', a

  semiminor_axis = SQRT(alpha*a)
  PRINT * , 'The semiminor axis for this transfer orbit is ', semiminor_axis

  bz=semiminor_axis/100000000.					!Puts semiminor axis into a workable value
  az = a /100000000.						!Puts semimajor axis into a workable value
  PRINT *, bz, az

  DO h = -100, 100, 1						!calculate points for the pgplot based on 
								!semimajor axis
	xpts(h) = az*h/100.
	ypts(h) = (bz)*SQRT(ABS(1. - ((az*h/100.)**2 / (az)**2)))
	negypts(h) = -(bz)*SQRT(ABS(1. - ((az*h/100.)**2 / (az)**2)))

  END DO

  PRINT * , ''

  IER = PGBEG(0, 'proj3.ps/PS', 1, 1)		!name postscript file

  IF (IER.NE.1) STOP

!jc  For these plots you should have used the option that forces the scale
!jc   to be the same on both axes.  That way it would have been more
!jc  visibily obvious how the shapes of the orbits change.  Right now 
!jc  the plots all look the same, it is just the scales that change.
!jc  You also might have just plotted one half of the orbit, since that
!jc  is the usual path used for transfer.

  CALL PGENV (REAL(-az), REAL(az), REAL(-bz), REAL(bz), 0, 1)	!Call subroutines for PGPlot
  CALL PGLAB ('1 x 10^8 m','1 x 10^8 m', 'Hohmann Transfer Orbit')		
  CALL PGPT (201, REAL(xpts(-100:100:1)), REAL(ypts(-100:100:1)), 20)
  CALL PGLINE (201, REAL(xpts(-100:100:1)), REAL(ypts(-100:100:1)))
  CALL PGPT (201, REAL(xpts(-100:100:1)), REAL(negypts(-100:100:1)), 20)
  CALL PGLINE (201, REAL(xpts(-100:100:1)), REAL(negypts(-100:100:1)))
  !jc  You don't need to include the step size.
  CALL PGEND

STOP

END PROGRAM

!jc  OK program.  The program could have used a few more comments,
!jc  as well as a bit more complexity.