SHARE
TWEET

Gas Giant pressure distribution

a guest Oct 25th, 2013 74 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. ! "Gas Giant" planetary physics analysis code
  2. !   Start with some central pressure, then this code integrates
  3. !  the system of differential equations to predict the gas pressure
  4. !  and gravitational field throughout the interior of that planet.
  5. !  
  6. !   Intended to produce numbers for the fictional world Virga, with
  7. !  two example cases where the central pressure is 1 atmosphere and
  8. !  then 3 atmospheres. Spatial data is outputted to files.
  9. !   Parameters are also put in place for Jupiter and the sun, but these
  10. !  cases have a problem with spatial resolution near the center.
  11. !  With 10000 spatial points, Jupiter can be resolved to some extent but
  12. !  not the sun.
  13.  
  14. module gas_values
  15.   double precision, parameter :: G = 6.67384e-11
  16.   double precision :: Rsp = 287.058
  17.   double precision :: T = 293.
  18. end module
  19.  
  20. program ball
  21.   use gas_values
  22.   implicit none
  23.   double precision :: r, r_max, r_step, rold
  24.   double precision, dimension(2) :: x, xold
  25.   integer :: theunit
  26.   double precision :: kB, amu
  27.   interface
  28.     function Pg(xp,tp)
  29.       implicit none
  30.       double precision, dimension(2), intent(in) :: xp
  31.       double precision, intent(in) :: tp
  32.       double precision, dimension(2) :: Pg
  33.     end function
  34.   end interface
  35.   interface
  36.     function PM(xp,tp)
  37.       implicit none
  38.       double precision, dimension(2), intent(in) :: xp
  39.       double precision, intent(in) :: tp
  40.       double precision, dimension(2) :: PM
  41.     end function
  42.   end interface
  43.  
  44.   open( unit=2, file="output_1atm.txt", status="replace" )
  45.   open( unit=3, file="output_3atm.txt", status="replace" )
  46.   open( unit=4, file="output_jupiter.txt", status="replace" )
  47.   open( unit=5, file="output_sun.txt", status="replace" )
  48.  
  49.   r_max = 1.0e8
  50.   r_step = r_max / 10000
  51.  
  52.   kB = 1.3806488e-23 ! boltzmann constant
  53.   amu = 1.66053892e-27 ! in kg
  54.   write(*,*) ' R specific for air '
  55.   write(*,*) ' using: ',Rsp
  56.   write(*,*) ' calc: ', kB / ( 28.97 * amu )
  57.  
  58.   ! normal case
  59.   x = (/ 1.0e5, 0.0 /)
  60.   theunit = 2
  61.   r = 100.0
  62.   write(theunit,*) '  radius     pressure     field '
  63.   do
  64.     write(theunit,*) r, x
  65.     xold = x
  66.     call rk4(Pg,xold,x,r,r_step)
  67.     r = r + r_step
  68.     if (r > r_max ) exit
  69.   end do
  70.  
  71.  
  72.   ! high pressure case
  73.   x = (/ 3.0e5, 0.0 /)
  74.   theunit = 3
  75.   r = 100.0
  76.   write(theunit,*) '  radius     pressure     field '
  77.   do
  78.     write(theunit,*) r, x
  79.     xold = x
  80.     call rk4(Pg,xold,x,r,r_step)
  81.     r = r + r_step
  82.     if (r > r_max ) exit
  83.   end do
  84.  
  85.   ! Jupiter case
  86.   Rsp = kB / ( (1.*0.898+4.*0.102+16.04*0.003)* amu )
  87.   write(*,*) ' jupiter Rsp ',Rsp
  88.   x = (/ 1.0e5*100.0e6, 0.0 /) ! 6 million atmospheres at center
  89.   theunit = 4
  90.   r = 100.0
  91.   write(theunit,*) '  radius     pressure     field '
  92.   do
  93.     write(theunit,*) r, x
  94.     xold = x
  95.     call rk4(Pg,xold,x,r,r_step)
  96.     r = r + r_step
  97.     if (r > r_max ) exit
  98.   end do
  99.  
  100.   ! Sun case
  101.   Rsp = kB / ( (1.*0.912+4.*0.087+16.*0.0097)* amu )
  102.   write(*,*) ' sun Rsp ',Rsp
  103.   x = (/ 1.0e5*2.5e11, 0.0 /)
  104.   T = 5500.
  105.   theunit = 5
  106.   r = 100.0
  107.   write(theunit,*) '  radius     pressure     field '
  108.   do
  109.     write(theunit,*) r, x
  110.     xold = x
  111.     call rk4(Pg,xold,x,r,r_step)
  112.     r = r + r_step
  113.     if (r > r_max ) exit
  114.   end do
  115.  
  116.   close(2)
  117.   close(3)
  118.   close(4)
  119.   close(5)
  120.  
  121. end program ball
  122.  
  123.  
  124. function PM(x,r)
  125.   use gas_values
  126.   implicit none
  127.   double precision, dimension(2), intent(in) :: x
  128.   double precision, intent(in) :: r
  129.   double precision, dimension(2) :: PM
  130.  
  131.   PM(1) = (G / (Rsp * T))*(r**2)*x(2)
  132.   PM(2) = -(4. * 3.141592654 / (Rsp * T))*x(1)*x(2)/(r**2)
  133. end function
  134.  
  135.  
  136. function Pg(x,r)
  137.   use gas_values
  138.   implicit none
  139.   double precision, dimension(2), intent(in) :: x
  140.   double precision, intent(in) :: r
  141.   double precision, dimension(2) :: Pg
  142.    
  143.   Pg(1) = -(G / (Rsp * T))*x(2)*x(1)/G
  144.   Pg(2) = (4. * 3.141592654 / (Rsp * T))*G*x(1) - 2*x(2) / r
  145. end function
  146.  
  147. subroutine rk4(F,x0,x,t0,delt)
  148.   implicit none
  149.   double precision, intent(in) :: x0(2)
  150.   double precision, intent(out) :: x(2)
  151.   double precision, intent(in) :: t0, delt
  152.   double precision :: k1(2), k2(2), k3(2), k4(2)
  153.   interface
  154.     function F(xp,tp)
  155.       implicit none
  156.       double precision, dimension(2), intent(in) :: xp
  157.       double precision, intent(in) :: tp
  158.       double precision, dimension(2) :: F
  159.     end function
  160.   end interface
  161.  
  162.   k1 = F(x0,t0)
  163.   k2 = F(x0 + k1*delt/(2.0),t0+0.5*delt)
  164.   k3 = F(x0 + k2*delt/(2.0),t0+0.5*delt)
  165.   k4 = F(x0 + k3*delt, t0 + delt)
  166. !   write(*,*) ' k ',k1,k2,k3,k4
  167.   x = x0 + (k1+2*k2+2*k3+k4)*delt/(6.0)
  168.  
  169. end subroutine rk4
RAW Paste Data
Pastebin PRO Summer Special!
Get 40% OFF on Pastebin PRO accounts!
Top