      program sun ! BB May 20, Nov 10 2017, Nov 20 2018.
c Solar system: Schwarzschild radius of sun, coordinate distance, ruler 
c distance and radar distance from earth to the surface of the sun.
      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,sr=695700.d03) ! Sun mass [kg], radius [m].
      parameter(r1=149.6d09) ! 1 [au] Distance earth-sun [m]. 
      ltest=.false.
      ltest=.true.
      write(iuo,'(/," Program sun.f:",)')
      write(iuo,'(/," Solar system: Schwarzschild radius, ruler",
     &              " and radar distance.")')
      write(iuo,'(/," c  [m/s]",10x,"=",g11.4,4x,"speed of light")') c
      write(iuo,'(" G  [m^3/(kg*s^2)] =",g11.4,
     &        4x,"gravitational constant")') G
      write(iuo,'(" y  [s]            =",g11.4,4x,"year")') y
      write(iuo,'(" ly [m]",12x,"=",g11.4,4x,"light year xly")') xly
      write(iuo,'(" sm [kg]",11x,"=",g11.4,4x,"mass  of the sun")') sm
      write(iuo,'(" sr [m]",12x,"=",g11.4,4x,"radius of the sun")') sr
      write(iuo,'(" r1 [m]",12x,"=",g11.4,4x,"distance earth-sun",
     &            " = 1 astronomical unit [au] ")') r1
      sg=G*sm/sr**2
      write(iuo,'(" sg [m/s^2]",8x,"=",f7.1,1x,
     &  "= G*sm/sr^2 sun surface gravitational acceleration",/)') sg
      rs=2*G*sm/c**2 ! Schwarzschild radius.
      write(iuo,'("(a) Sun Schwarzschild radius rs [m] =",1g12.5)') rs
      xmor=rs/sr
      write(iuo,'(4x,"rs/sr [dimensionless]",11x,"=",g12.4,/)') xmor
c
      write(iuo,'(" Units of coice c=G=1 and everything in seconds:")') 
      s=1.d00
      xm=s/c
      write(iuo,'(/," meter [s]",10x,"=",g11.4,4x,"= xm [s]")') xm
      xkg=G*xm**3/s**2
      write(iuo,'(" kg [s]   ",10x,"=",g11.4,4x,"= xkg [s]")') xkg
      sg=(sm*xkg)/(sr*xm)**2
      cg=sg*c ! To convert from [1/s] to [m/s^2].
      write(iuo,'(" sg [s]",13x,"=",g11.4,4x,"grav acceleration")') sg
      write(iuo,'(" cg [m/s^2]",9x,"=",g11.4,4x,"back conversion")') cg
      sms=(sm*xkg)
      write(iuo,'(" sm [s]",13x,"=",g11.4,4x,"sun mass sms")') sms
      srs=(sr*xm)
      write(iuo,'(" sr [s]",13x,"=",g11.4,4x,"sun radius srs")') srs
      xmor=2*sms/srs
      write(iuo,'(" 2*sms/srs [dimensionless]",1x,"=",g12.4,/)') xmor
c
      write(iuo,'("(b) Coordinate and ruler distance:",/,/,4x,
     &  "dr coordinate distance earth to radius of the sun,",/,4x,
     &  "dl ruler distance earth to radius of the sun.",/)') 
      dr=(r1-sr)*xm
      drul=sms*log(r1/sr)
      dl=dr+drul
      write(iuo,'(" dr [s]",13x,"=",g17.10," coordinate distance")') dr
      write(iuo,'(" dl [s]",13x,"=",g17.10,1x,"ruler distance")') dl
      write(iuo,'(" drul [s]",11x,"=",f13.7,5x,"= dl-dr =",g10.3)') 
     &              drul,drul
      write(iuo,'(" Conversion to meter:")') 
      smsmeter=sms/xm
      xlnrat=log(r1/sr)
      print*,"smsmeter, xlnrat =",smsmeter,xlnrat
      drm=dr/xm
      dlm=dl/xm
      drulm=drul/xm
      write(iuo,'(" dr [m]",13x,"=",g17.10," coordinate distance")') drm
      write(iuo,'(" dl [m]",13x,"=",g17.10,1x,"ruler distance")') dlm
      write(iuo,'(" drul [m]",11x,"=",f7.1,11x,"= dl-dr")') drulm
c
      write(iuo,'(/,"(c) Radar distance:",/)') 
      tdil=-sms/(r1*xm)
      write(iuo,'(" tdil [s]",11x,"=",g12.4,5x," dilation fraction")')
     &              tdil
      fdil=1.d00+tdil
      write(iuo,'(" fdil [s]",11x,"=",g18.10,"dilation factor")') fdil
      write(iuo,'(" 2*drul [s]",9x,"=",f13.7)') 2*drul
      R=fdil*(dr+2*drul)
      write(iuo,'(" R [s]",14x,"=",g17.10," radar distance")') R
      drad=R-dr
      write(iuo,'(" drad [s]",11x,"=",f13.7,6x,"= R-dr =",g10.3)') 
     &              drad,drad
      write(iuo,'(" fdil dilation factor accounts for 2*drul-drad",1x,
     &                          "=",g10.3)') 2*drul-drad
      dradl=R-dl
      write(iuo,'(" dradl [s]",10x,"=",f13.7,6x,"= R-dl =",g10.3)') 
     &              dradl,dradl
      write(iuo,'(/," radar distance R is larger",
     &  " than the other distance definitions.",/)')
      write(iuo,'(" Conversion to meter:")') 
      R=R/xm
      drad=drad/xm
      dradl=dradl/xm
      write(iuo,'(" R [m]",14x,"=",g17.10,1x,"radar distance")') R
      write(iuo,'(" drad [m]",11x,"=",f8.1,10x,"= R-dr")') drad
      write(iuo,'(" dradl [m]",10x,"=",f8.1,10x,"= R-dl")') dradl
c
      if(ltest) stop "sun: ltest."
      stop "sun: all done."
      end
