c Mahc_Init_Data: Various Initial Data Routines for Mahc_Evolve c Copyright (C) 2001 Mark Miller /*@@ @file tov_isotropic.F @date December 1999 @author Mark Miller @desc Solve the 1D TOV equations, with various utilities @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_DefineThorn.h" /*@@ @routine solvetov_iso @date December 1999 @author Mark Miller @desc This routine solves the TOV equations for a given initial central pressure in isotropic coordinates. @enddesc @calls @calledby @history @endhistory @@*/ subroutine solvetov_iso(press,rho,eps,m,phi,r_schwar,r_iso,ntov, . imaxtov,centdenloc,gamma,eosk,rmaxtov,mahc_init_eos_handle) implicit none integer ntov,imaxtov,mahc_init_eos_handle CCTK_REAL press(ntov),rho(ntov),eps(ntov),m(ntov),phi(ntov) CCTK_REAL r_schwar(ntov),r_iso(ntov) CCTK_REAL centdenloc,gamma,eosk,rmaxtov DECLARE_CCTK_PARAMETERS c LOCAL VARS integer i,j,irad,first,CCTK_Equals CCTK_REAL dr CCTK_REAL temp,term1,term2,term3,f,fp CCTK_REAL rho0,rho1,r1,pi,h CCTK_REAL dpress1,dpress2,dpress3,dpress4,dm1,dm2,dm3,dm4 CCTK_REAL mrest_old,mrest_new,dmrest1,dmrest2,dmrest3,dmrest4 CCTK_REAL drschwar1,drschwar2,drschwar3,drschwar4 CCTK_REAL dphi1,dphi2,dphi3,dphi4,s_press,s_m,s_mrest,s_phi,s_rschwar CCTK_REAL presssend,msend,rschwarsend,risosend,riso_scale,phi_const CCTK_REAL rest_mass, energy_mass, proper_radius logical mahc_init_eos_tables #ifdef CACTUSEOS_EOS_1D_BASE #include "EOS_1D_Base.inc" #endif pi = 4.0d0 * atan(1.0d0) mahc_init_eos_tables = CCTK_EQUALS(mahc_init_eos,'tables') first = 0 dr = rmaxtov/dble(ntov-1) write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' write(*,*) '*/*/ solvetov: initial isotropic dr is ',dr write(*,*) '*/*/ central rest mass density is ',centdenloc do i = 1,ntov r_iso(i) = (rmaxtov/dble(ntov-1)) * dble(i-1) enddo c if central density is 0, just return flat space with very c small density if(centdenloc .le. 0.0d0) then write(*,*) '*/*/ Central density is 0: setting TOV spacetime to' write(*,*) '*/*/ be flat with very small density' imaxtov = 1 do i=1,ntov rho(i) = 1.0e-15 press(i) = 0.0d0 m(i) = 0.0d0 phi(i) = 0.0d0 enddo return endif if(mahc_init_eos_tables) then #ifdef CACTUSEOS_EOS_1D_BASE rho(1) = centdenloc press(1) = EOS_1D_pressfromrho(mahc_init_eos_handle,centdenloc) eps(1) = EOS_1D_epsfromrho(mahc_init_eos_handle,centdenloc) #else write(*,*) 'if you want to use tables in tov solvers, you' write(*,*) 'need to compile with EOS_1D_Base' call CCTK_WARN(0,"tov_isotropic") #endif else rho(1) = centdenloc press(1) = eosk * (centdenloc ** gamma) eps(1) = press(1)/((gamma - 1.0d0)*rho(1)) endif m(1) = 0.0d0 mrest_old = 0.0d0 mrest_new = 0.0d0 s_mrest = 0.0d0 phi(1) = 0.0d0 r_schwar(1) = 0.0d0 write(*,*) '*/*/ tov_rk: integrating the TOV equations: ' c integrate the tov equations out do i= 1,ntov-1 h = r_iso(i+1)-r_iso(i) call press_source_iso(press(i),m(i),r_schwar(i),r_iso(i), . eosk,gamma,pi,s_press,mahc_init_eos_tables,mahc_init_eos_handle) call m_source_iso(press(i),m(i),r_schwar(i),r_iso(i), . eosk,gamma,pi,s_m,mahc_init_eos_tables,mahc_init_eos_handle) call mrest_source_iso(press(i),m(i),r_schwar(i),r_iso(i), . eosk,gamma,pi,s_mrest,mahc_init_eos_tables,mahc_init_eos_handle) call phi_source_iso(press(i),m(i),r_schwar(i),r_iso(i), . eosk,gamma,pi,s_phi) call rschwar_source_iso(m(i),r_schwar(i),r_iso(i),s_rschwar) dpress1 = h * s_press dm1 = h * s_m dmrest1 = h * s_mrest dphi1 = h * s_phi drschwar1 = h * s_rschwar presssend = press(i)+dpress1 / 2.0d0 msend = m(i)+dm1 / 2.0d0 rschwarsend = r_schwar(i) + drschwar1/2.0d0 risosend = r_iso(i)+h/2.0d0 call press_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_press,mahc_init_eos_tables,mahc_init_eos_handle) call m_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_m,,mahc_init_eos_tables,mahc_init_eos_handle) call mrest_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_mrest,mahc_init_eos_tables,mahc_init_eos_handle) call phi_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_phi) call rschwar_source_iso(msend,rschwarsend,risosend,s_rschwar) dpress2 = h * s_press dm2 = h * s_m dmrest2 = h * s_mrest dphi2 = h * s_phi drschwar2 = h * s_rschwar presssend = press(i)+dpress2 / 2.0d0 msend = m(i)+dm2 / 2.0d0 rschwarsend = r_schwar(i) + drschwar2/2.0d0 risosend = r_iso(i)+h/2.0d0 call press_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_press,mahc_init_eos_tables,mahc_init_eos_handle) call m_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_m,mahc_init_eos_tables,mahc_init_eos_handle) call mrest_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_mrest,mahc_init_eos_tables,mahc_init_eos_handle) call phi_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_phi) call rschwar_source_iso(msend,rschwarsend,risosend,s_rschwar) dpress3 = h * s_press dm3 = h * s_m dmrest3 = h * s_mrest dphi3 = h * s_phi drschwar3 = h * s_rschwar presssend = press(i)+dpress3 msend = m(i)+dm3 rschwarsend = r_schwar(i) + drschwar3 risosend = r_iso(i)+h call press_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_press,mahc_init_eos_tables,mahc_init_eos_handle) call m_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_m,mahc_init_eos_tables,mahc_init_eos_handle) call mrest_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_mrest,mahc_init_eos_tables,mahc_init_eos_handle) call phi_source_iso(presssend,msend,rschwarsend,risosend, . eosk,gamma,pi,s_phi) call rschwar_source_iso(msend,rschwarsend,risosend,s_rschwar) dpress4 = h * s_press dm4 = h * s_m dmrest4 = h * s_mrest dphi4 = h * s_phi drschwar4 = h * s_rschwar press(i+1) = press(i) + ((dpress2 + dpress3)/3.0d0 + . (dpress1 + dpress4)/6.0d0) m(i+1) = m(i) + ((dm2+dm3)/3.0d0 + (dm1+dm4)/6.0d0) mrest_old = mrest_new mrest_new = mrest_old + ((dmrest2+dmrest3)/3.0d0 + . (dmrest1+dmrest4)/6.0d0) phi(i+1) = phi(i) +((dphi2+dphi3)/3.0d0 + (dphi1+dphi4)/6.0d0) r_schwar(i+1) = r_schwar(i) + ((drschwar2 + drschwar3)/3.0d0 + . (drschwar1 + drschwar4)/6.0d0) if( press(i+1) .le. 0.0d0) then if(first .eq. 0) then first = 1 imaxtov = i press(i+1) = 0.0d0 rho(i+1) = 0.0d0 eps(i+1) = 0.0d0 m(i+1) = m(i) write(*,*) '/*/* tov: ADM mass of star: ',m(imaxtov) c That aint the rest mass, stupid! c write(*,*) '/*/* tov: rest mass of star: ',mrest_old endif endif if(mahc_init_eos_tables) then #ifdef CACTUSEOS_EOS_1D_BASE rho(i+1) = EOS_1D_rhofrompress(mahc_init_eos_handle,press(i+1)) eps(i+1) = EOS_1D_epsfrompress(mahc_init_eos_handle,press(i+1)) #endif else rho(i+1) = (press(i+1)/eosk)**(1.0d0/gamma) if( rho(i+1) .le. 0.0d0) then eps(i+1) = 0.0d0 else eps(i+1) = press(i+1)/((gamma-1.0d0)*rho(i+1)) endif endif if(first .eq. 1) exit if(first .eq. 1) then write(*,*) 'tov: I am not exiting!!!!' endif enddo if(first .eq. 0) then write(*,*) 'solvetov_iso: you ran out of 1-d array for the tov' write(*,*) ' solve. Make rmaxtov bigger and try again!' STOP endif if(imaxtov + 2 .ge. ntov) then write(*,*) 'solvetov_iso: you ran out of 1-d array for the tov' write(*,*) ' solve. Make rmaxtov bigger and try again!' STOP endif c match the solution to schwarzschild exterior c first, scale the isotropic radial coordinate to match vacuum c metric riso_scale = (r_schwar(imaxtov) - m(imaxtov) + . sqrt(r_schwar(imaxtov) * (r_schwar(imaxtov) - 2.0d0*m(imaxtov))))/ . (2.0d0 * r_iso(imaxtov)) do i=1,ntov r_iso(i) = riso_scale * r_iso(i) enddo write(*,*) '/*/* tov_iso: m/(2 r_iso) at stellar surfaces is ', . m(imaxtov)/(2.0d0 * r_iso(imaxtov)) c now, add const to phi to match the lapse at the surface phi_const = 0.5d0*log(((1.0d0 - m(imaxtov)/(2.0d0*r_iso(imaxtov)))**2)/ . ((1.0d0 + m(imaxtov)/(2.0d0*r_iso(imaxtov)))**2))- . phi(imaxtov) do i=1,imaxtov phi(i) = phi_const + phi(i) enddo c fix up phi at (imaxtov+1) and (imaxtov+2) phi(imaxtov+1) = 0.5d0* . log(((1.0d0 - m(imaxtov)/(2.0d0*r_iso(imaxtov+1)))**2)/ . ((1.0d0 + m(imaxtov)/(2.0d0*r_iso(imaxtov+1)))**2)) phi(imaxtov+2) = 0.5d0* . log(((1.0d0 - m(imaxtov)/(2.0d0*r_iso(imaxtov+2)))**2)/ . ((1.0d0 + m(imaxtov)/(2.0d0*r_iso(imaxtov+2)))**2)) c fix up m at (imaxtov+1) and (imaxtov+2) m(imaxtov+1) = m(imaxtov) m(imaxtov+2) = m(imaxtov) c fix up press at (imaxtov+1) and (imaxtov+2) press(imaxtov+1) = 0.0d0 press(imaxtov+2) = 0.0d0 c fix up rho and eps at (imaxtov+1) and (imaxtov+2) rho(imaxtov+1) = 0.0d0 rho(imaxtov+2) = 0.0d0 eps(imaxtov+1) = 0.0d0 eps(imaxtov+2) = 0.0d0 c fix up r_schwar at (imaxtov+1) and (imaxtov+2) r_schwar(imaxtov+1) = r_iso(imaxtov+1) * . (1.0d0 + m(imaxtov)/(2.0d0*r_iso(imaxtov+1)))**2 r_schwar(imaxtov+2) = r_iso(imaxtov+2) * . (1.0d0 + m(imaxtov)/(2.0d0*r_iso(imaxtov+2)))**2 write(*,*) '/*/* tov_iso: phi at stellar surface is ', . phi(imaxtov) write(*,*) '/*/* tov_iso: final isotropic radius of star: ', . r_iso(imaxtov) write(*,*) '/*/* conformal factor at center of star: ', . sqrt(r_schwar(2)/r_iso(2)) c now, calculate the rest mass and total energy (minus gravity) c and proper_radius rest_mass = 0.5d0*rho(imaxtov)*(r_schwar(imaxtov)**3)/ . r_iso(imaxtov) if(mahc_init_eos_tables) then #ifdef CACTUSEOS_EOS_1D_BASE energy_mass = 0.5d0*((r_schwar(imaxtov)**3)/r_iso(imaxtov))* . rho(imaxtov)* . (1.0d0 + EOS_1D_epsfromrho(mahc_init_eos_handle,rho(imaxtov))) #endif else energy_mass = 0.5d0*((r_schwar(imaxtov)**3)/r_iso(imaxtov))* . (rho(imaxtov) + eosk*(rho(imaxtov)**gamma)/ . (gamma - 1.0d0)) endif proper_radius = 0.5d0*(1.0d0 + r_schwar(imaxtov)/r_iso(imaxtov)) do i=2,imaxtov-1 rest_mass = rest_mass + ((r_schwar(i)**3)/r_iso(i)) * . rho(i) if(mahc_init_eos_tables) then #ifdef CACTUSEOS_EOS_1D_BASE energy_mass = energy_mass + ((r_schwar(i)**3)/r_iso(i)) * . rho(i)* . (1.0d0 + EOS_1D_epsfromrho(mahc_init_eos_handle,rho(i))) #endif else energy_mass = energy_mass + ((r_schwar(i)**3)/r_iso(i)) * . (rho(i) + eosk * (rho(i)**gamma)/(gamma - 1.0d0)) endif proper_radius = proper_radius + r_schwar(i)/r_iso(i) enddo rest_mass = rest_mass * 4.0d0 * pi * r_iso(2) energy_mass = energy_mass * 4.0d0 * pi * r_iso(2) proper_radius = proper_radius * r_iso(2) write(*,*) '/*/* rest mass of star: ', . rest_mass write(*,*) '/*/* rest mass energy plus internal energy of star: ', . energy_mass write(*,*) '/*/* proper radius of star: ', . proper_radius write(*,*) '*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/' return end c**************************************************************** subroutine press_source_iso(press,m,r_schwar,r_iso,eosk,gamma,pi,retval, . mahc_init_eos_tables,mahc_init_eos_handle) implicit none CCTK_REAL press,m,r_schwar,r_iso,eosk,gamma,pi,retval integer mahc_init_eos_handle logical mahc_init_eos_tables CCTK_REAL tiny,rho,rhoeps c have to include EOS_1D_Base.inc due to the ifndef-define crap #ifdef CACTUSEOS_EOS_1D_BASE CCTK_REAL EOS_1D_pressfromrho CCTK_REAL EOS_1D_pressfromeps CCTK_REAL EOS_1D_epsfromrho CCTK_REAL EOS_1D_epsfrompress CCTK_REAL EOS_1D_rhofrompress CCTK_REAL EOS_1D_rhofromeps #endif if(press > 0.0d0) then tiny = 1.0d-30 if(mahc_init_eos_tables) then #ifdef CACTUSEOS_EOS_1D_BASE rho = EOS_1D_rhofrompress(mahc_init_eos_handle,press) rhoeps = rho*EOS_1D_epsfrompress(mahc_init_eos_handle,press) #endif else rho = (press/eosk)**(1.0d0/gamma) rhoeps = press / (gamma - 1.0d0) endif retval = - (rho + rhoeps + press)*(m + 4.0d0*pi*press*r_schwar**3)/ . ( r_iso * sqrt(r_schwar*(r_schwar - 2.0d0*m)) + tiny) else retval = -1.0d0 endif return end c**************************************************************** subroutine m_source_iso(press,m,r_schwar,r_iso,eosk,gamma,pi,retval, . mahc_init_eos_tables,mahc_init_eos_handle) implicit none CCTK_REAL press,m,r_schwar,r_iso,eosk,gamma,pi,retval integer mahc_init_eos_handle logical mahc_init_eos_tables CCTK_REAL tiny,rho,rhoeps c have to include EOS_1D_Base.inc due to the ifndef-define crap #ifdef CACTUSEOS_EOS_1D_BASE CCTK_REAL EOS_1D_pressfromrho CCTK_REAL EOS_1D_pressfromeps CCTK_REAL EOS_1D_epsfromrho CCTK_REAL EOS_1D_epsfrompress CCTK_REAL EOS_1D_rhofrompress CCTK_REAL EOS_1D_rhofromeps #endif if(press > 0.0d0) then tiny = 1.0d-30 if(mahc_init_eos_tables) then #ifdef CACTUSEOS_EOS_1D_BASE rho = EOS_1D_rhofrompress(mahc_init_eos_handle,press) rhoeps = rho*EOS_1D_epsfrompress(mahc_init_eos_handle,press) #endif else rho = (press/eosk)**(1.0d0/gamma) rhoeps = press / (gamma - 1.0d0) endif retval = 4.0d0 * pi * (rho + rhoeps) * r_schwar * r_schwar * . sqrt(r_schwar*(r_schwar - 2.0d0*m))/(r_iso + tiny) else retval = 0.0d0 endif return end c**************************************************************** subroutine mrest_source_iso(press,m,r_schwar,r_iso,eosk,gamma,pi,retval, . mahc_init_eos_tables,mahc_init_eos_handle) implicit none CCTK_REAL press,m,r_schwar,r_iso,eosk,gamma,pi,retval integer mahc_init_eos_handle logical mahc_init_eos_tables CCTK_REAL tiny,rho c have to include EOS_1D_Base.inc due to the ifndef-define crap #ifdef CACTUSEOS_EOS_1D_BASE CCTK_REAL EOS_1D_pressfromrho CCTK_REAL EOS_1D_pressfromeps CCTK_REAL EOS_1D_epsfromrho CCTK_REAL EOS_1D_epsfrompress CCTK_REAL EOS_1D_rhofrompress CCTK_REAL EOS_1D_rhofromeps #endif if(press > 0.0d0) then tiny = 1.0d-30 if(mahc_init_eos_tables) then #ifdef CACTUSEOS_EOS_1D_BASE rho = EOS_1D_rhofrompress(mahc_init_eos_handle,press) #endif else rho = (press/eosk)**(1.0d0/gamma) endif retval = 4.0d0 * pi * (rho) * r_schwar * r_schwar * . sqrt(r_schwar*(r_schwar - 2.0d0*m))/(r_iso + tiny) else retval = 0.0d0 endif return end c**************************************************************** subroutine phi_source_iso(press,m,r_schwar,r_iso,eosk,gamma,pi,retval) implicit none CCTK_REAL press,m,r_schwar,r_iso,eosk,gamma,pi,retval CCTK_REAL tiny if (press > 0.0d0) then tiny = 1.0d-30 retval = (m + 4.0d0 * pi * press * r_schwar**3)/ . (r_iso*sqrt(r_schwar*(r_schwar - 2.0d0*m)) + tiny) else retval = 0.0d0 endif return end c**************************************************************** subroutine rschwar_source_iso(m,r_schwar,r_iso,retval) implicit none CCTK_REAL m,r_schwar,r_iso,retval if(r_iso .le. 1.0d-16) then retval = 1.0d0 else retval = sqrt(r_schwar*(r_schwar - 2.0d0*m))/(r_iso ) endif return end c**************************************************************** c GETTOV_ISO c**************************************************************** c This routine returns the TOV information at a certain input c coordinate for isotropic tov information. (input coordinate c value is r_iso). c**************************************************************** subroutine gettov_iso(press,rho,eps,m,phi,r_schwar,r_iso, . presstov,rhotov,epstov,mtov,phitov,r_schwartov,r_isotov,ntov, . imaxtov,rmaxtov) implicit none integer ntov,imaxtov CCTK_REAL presstov(ntov),rhotov(ntov),epstov(ntov),mtov(ntov) CCTK_REAL phitov(ntov),r_schwartov(ntov),r_isotov(ntov) CCTK_REAL press,rho,eps,m,phi,r_schwar,r_iso CCTK_REAL rmaxtov c LOCAL VARS integer i,ilow,ihigh if(abs(r_iso) .ge. r_isotov(imaxtov+1)) then rho = 0.0d0 press = 0.0d0 eps = 0.0d0 m = mtov(imaxtov) r_schwar = r_iso * (1.0d0 + m/(2.0d0*r_iso))**2 phi = 0.5d0 * . log( ( (1.0d0 - m/(2.0d0*r_iso))**2) / . ( (1.0d0 + m/(2.0d0*r_iso))**2) ) else i = ( abs(r_iso) / r_isotov(2) + 1.5d0) ilow = i if(r_isotov(ilow) .le. r_iso) then ihigh = ilow+1 if(r_isotov(ihigh) .lt. r_iso) then write(*,*) 'gettov_iso: alg error 1 ' write(*,*) ' r_iso: ',r_iso write(*,*) ' ilow, r(ilow): ',ilow,r_isotov(ilow) write(*,*) ' ihigh, r(ihigh): ',ihigh,r_isotov(ihigh) STOP endif else ihigh = ilow ilow = ilow-1 if(r_iso .lt. r_isotov(ilow)) then write(*,*) 'gettov_iso: alg error 2 ' write(*,*) ' r_iso: ',r_iso write(*,*) ' ilow, r(ilow): ',ilow,r_isotov(ilow) write(*,*) ' ihigh, r(ihigh): ',ihigh,r_isotov(ihigh) write(*,*) ' imaxtov: ',imaxtov write(*,*) ' r_isotov(imaxtov+1): ',r_isotov(imaxtov+1) STOP endif endif if(i .eq. 1) then press = presstov(i) rho = rhotov(i) eps = epstov(i) m = mtov(i) phi = phitov(i) r_schwar = r_iso * (r_schwartov(2)/r_isotov(2)) else press = ( presstov(ihigh)*(r_iso - r_isotov(ilow)) + . presstov(ilow) *(r_isotov(ihigh) - r_iso))/ . (r_isotov(ihigh) - r_isotov(ilow)) rho = ( rhotov(ihigh)*(r_iso - r_isotov(ilow)) + . rhotov(ilow) *(r_isotov(ihigh) - r_iso))/ . (r_isotov(ihigh) - r_isotov(ilow)) eps = ( epstov(ihigh)*(r_iso - r_isotov(ilow)) + . epstov(ilow) *(r_isotov(ihigh) - r_iso))/ . (r_isotov(ihigh) - r_isotov(ilow)) m = ( mtov(ihigh)*(r_iso - r_isotov(ilow)) + . mtov(ilow) *(r_isotov(ihigh) - r_iso))/ . (r_isotov(ihigh) - r_isotov(ilow)) phi = ( phitov(ihigh)*(r_iso - r_isotov(ilow)) + . phitov(ilow) *(r_isotov(ihigh) - r_iso))/ . (r_isotov(ihigh) - r_isotov(ilow)) r_schwar = ( r_schwartov(ihigh)*(r_iso - r_isotov(ilow)) + . r_schwartov(ilow) *(r_isotov(ihigh) - r_iso))/ . (r_isotov(ihigh) - r_isotov(ilow)) endif endif return end c**************************************************************** c GETTOV_CART_ISO c**************************************************************** c This routine returns the TOV information at a certain input c isotropic coordinate. c**************************************************************** subroutine gettov_cart_iso(xx,yy,zz,press,rho,eps,gxx,gxy,gxz,gyy,gyz,gzz, . dxgxx,dxgxy,dxgxz,dxgyy,dxgyz,dxgzz, . dygxx,dygxy,dygxz,dygyy,dygyz,dygzz, . dzgxx,dzgxy,dzgxz,dzgyy,dzgyz,dzgzz, . alp,ax,ay,az, . eos_gammaloc,eos_kloc,centdenloc, . presstov,rhotov,epstov,mtov,phitov,r_schwartov, . r_isotov,ntov,imaxtov,rmaxtov,mahc_init_eos_handle) implicit none DECLARE_CCTK_PARAMETERS integer ntov,imaxtov,mahc_init_eos_handle CCTK_REAL presstov(ntov),rhotov(ntov),epstov(ntov),mtov(ntov) CCTK_REAL phitov(ntov),r_schwartov(ntov),r_isotov(ntov) CCTK_REAL xx,yy,zz,press,rho,eps,gxx,gxy,gxz,gyy,gyz,gzz CCTK_REAL dxgxx,dxgxy,dxgxz,dxgyy,dxgyz,dxgzz CCTK_REAL dygxx,dygxy,dygxz,dygyy,dygyz,dygzz CCTK_REAL dzgxx,dzgxy,dzgxz,dzgyy,dzgyz,dzgzz CCTK_REAL alp,ax,ay,az,eos_gammaloc,eos_kloc,centdenloc,rmaxtov c LOCAL VARS CCTK_REAL rcoord,mval,phival,grr,dgrr,rr,tiny,pi,dalpha integer i,ilow,ihigh,CCTK_Equals CCTK_REAL g,dg,lap,dlap,r_schwarval,dphidr_a,droverrdr_a CCTK_REAL dphidr_b,dphidr,droverrdr_b,droverrdr #ifdef CACTUSEOS_EOS_1D_BASE CCTK_REAL EOS_1D_epsfromrho #endif pi = 3.14159265358979d0 tiny = 1.0d-30 rcoord = sqrt(xx*xx + yy*yy + zz*zz) rr = rcoord if(abs(rcoord) .ge. r_isotov(imaxtov+1)) then c hey, at least we do not have to worry about r=0! g = (1.0d0 + mtov(imaxtov)/(2.0d0*rr))**4 dg = - 2.0d0 * mtov(imaxtov) * . ((1.0d0 + mtov(imaxtov)/(2.0d0*rr))**3)/(rr**3) lap = (1.0d0 - mtov(imaxtov)/(2.0d0*rr)) / . (1.0d0 + mtov(imaxtov)/(2.0d0*rr)) dlap = (mtov(imaxtov) / (rr + mtov(imaxtov)/2.0d0)**2)/rr c hydro vars press = 0.0d0 rho = 0.0d0 eps = 0.0d0 c metric gxx = g gxy = 0.0d0 gxz = 0.0d0 gyy = g gyz = 0.0d0 gzz = g c derivatives of metric divided by GDMFSOB factor of 2. dxgxx = xx * dg / 2.0d0 dygxx = yy * dg / 2.0d0 dzgxx = zz * dg / 2.0d0 dxgyy = xx * dg / 2.0d0 dygyy = yy * dg / 2.0d0 dzgyy = zz * dg / 2.0d0 dxgzz = xx * dg / 2.0d0 dygzz = yy * dg / 2.0d0 dzgzz = zz * dg / 2.0d0 dxgxy = 0.0d0 dygxy = 0.0d0 dzgxy = 0.0d0 dxgxz = 0.0d0 dygxz = 0.0d0 dzgxz = 0.0d0 dxgyz = 0.0d0 dygyz = 0.0d0 dzgyz = 0.0d0 alp = lap ax = xx * dlap / lap ay = yy * dlap / lap az = zz * dlap / lap return endif i = ( abs(rr) / r_isotov(2) + 1.5d0) ilow = i if(r_isotov(ilow) .le. rr) then ihigh = ilow+1 if(r_isotov(ihigh) .lt. rr) then write(*,*) 'gettov_iso: alg error 1 ' write(*,*) ' r_iso: ',rr write(*,*) ' ilow, r(ilow): ',ilow,r_isotov(ilow) write(*,*) ' ihigh, r(ihigh): ',ihigh,r_isotov(ihigh) STOP endif else ihigh = ilow ilow = ilow-1 if(rr .lt. r_isotov(ilow)) then write(*,*) 'gettov_iso: alg error 2 ' write(*,*) ' r_iso: ',rr write(*,*) ' ilow, r(ilow): ',ilow,r_isotov(ilow) write(*,*) ' ihigh, r(ihigh): ',ihigh,r_isotov(ihigh) write(*,*) ' imaxtov: ',imaxtov write(*,*) ' r_isotov(imaxtov+1): ',r_isotov(imaxtov+1) STOP endif endif c get interpolated values rr = rr + tiny if(ilow .eq. 1) then r_schwarval = rr * (r_schwartov(2)/r_isotov(2)) dphidr_a = 0.0d0 droverrdr_a = (r_schwartov(2)/r_isotov(2)) droverrdr_b = (r_schwartov(2)/r_isotov(2)) else r_schwarval = ( r_schwartov(ihigh)*(rr - r_isotov(ilow)) + . r_schwartov(ilow) *(r_isotov(ihigh) - rr))/ . (r_isotov(ihigh) - r_isotov(ilow)) dphidr_a = (phitov(ilow+1)-phitov(ilow-1))/ . (r_isotov(ilow+1) - r_isotov(ilow-1)) if(ilow .eq. 2) then droverrdr_a = (r_schwartov(2)/r_isotov(2)) else droverrdr_a = (r_schwartov(ilow+1)/r_isotov(ilow+1) - . r_schwartov(ilow-1)/r_isotov(ilow-1))/ . (r_isotov(ilow+1) - r_isotov(ilow-1)) endif droverrdr_b = (r_schwartov(ihigh+1)/r_isotov(ihigh+1) - . r_schwartov(ihigh-1)/r_isotov(ihigh-1))/ . (r_isotov(ihigh+1) - r_isotov(ihigh-1)) endif dphidr_b = (phitov(ihigh+1)-phitov(ihigh-1))/ . (r_isotov(ihigh+1) - r_isotov(ihigh-1)) press = ( presstov(ihigh)*(rr - r_isotov(ilow)) + . presstov(ilow) *(r_isotov(ihigh) - rr))/ . (r_isotov(ihigh) - r_isotov(ilow)) rho = ( rhotov(ihigh)*(rr - r_isotov(ilow)) + . rhotov(ilow) *(r_isotov(ihigh) - rr))/ . (r_isotov(ihigh) - r_isotov(ilow)) eps = ( epstov(ihigh)*(rr - r_isotov(ilow)) + . epstov(ilow) *(r_isotov(ihigh) - rr))/ . (r_isotov(ihigh) - r_isotov(ilow)) mval = ( mtov(ihigh)*(rr - r_isotov(ilow)) + . mtov(ilow) *(r_isotov(ihigh) - rr))/ . (r_isotov(ihigh) - r_isotov(ilow)) phival = ( phitov(ihigh)*(rr - r_isotov(ilow)) + . phitov(ilow) *(r_isotov(ihigh) - rr))/ . (r_isotov(ihigh) - r_isotov(ilow)) dphidr = (dphidr_b*(rr - r_isotov(ilow)) + . dphidr_a*(r_isotov(ihigh) - rr))/ . (r_isotov(ihigh) - r_isotov(ilow)) droverrdr = (droverrdr_b*(rr - r_isotov(ilow)) + . droverrdr_a*(r_isotov(ihigh) - rr))/ . (r_isotov(ihigh) - r_isotov(ilow)) g = (r_schwarval/rr)**2 dg = 2.0d0 * r_schwarval * droverrdr / (rr**2) lap = exp(phival) dlap = lap * dphidr / rr c hydro values: rho is already set, use EOS here to get eps c from press and rho if(.not.(CCTK_EQUALS(mahc_init_eos,'tables'))) then eps = press / ((eos_gammaloc - 1.0d0)*rho + tiny) else #ifdef CACTUSEOS_EOS_1D_BASE eps = EOS_1D_epsfromrho(mahc_init_eos_handle,rho + tiny) #else write(*,*) 'if you want to use tables in mahc_init_data, you' write(*,*) 'need to compile in with EOS_1D_Base' call CCTK_WARN(0,"tov_isotropic.F") #endif endif c metric gxx = g gxy = 0.0d0 gxz = 0.0d0 gyy = g gyz = 0.0d0 gzz = g c derivatives of metric divided by GDMFSOB factor of 2. dxgxx = xx * dg / 2.0d0 dygxx = yy * dg / 2.0d0 dzgxx = zz * dg / 2.0d0 dxgyy = xx * dg / 2.0d0 dygyy = yy * dg / 2.0d0 dzgyy = zz * dg / 2.0d0 dxgzz = xx * dg / 2.0d0 dygzz = yy * dg / 2.0d0 dzgzz = zz * dg / 2.0d0 dxgxy = 0.0d0 dygxy = 0.0d0 dzgxy = 0.0d0 dxgxz = 0.0d0 dygxz = 0.0d0 dzgxz = 0.0d0 dxgyz = 0.0d0 dygyz = 0.0d0 dzgyz = 0.0d0 alp = lap ax = xx * dlap / lap ay = yy * dlap / lap az = zz * dlap / lap return end