      program spacetravel ! Bernd Berg Apr 1 2010.
C Plot time, distance, received time as function of proper time tau.
      implicit real*8 (a-h,o-z)
      parameter(iuo=6,iud=9,ngnu=200, zero=0.d0,one=1.d0)
      parameter(g=9.81d0,year=365*24*3600.d0,c=299792458.d0,alpha=g/c)
      parameter(tau0p25=5.d0,tau0p50=2*tau0p25,tau0p75=3*tau0p25)
      parameter(tau1p00=4*tau0p25, t0=0.d0)
      write(iuo,'(/,"Rocket Parameters used:",/)')
      write(iuo,'(2x,"alpha =",g15.6,/)') alpha
      write(iuo,'(6x,"tau",7x,"t",9x,"x",7x,"tnews",/)')
      open(iud,file="spacetravel.d",form='formatted',status='unknown')
      do ignu=0,ngnu
        tau=ignu*tau1p00/ngnu
        call sspacetravel(tau,t,x,tnews)
        write(iud,'(4F10.4)') tau,t,x,tnews
        if(ignu.eq.ngnu/4) write(iuo,'(4F10.4)') tau,t,x,tnews
        if(ignu.eq.ngnu/2) write(iuo,'(4F10.4)') tau,t,x,tnews
        if(ignu.eq.3*ngnu/4) write(iuo,'(4F10.4)') tau,t,x,tnews
        if(ignu.eq.ngnu) write(iuo,'(4F10.4,/)') tau,t,x,tnews
      enddo
      close(iud)
      stop "spacetravel: All done."
      end


      subroutine sspacetravel(tau,t,x,tnews) ! Bernd Berg Mar 1 2009.
C Input:   Proper time tau in units of years.
C Output:  t time on earth.
C          x distance as seen at time t from earth.
C          tnews time of news received from earth via light signal.
      implicit real*8 (a-h,o-z)
      parameter(iuo=6,iud=9,ngnu=200, zero=0.d0,one=1.d0)
      parameter(g=9.81d0,year=365*24*3600.d0,c=299792458.d0,alpha=g/c)
      parameter(tau0p25=5.d0,tau0p50=2*tau0p25,tau0p75=3*tau0p25)
      parameter(tau1p00=4*tau0p25, t0=0.d0)
      if(tau.lt.zero) stop "sspacetravel: tau.lt.zero"
c
      if(tau.le.tau0p25) then
        zeta=alpha*tau*year
        t=sinh(zeta)/alpha/year
        x=(cosh(zeta)-one)/alpha/year
      else
        zeta0=alpha*tau0p25*year
        t01=sinh(zeta0)/alpha/year
        x01=(cosh(zeta0)-one)/alpha/year
      endif
c
      if(tau.gt.tau0p25) then
        if(tau.le.tau0p50) then
          zeta=zeta0-alpha*(tau-tau0p25)*year
          t=t01+(t01-sinh(zeta)/alpha/year)
          x=x01+(x01-(cosh(zeta)-one)/alpha/year)
        else
          t02=2*t01
          x02=2*x01
        endif
      endif
c
      if(tau.gt.tau0p50) then
        if(tau.le.tau0p75) then
          zeta=alpha*(tau0p50-tau)*year
          t=t02-sinh(zeta)/alpha/year
          x=x02-(cosh(zeta)-one)/alpha/year
        else
          zeta0=-alpha*tau0p25*year
          t03=3*t01
        endif
      endif
c
      if(tau.gt.tau0p75) then
        zeta=zeta0+alpha*(tau-tau0p75)*year
        t=t03+(t01+sinh(zeta)/alpha/year)
        x=x01-(x01-(cosh(zeta)-one)/alpha/year)
      endif
c
      if(tau.gt.tau1p00) stop "sspacetravel: tau.gt.tau1p00."
      t=t+t0
      tnews=t-x
      return
      end
