      program sden ! BB May 22 2017.
c Schwarzschild density as function of Schwarzschild radius.
      implicit real*8 (a-h,o-z)
      logical ltest
      parameter(c=2.998d+08) ! 2.998d+08 [m/s] = 1 -> c.
      parameter(G=6.674d-11) ! 6.674d-11 [m^3 kg^{-1} s^{-2}] = 1 -> G.
      parameter(y=365.d00*24.d00*3600.d0) ! Year in [s].
      parameter(iuo=6,xly=y*c)            ! Light year in [m].
      parameter(sm=1.9886d30,em=5.972d24) ! Sun mass [kg], Earth mass [kg].
      parameter(um=3.2d54) ! Mass of universe including dark matter.
      ltest=.true.
      ltest=.false.
      if(ltest) stop "sden: ltest."
c
      pi=acos(-1.d0)
      write(iuo,'(/," g [m/s^2], y [s], ly [m] =",2g12.4)') y,xly
      srs=2*G*sm/c**2 ! Schwarzschild radius sun.
      write(iuo,'(/," Schwarzschild radius sun [m]:",g19.4)') srs
      vols=4*pi*srs**3/3.d00
      write(iuo,'(" Schwarzschild volume sun [m^3]:",g18.4)') vols
      dens=sm/vols
      write(iuo,'(" Schwarzschild density sun [kg/m^3]:",g14.4)') dens
c
      sre=2*G*em/c**2 ! Schwarzschild radius earth.
      write(iuo,'(/," Schwarzschild radius earth [m]:",g18.4)') sre
      vole=4*pi*sre**3/3.d00
      write(iuo,'(" Schwarzschild volume earth [m^3]:",g16.4)') vole
      dene=em/vole
      write(iuo,'(" Schwarzschild density earth [kg/m^3]:",g12.4)')dene
c
      sru=2*G*um/c**2 ! Schwarzschild radius universe.
      write(iuo,'(/," Schwarzschild radius universe [m]:",g15.4)') sru
      write(iuo,'(30x,"[ly]:",g15.4)') sru/xly
c
      write(iuo,'(/," Density as function of Schwarschild radius:")')
      pm=1.67d-27
      write(iuo,'(/," Proton weight in [kg]:       ",1g12.4)') pm
      r=sre
      den=3*c**2/(8*pi*G*r**2)
      write(iuo,'(/," E Density [kg/m^3], Radius [m]:",2g12.4)') den,r
      r=srs
      den=3*c**2/(8*pi*G*r**2)
      write(iuo,'(" S Density [kg/m^3], Radius [m]:",2g12.4,/)') den,r
c
      rly=15.d09 ! Radius of the known universe.
      r=rly*xly ! Radius of the known universe.
      den=3*c**2/(8*pi*G*r**2)
      write(iuo,'(" U Density [kg/m^3], Radius [m]:",3G12.4," [ly]")') 
     &   den,r,rly
      xn=den/pm
      write(iuo,'(" So",f6.3," < 5 protons weight per [m^3].")') xn
      volu=4*pi*r**3/3.d00
      denu=um/volu
      write(iuo,'(" Density universe [kg/m^3]:",g12.4,/)') denu
c
      rly=46.d09 ! Radius of the known universe.
      r=rly*xly ! Radius of the known universe.
      den=3*c**2/(8*pi*G*r**2)
      write(iuo,'(" U Density [kg/m^3], Radius [m]:",3G12.4," [ly]")') 
     &   den,r,rly
      xn=den/pm
      write(iuo,'(" So",f6.3," < 1 protons weight per [m^3].")') xn
      volu=4*pi*r**3/3.d00
      denu=um/volu
      write(iuo,'(" Density universe [kg/m^3]:",g12.4)') denu
      xn=denu/pm
      write(iuo,'(" So",f6.3," < 6 protons weight per [m^3].",/)') xn
c
      stop "sden: all done."
      end
