Guest User

Gas Giant pressure distribution

a guest
Oct 25th, 2013
125
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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×