equation_peos.f90

Go to the documentation of this file.
00001 ! --- tov equations are in this routine.
00002 !
00003 SUBROUTINE equation(x, y, f)
00004   use phys_constant, only : pi
00005   implicit none
00006 !
00007   integer, parameter :: neq = 6
00008   real(8), intent(inout) :: y(neq), f(neq), x
00009   real(8) :: rr, rma, hh, bma, tma, psi, adm, 
00010             rho0, pre, ene, q
00011 !
00012   rr  = x
00013   rma = y(1)
00014   hh  = y(2)
00015   bma = y(3)
00016   tma = y(4)
00017   psi = y(5)
00018   adm = y(6)
00019 !  if (q <= 0.0d0) then
00020   if (hh <= 1.0d0) then
00021     q = 0.0d0 !!!
00022     pre = 0.0d0
00023     rho0= 0.0d0
00024     ene = 0.0d0
00025   else
00026     call peos_h2qprho(hh, q, pre, rho0, ene)
00027   end if
00028 !
00029   if (rr == 0.0d0) then
00030     f(1:6) = 0.0d0
00031   else
00032     f(1) = 4.0d0*pi*rr**2*ene
00033 !c      f(2) = - (hh + 0.93d0)*(rma + 4.0d0*pi*pre*rr**3)
00034 !c      f(2) = - (hh + 1.0d0)*(rma + 4.0d0*pi*pre*rr**3)
00035     f(2) = - hh*(rma + 4.0d0*pi*pre*rr**3) &
00036       &    /(rr**2 - 2.0d0*rma*rr)
00037     f(3) = 4.0d0*pi*rr**2*rho0/(1.0d0-2.0d0*rma/rr)**0.5d0
00038     f(4) = 4.0d0*pi*rr**2*ene/(1.0d0-2.0d0*rma/rr)**0.5d0
00039     f(5) = 0.5d0*psi/rr*(1.0d0-1.0d0/(1.0d0-2.0d0*rma/rr)**0.5d0)
00040     f(6) = 4.0d0*pi*rr**2*ene/(1.0d0-2.0d0*rma/rr)**0.5d0/psi
00041   end if
00042 !
00043   return
00044 END SUBROUTINE equation