Starfire Kickstarter Calculations

! Starfire kickstarter claculations (suborbital space cannon)
!   Does simulation of trajectory including air resistance
!   Assumes gun fired directly up
!   Does not include wind
!   Several specific cases inputed  &  mass vs. altitude table produced
!
! This calculates full version of calculations outlined here:
!   http://space.stackexchange.com/questions/3667/could-any-existing-gun-reach-the-karman-line
!
! rk4 approach inspired by
! http://www.mathworks.com/matlabcentral/fileexchange/2181-numerical-methods-using-matlab-2e/content/edition2/matlab/chap_9/rk4.m
!
! Pressure/density equation:
!  P(r) = P0 exp(H' (1/r - 1/r0) )
!   H' = r0^2 / H
!   H = 8.6 km     or maybe    7.4 km  ...depends on temperature
!   r0 = 6,378.1 km

module problem_values
  implicit none
  real, parameter :: Hchar = 7.4, &
                    & r0 = 6378.1, &
                    & pi = 3.14159265359, &
                    & tube_length = 13.716, & ! 45 feet
!                     & tube_length = 19.812, & ! 65 feet
                    & tube_diam = 0.2032, & ! 8 inches
                    & tube_P = 4.13685e8 ! Pascals

  ! case specific
  real :: v0, diam, mass, Cd, area, rho0
  logical :: toprint

end module problem_values

program rk4_shell
  use problem_values
  implicit none
  integer, parameter :: N = 2
  real, dimension(2) :: xp_temp, x_temp, x0
  interface
    function x_prime(x,t)
      implicit none
      real, dimension(2) :: x_prime
      real, dimension(2), intent(in) :: x
      real, intent(in) :: t
    end function
  end interface

  toprint = .false.  ! turn back on to print the trajectory to file
                     ! but comment out the cases you don't want, only prints once

  write(*,*) '     General calculations: '
  write(*,*) ' maximum KE with tube = ',(tube_diam/2.)**2*tube_length*tube_P,' J '
  write(*,*) '     = ',(tube_diam/2.)**2*tube_length*tube_P/(1.e6),' MJ '
  write(*,*) ' '

  !  case-specific data

  ! tank gun data
  v0 = 1750.0
  diam = 0.12
  mass = 21.0
  Cd = 0.47
  rho0 = 1.3

  ! starfire basic
  v0 = 7010.0
  diam = 0.2032
  mass = 34.05
  Cd = 0.47
  rho0 = 1.0

  ! starfire with revised missile shape, high altitude
  v0 = 1500.
  diam = 0.2032 / 2.
  mass = 65.9 ! (2 inches)^2*Pi*(8000 kg/m^3)*(40 inches)
  Cd = 0.75
  rho0 = 0.447*1.3
      ! image is 38 px high, 596 long, aspect ratio of 13

  ! maximum allowable mass
  mass = 163.5
  diam = 0.2032
  Cd = 0.5

  ! generous Cd value
  rho0 = 1.3
  Cd = 0.1

  write(*,*) '    Parameters: '
  write(*,*) ' mass            = ',mass,' kg'
  write(*,*) ' muzzle velocity = ',v0,' m/s'
  write(*,*) ' tube diameter   = ',diam,' meters '
  write(*,*) ' drag coefficient= ',Cd,' unitless '
  write(*,*) ' air density     = ',rho0,' kg/m3 '
  write(*,*) ' altitude        = ',-log(rho0/1.3)*Hchar + 228.,' m'
  write(*,*) ' '

  area = (diam/2.)**2*pi

  x0(1) = 0.   ! state vector
  x0(2) = v0
  ! composed of  ( position , velocity )

  write(*,'(a)',advance='no') '  Maximum altitude (km) = '
  call rk4(x_prime,N,0.,1000.,x0,10000,x_temp)
  write(*,*) ' '
  write(*,*) ' '

  write(*,*) '  Iteration on mass values with calculated muzzle velocity '
  diam = 0.2032
  rho0 = 1.3
  Cd = 0.3

  mass = 5.0
  write(*,*) ' '
  write(*,*) '     M         v0       max_height '
  do
    v0 = sqrt( 2. * (tube_diam/2.)**2*tube_length*tube_P / mass )
    x0(1) = 0.   ! state vector
    x0(2) = v0
    write(*,'(2F15.5)',advance='no') mass, v0
    call rk4(x_prime,N,0.,1000.,x0,10000,x_temp)
    if (mass > 1000.) exit
    mass = mass * 1.2
  end do


end program rk4_shell

function x_prime(x,t)
  use problem_values
  implicit none
  real, dimension(2) :: x_prime
  real, dimension(2), intent(in) :: x
  real, intent(in) :: t
  real :: Hp, rho

  Hp = r0**2 / Hchar
  rho = rho0*exp(Hp*(1/(r0+x(1)/1000.)-1/r0))  ! equation for density with height
  !   otherwise, rho = rho0 * exp( -h/H )
!   write(*,*) ' (h,rho) ',x(1),rho

  x_prime(1) = x(2)
  x_prime(2) = -0.5*Cd*x(2)*abs(x(2))*area*rho / mass - 9.8

end function x_prime

subroutine rk4(f,n,a,b,x0,m,x)
  use problem_values
  implicit none
  integer, intent(in) :: n,m   ! n dim of vector,  m steps to use
  interface
    function f(x,t)
      implicit none
      real, dimension(2) :: f
      real, dimension(2), intent(in) :: x
      real, intent(in) :: t
    end function
  end interface
  real, intent(in) :: a,b
  real, dimension(2), intent(in) :: x0
  real, dimension(2), intent(out) :: x
  real :: h, t, h_max
  integer :: j
  real, dimension(n) :: k1,k2,k3,k4

  h = (b-a)/m
  h_max = h
  t = a
  x = x0
  if (toprint) then
    open( unit = 29 , file= "trajectory.dat" , status = "replace" )
    write(29,*) '   altitude     velocity'
  end if
  do j = 1, m
    if (toprint) then
      if (mod(j,10)==0) then
        write(29,*) x
      end if
    end if
    k1 = h*f(x,t)
    k2 = h*f(x+k1/2.,t+h/2.)
    k3 = h*f(x+k2/2.,t+h/2.)
    k4 = h*f(x+k3,t+h)
    t = t + h
    x = x + (k1+2*k2+2*k3+k4)/6.
    if (x(1) > h_max) h_max = x(1)
    if (x(1) < 0.) exit
  end do
!   write(*,*) ' '
!   write(*,*) ' max_h ',h_max,' m'
!   write(*,*) '       ',h_max/1000.,' km '
  write(*,*) h_max/1000.
  if (toprint) close(29)


end subroutine rk4