!!****if* source/Simulation/SimulationMain/Test/Simulation_init
!!
!! NAME
!!
!!  Simulation_init
!!
!!
!! SYNOPSIS
!!
!!  Simulation_init()
!!
!!
!! DESCRIPTION
!!
!!  Initializes all the parameters needed for the Sod shock tube
!!  problem
!!
!! ARGUMENTS
!!
!!  
!!
!! PARAMETERS
!!
!!  sim_rhoTarg    Density in the left part of the grid
!!  sim_rhoAtm   Density in the right part of the grid
!!  sim_pTarg      Pressure  in the left part of the grid
!!  sim_pAtm     Pressure  in the righ part of the grid
!!  sim_uTarg      fluid velocity in the left part of the grid
!!  sim_uAtm     fluid velocity in the right part of the grid
!!  sim_xangle     Angle made by diaphragm normal w/x-axis (deg)
!!  sim_ yangle    Angle made by diaphragm normal w/y-axis (deg)
!!  sim_posSurf     Point of intersection between the shock plane and the x-axis
!!  sim_rhoObj     Density of impactor
!!  sim_radiusObj  Radius of impactor
!!
!!  left is bottom, right is top
!!***

subroutine Simulation_init()
  
  use Simulation_data
  use Driver_interface, ONLY : Driver_getMype, Driver_abortFlash
  use RuntimeParameters_interface, ONLY : RuntimeParameters_get
  use Logfile_interface, ONLY : Logfile_stamp
  use Multispecies_interface, ONLY : Multispecies_setProperty
  implicit none
#include "constants.h"
#include "Flash.h"
#include "Multispecies.h"

  
  call Driver_getMype(MESH_COMM, sim_meshMe)


  call RuntimeParameters_get('smallp', sim_smallP)
  call RuntimeParameters_get('smallx', sim_smallX) 
  
  call RuntimeParameters_get('gamma', sim_gamma)
  
  call RuntimeParameters_get('sim_rhoTarg', sim_rhoTarg)
  call RuntimeParameters_get('sim_rhoAtm', sim_rhoAtm)
  
  call RuntimeParameters_get('sim_pTarg', sim_pTarg)
  call RuntimeParameters_get('sim_pAtm', sim_pAtm)
  
  call RuntimeParameters_get('sim_uTarg', sim_uTarg)
  call RuntimeParameters_get('sim_uAtm', sim_uAtm)
  
  call RuntimeParameters_get('sim_xangle', sim_xAngle)
  call RuntimeParameters_get('sim_yangle', sim_yAngle)
  
  call RuntimeParameters_get('sim_posSurf', sim_posSurf)
  call RuntimeParameters_get('sim_rhoObj', sim_rhoObj)
  call RuntimeParameters_get('sim_radiusObj', sim_radiusObj)
  call RuntimeParameters_get('sim_xObj', sim_xObj)
  call RuntimeParameters_get('sim_yObj', sim_yObj)
  call RuntimeParameters_get('sim_uObj', sim_uObj)

#ifdef SIMULATION_TWO_MATERIALS
  call RuntimeParameters_get('sim_abarTarg', sim_abarTarg)
  call RuntimeParameters_get('sim_zbarTarg', sim_zbarTarg)
  call RuntimeParameters_get('sim_abarAtm', sim_abarAtm)
  call RuntimeParameters_get('sim_zbarAtm', sim_zbarAtm)
#endif

  call Logfile_stamp( "initializing Sod problem",  &
       "[Simulation_init]")
     
#ifdef SIMULATION_TWO_MATERIALS
  call Multispecies_setProperty(LEFT_SPEC, A, sim_abarTarg)
  call Multispecies_setProperty(LEFT_SPEC, Z, sim_zbarTarg)
  call Multispecies_setProperty(RGHT_SPEC, A, sim_abarAtm)
  call Multispecies_setProperty(RGHT_SPEC, Z, sim_zbarAtm)
#endif

#ifdef FLASH_3T
  call RuntimeParameters_get('sim_pionTarg' , sim_pionTarg )
  call RuntimeParameters_get('sim_pionAtm', sim_pionAtm)
  call RuntimeParameters_get('sim_peleTarg' , sim_peleTarg )
  call RuntimeParameters_get('sim_peleAtm', sim_peleAtm)
  call RuntimeParameters_get('sim_pradTarg' , sim_pradTarg )
  call RuntimeParameters_get('sim_pradAtm', sim_pradAtm)

  if (sim_pionTarg  < 0.0) call Driver_abortFlash("Must specify sim_pionTarg" )
  if (sim_pionAtm < 0.0) call Driver_abortFlash("Must specify sim_pionAtm")
  if (sim_peleTarg  < 0.0) call Driver_abortFlash("Must specify sim_peleTarg" )
  if (sim_peleAtm < 0.0) call Driver_abortFlash("Must specify sim_peleAtm")
  if (sim_pradTarg  < 0.0) call Driver_abortFlash("Must specify sim_pradTarg" )
  if (sim_pradAtm < 0.0) call Driver_abortFlash("Must specify sim_pradAtm")
  
  call RuntimeParameters_get('gammaEle', sim_gammaEle)
  call RuntimeParameters_get('gammaIon', sim_gammaIon) !This may have to be 5./3 when using multiTemp/Gamma - KW
#endif

  ! convert the shock angle paramters
  sim_xAngle = sim_xAngle * 0.0174532925 ! Convert to radians.
  sim_yAngle = sim_yAngle * 0.0174532925

  sim_xCos = cos(sim_xAngle)
  
  if (NDIM == 1) then
     sim_xCos = 1.
     sim_yCos = 0.
     sim_zCos = 0.
     
  elseif (NDIM == 2) then
     sim_yCos = sqrt(1. - sim_xCos*sim_xCos)
     sim_zCos = 0.
     
  elseif (NDIM == 3) then
     sim_yCos = cos(sim_yAngle)
     sim_zCos = sqrt( max(0., 1. - sim_xCos*sim_xCos - sim_yCos*sim_yCos) )
  endif
end subroutine Simulation_init