equation_polytrope.f90

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