The Museum of HP Calculators


Einstein's Twin Paradox for the HP-41

This program is Copyright © 2005 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

-The following programs use Einstein's special relativity to study uniformly ( and non-uniformly ) accelerated linear motion of a spacecraft.
-In special relativity, the speed is also  v = dx/dt , but acceleration is  a = ( dv/dt ) ( 1 - v2/c2 )-3/2   where  t = time passed on Earth
-The proper Time T ( passed on board ship ) verifies  dT = dt ( 1 - v2/c2 )1/2

-In all these programs, we use the following units:

    -Times ( t , T ) are expressed in Julian years ( 365.25 days )
    -Abcissas and distances ( x , D ) are expressed in light-years ( one light-year = 9.460,730,472,580,8 1015 m. )
    -The unit of accelerations ( a ) is the "standard" gravity on Earth  g  = 9.80665 m/s2

    -The velocity of light = c = 299,792,458 m/s  ( from the definition of 1 meter )

Remark:   Integrals are denoted  §
 

1°) A Trip to the Stars
 

-A spaceship travels to a star ( situated at a distance D from the Earth ) at a constant acceleration "a" , in order to simulate an artificial gravity.
-Halfway to the star, the constant acceleration "a" becomes a constant deceleration "-a"
-The return journey goes off similarly.

                                                    Return Journey
              <<<<<Deceleration "a" <<<<<  |  <<<<<Acceleration "-a" <<<<<<<
 Earth=O-------------------------------D/2------------------------------------D=(Star)--------------->   x
              >>>>>Acceleration "a" >>>>>  |  >>>>>Deceleration "-a" >>>>>>>
                                                  Outward Journey

-The graphic representation of  x(t) looks like this:

        x
         |
         | D                                              *
         |                                    *                         *
         |                             *                                       *
         | D/2                 *                                                  *
         |                  *                                                              *
         |           *                                                                            *
         |*---------------|----------------|-----------------|-----------------*-----  t
    Takeoff                          Spacecraft reaches the star                        Landing
 

  Formulae:   Proper Time on board ship = TS = (4c/a) Arccosh ( 1 + a.D/(2c2) )
                              Time passed on Earth = tE = (4c/a) ( aD/c2 + a2D2/(4c4) )1/2
        Maximum speed ( when x = D/2 ) = vmax/c = Tanh µ   where µ = a.TS/(4c)

 -Times concern the whole journey:  Earth-Star-Earth.

-Synthetic register M may be replaced by register R00.

01  LBL "ETP"
02  X<>Y
03  .2580738189                    this constant = 365.25*86400*9.80665/(299792458*4)
04  *                                       the "standard" gravity 9.80665 doesn't seem quite "standardized"
05  STO M                             for instance, I've seen here and there the value 9.80616
06  *                                       in this case, replace line 02 by  0.2580609239
07  ST+ X
08  ENTER^
09  ENTER^
10  X^2
11  LASTX
12  ST+ X
13  +
14  SQRT
15  STO Z
16  +
17  LN1+X
18  0
19  X<> M
20  ST/ Z
21  /
22  R^
23  1
24  +
25  1/X
26  ASIN
27  COS
28  2
29  ST/ L
30  X<> L
31  SIN
32  X^2
33  ST+ X
34  X<>Y
35  R^
36  R^
37  END

( 64 bytes / SIZE 000 )
 
 
      STACK        INPUTS      OUTPUTS
           T             /        1-vmax/c
           Z             /          vmax/c
           Y             a             tE
           X             D            TS

 where  TS = time passed on board ship
             tE = time passed on Earth.
            vmax = maximum speed  ( when x = D/2 )
 

Example1:    On 2007/12/25  immediate boarding for the first interstellar flight in the direction of Toliman ( alpha Centauri ) Distance D = 4.343 light-years.
                     Fox Mulder remains on Earth but Dana Scully's flying saucer takes off at 10h13m TT , acceleration = 0.7g

   0.7    ENTER^
 4.343  XEQ "ETP"   >>>>      TS    =   8.837390060 years
                                 RDN        tE    = 13.09998313  years
                                 RDN    vmax/c  =  0.9211383858              vmax  = 92.11383858% * 299792458 m/s = 276150341 m/s
                                 RDN  1-vmax/c = 0.07886161418

-When she steps from the ship, Dana's watch shows 2016/10/26 , 6h46m41s
  wheras according to Mulder, it's 4h40m08s o'clock and we are on 2021/01/30
 

Example2:   With a = 1g and D = 2,200,000 light-years   ( approximately the distance of the Andromeda galaxy )

       1         ENTER^
 2200000      R/S        >>>>  56.71150062    years
                                   RDN  4,400,003.874  years
                                   RDN       vmax/c = 1
                                   RDN   1-vmax/c = 3.877715920 E-13  whence  vmax/c  was actually  0.999,999,999,999,612,228,408

-Caculating   1-vmax/c is only useful if the maximum  v/c  is very close to 1
-If you don't want to compute this value, replace lines 28 to 36 by   STO T   RDN
-If you prefer (1-v2max/c2)1/2 , replace lines 26 to 34 by   ENTER^   ASIN   COS
-If you don't want to know   vmax/c  and  1-vmax/c , delete lines 22 to 36.
 

2°) Uniformly Accelerated Motion
 

-This program converts t ( time on Earth ) , T ( time on the spacecraft ) and x ( distance from the Earth ) for a given constant acceleration "a"
-It also calculates the speed  v/c  ( and 1-v/c ).

   Formulae:   v/c = Tanh µ     where  µ = a.T/c  ,   v/c = ( at/c ).( 1 + a2t2/c2 )-1/2
                         t  = ( c/a ) Sinh µ
                         x = ( c2/a ) ( Cosh µ - 1 ) = ( c2/a ) ( ( 1 + a2t2/c2 )1/2 - 1 )

-Due to these formulae, uniformly accelerated motion is often called "hyperbolic motion"

-Synthetic register M may be replaced by register R00.
 

01  LBL "ETP2"
02  LBL A
03  SF 00
04  GTO 01
05  LBL B
06  CF 01
07  GTO 00
08  LBL C
09  SF 01
10  LBL 00
11  CF 00
12  LBL 01
13  X<>Y
14  1.032295276              this constant is 4 times the value appearing line 02 of the "ETP" listing above.
15  *
16  STO M
17  *
18  FS? 00
19  GTO 00
20  FS? 01
21  GTO 01
22  X^2
23  STO Y
24  1
25  +
26  SQRT
27  1
28  +
29  /
30  LBL 00
31  ENTER^
32  ENTER^
33  X^2
34  LASTX
35  ST+ X
36  +
37  SQRT
38  FS? 00
39  STO Z
40  +
41  LN1+X
42  RCL M
43  ST/ Z
44  GTO 02
45  LBL 01
46  E^X-1
47  LASTX
48  CHS
49  E^X-1
50  +
51  2
52  /
53  STO Y
54  RCL M
55  /
56  RCL Y
57  2
58  +
59  R^
60  *
61  SQRT
62  RCL M
63  LBL 02
64  /
65  R^
66  1
67  +
68  1/X
69  ASIN
70  COS
71  2
72  ST/ L
73  X<> L
74  SIN
75  X^2
76  ST+ X
77  X<>Y
78  R^
79  R^
80  CLA
81  END

( 128 bytes / SIZE 000 )
 

 1°)  a , x  >>>>  T , t , v/c
 
 
      STACK        INPUTS        XEQ A      OUTPUTS
           T             /         1-v/c
           Z             /           v/c
           Y             a            t
           X             x            T

Example:      0.7  ENTER^
                      10   XEQ "ETP2" or  XEQ A or [A] in user mode    >>>>    T = 3.870348933
                                                                                                      RDN      t = 11.29945015
                                                                                                      RDN    v/c = 0.992583500
                                                                                                      RDN 1-v/c = 0.007416499860

 2°)  a , t  >>>>  T , x , v/c
 
 
      STACK        INPUTS        XEQ B      OUTPUTS
           T             /         1-v/c
           Z             /           v/c
           Y             a            x
           X             t            T

Example:      0.7  ENTER^
                      10   XEQ B  or  [B] in user mode  >>>>    T = 3.702700069
                                                                           RDN     x  = 8.711423203
                                                                           RDN    v/c = 0.990559778
                                                                           RDN 1-v/c = 0.009440221722

 3°)  a , T  >>>>  t , x , v/c
 
 
      STACK        INPUTS        XEQ C      OUTPUTS
           T             /         1-v/c
           Z             /           v/c
           Y             a            x
           X             T            t

Example:      0.7  ENTER^
                      10   XEQ C  or  [C] in user mode  >>>>     t = 951.2809258
                                                                            RDN     x = 949.8980538
                                                                            RDN    v/c =  0.999998942
                                                                            RDN 1-v/c = 0.000001058151320

-If you don't want to compute 1-v/c, replace lines 71 to 79 by   STO T   RDN
-If you prefer (1-v2max/c2)1/2 , replace lines 69 to 77 by   ENTER^   ASIN   COS
-If you don't want to know   vmax/c  and  1-vmax/c , delete lines 65 to 79.
 

3°) An Example of non-uniformly Accelerated Motion
 

-The spaceship acceleration  a(T)  is now supposed to be known as a function of the proper time T ( on board )
-We have:

   a(T) = ( dv/dt ) ( 1 - v2/c2 )-3/2  =  ( dv/dT ) ( 1 - v2/c2 )-1   since   dT = dt ( 1 - v2/c2 )1/2     therefore,   dv/( 1 - v2/c2 ) = a(T).dT

 whence   v/c = Tanh ( (1/c) F(T) )        where   F(T) is the primitive of  a(T)  such that F(0) = 0    ( we assume  x = 0 & v = 0 when t = T = 0 )

-Likewise,  dt = dT ( 1 - v2/c2 )-1/2  =  dT ( 1 - Tanh2(F(T)) )-1/2 = dT Cosh(F(T))         whence  t = §0T  Cosh((1/c)F(u)).du          ( § = Integral )
  and  dx = v.dt = c.Tanh(F(T)).Cosh(F(T)).dT = c. Sinh(F(T)).dT                                  whence  x = c §0T  Sinh((1/c)F(u)).du

-"ETP3" evaluates  t ( time on Earth ), x ( abscissa ) and v/c ( speed of the spacecraft ) at a given proper time T,
  assuming the primitive F(T) may be expressed by elementary functions
 

Data Registers:              R00:  F name                                                  ( Register R00 is to be initialized before executing "ETP3" )
                               R01 = 0        R04 thru R08 are used by "GL3"
                               R02 = T        R12 , R13: temp
                               R03 = n = number of subintervals used in Gaussian integration.
           when the program stops:   R09 = t , R10 = x , R11 = v/c

Flags:  /
Subroutines:   "GL3"  ( cf "Numerical Integration for the HP-41" )
                   and a program which calculates F(T) assuming T is in X-register upon entry.
 

01  LBL "ETP3"
02  STO 02
03  X<>Y
04  STO 03
05  X<>Y
06  XEQ IND 00
07  1.032295276
08  STO 12
09  *
10  ST+ X
11  E^X-1
12  STO Y
13  2
14  +
15  /
16  STO 11
17  CLX
18  STO 01
19  RCL 03
20  LBL 00
21  STO 03
22  RCL 00
23  STO 13
24  "CT"
25  ASTO 00
26  XEQ "GL3"
27  STO 09
28  "ST"
29  ASTO 00
30  XEQ "GL3"
31  STO 10
32  RCL 13
33  STO 00
34  RCL 11
35  RCL 10
36  RCL 09
37  CLA
38  RTN
39  GTO 00
40  LBL "ST"
41  XEQ IND 13
42  RCL 12
43  *
44  E^X-1
45  LASTX
46  CHS
47  E^X-1
48  -
49  2
50  /
51  RTN
52  LBL "CT"
53  XEQ IND 13
54  RCL 12
55  *
56  E^X
57  LASTX
58  CHS
59  E^X
60  +
61  2
62  /
63  RTN
64  END

( 113 bytes / SIZE 014 )
 
 
      STACK        INPUTS      OUTPUTS
           Z             /           v/c
           Y             n            x
           X             T            t

 where  n is the number of subintervals over which the 3-point Gauss-Legendre formula is applied

Example:     acceleration = a(T) = 0.7 cos ( 360°(T/40) )    the whole journey will take 40 years on board.

-Here, we can easily find the required primitive:    F(T) = (14/pi) sin ( 360°(T/40) ) = (14/pi) sin ( 9°T )   So we load the following subroutine

 01  LBL "FF"   ( global label , at most 6 characters )
 02  9
 03  *
 04  SIN
 05  14
 06  *
 07  PI
 08  /
 09  RTN

- "FF" ASTO 00
-We set the HP-41 in DEG mode and if we choose  T = 1 ( year )  and  n = 2

  2  ENTER^
  1  XEQ "ETP3"  >>>>     t  =  1.088871751 year = R09    ( the difference between Earth time and ship time is already about 32 days )
                             RDN     x  = 0.376425754 light-year = R10
                             RDN    v/c = 0.616685490 = R11

-With T = 10 and n = 4

   4  ENTER^
  10 XEQ "ETP3"  >>>>   t = 190.9695910  ( in 1 minute )
                              RDN   x = 189.5025503
                              RDN  v/c = 0.999798046

-To have another approximation for the same T but a different n, simply key in the new n-value and press R/S , for instance:

    8  R/S  >>>>   t = 190.9695910
                RDN   x = 189.5025369
                RDN  v/c = 0.999798046

-Note that v/c is obtained exactly since no approximate integration formula is used.
-In the above example, the whole journey Earth-Star-Earth required 40 years on board ship but almost 764 years have passed on Earth!
-The distance Earth-Star was approximately 379 light-years.

-Results may be wrong if F(T) or its derivatives have singularities.
 

4°) Non-uniformly Accelerated Motion, a more general case
 

-We study the same problem as in paragraph 3 , but we suppose the primitive difficult ( or impossible ) to express analytically. We have to evaluate:

      t = §0T  Cosh( (1/c) §0u a(v).dv ).du
      x = c §0T  Sinh( (1/c) §0u a(v).dv ).du
     v/c = Tanh ( (1/c) §0T a(u).du )

Data Registers:              R00:  "a"  name                                                  ( Register R00 is to be initialized before executing "ETP4" )
                               R01 = 0        R04 thru R17 are used by "DIN2"
                               R02 = T        R18 , R22: temp
                               R03 = n = number of subintervals used in Gaussian integration.
           when the program stops:   R19 = t , R20 = x , R21 = v/c

Flags:  F04
Subroutines:   "DIN2" ( see the bottom of this page )
                          "GL3"  ( cf "Numerical Integration for the HP-41" )
                   and a program which calculates a(T) assuming T is in X-register upon entry.
 

01  LBL "ETP4"
02  1.032295276
03  STO 18
04  CLX
05  STO 01
06  RDN
07  STO 02
08  X<>Y
09  LBL 00
10  STO 03
11  RCL 00
12  STO 22
13  "T"
14  ASTO 00
15  "U"
16  ASTO 04
17  "V"
18  ASTO 05
19  "CH"
20  ASTO 17
21  SF 04
22  XEQ "DIN2"
23  STO 19
24  "SH"
25  ASTO 17
26  XEQ "DIN2"
27  STO 20
28  RCL 22
29  STO 00
30  XEQ "GL3"
31  RCL 18
32  *
33  ST+ X
34  E^X-1
35  STO Y
36  2
37  +
38  /
39  STO 21
40  RCL 20
41  RCL 19
42  CLA
43  RTN
44  GTO 00
45  LBL "T"
46  X<>Y
47  XEQ IND 22
48  RCL 18
49  *
50  RTN
51  LBL "U"
52  CLX
53  RTN
54  LBL "V"
55  RTN
56  LBL "SH"
57  E^X-1
58  LASTX
59  CHS
60  E^X-1
61  -
62  2
63  /
64  RTN
65  LBL "CH"
66  E^X
67  LASTX
68  CHS
69  E^X
70  +
71  2
72  /
73  RTN
74  END

( 156 bytes / SIZE 023 )
 
 
      STACK        INPUTS      OUTPUTS
           Z             /           v/c
           Y             n            x
           X             T            t

 where  n is the number of subintervals over which the 3-point Gauss-Legendre formula is applied

Example:     Let's take again  a(T) = 0.7 cos ( 360°(T/40) )   as if we couldn't calculate the primitive ( it will be a test case for our program )

-We load the following subroutine

 01  LBL "AA"   ( global label , at most 6 characters )
 02  9
 03  *
 04  COS                   Store 9 and 0.7 into registers  R23 and R24 and replace lines 02 & 05 by RCL 23 & RCL 24 respectively
 05  .7                        It will speed up execution.
 06  *
 07  RTN

- "AA" ASTO 00
-We set the HP-41 in DEG mode and with  T = 1 and  n = 2

  2  ENTER^
  1  XEQ "ETP4"  >>>>     t  =  1.088871751
                             RDN     x  = 0.376425755
                             RDN    v/c = 0.616685490

-With T = 10 and n = 4

   4  ENTER^
  10 XEQ "ETP4"  >>>>   t = 190.9695915      ( in 8mn17s )
                              RDN   x = 189.5025508
                              RDN  v/c = 0.999798046

-To have another approximation for the same T but a different n, simply key in the new n-value and press R/S , for instance:

    8  R/S  >>>>   t = 190.9695911
                RDN   x = 189.5025369
                RDN  v/c = 0.999798046

-Unlike "ETP3" , "ETP4" re-calculates v/c for each new n-value.
-When n is multiplied by 2 , execution time is approximately multiplied by 4
 

5°) Double Integrals   §ab G ( §u(x)v(x)  f(x,y) dy ) dx

-"ETP4" requires to evaluate integrals of this form, where G is a function of 1 variable.
-Modify the program listed in "Double Integrals for the HP-41" as follows:

-Replace line 01 by  LBL "DIN2"
-Delete line 96
-Add  3.6   /   FS? 04   GTO IND 17   after line  89

-In other words, the new listing looks like this:
 

  01  LBL "DIN2"
   .......................
  the same lines as "DIN"
   .......................

  87  RCL 14
  88  RCL 11
  89  *
  90  3.6
  91  /
  92  FS? 04
  93  GTO IND 17
  94  RTN
  95  LBL 04
  96  RCL 06
  97  RCL 09
  98  *
  99  3.6
100  /
101  STO 06
102  END

( 148 bytes / SIZE 018 )
 

-Set flag F04, store G name in register R17, and follow the same instructions as "DIN"

Example:    Compute  §12  Ln( 1 + §xx^2  dy/(x+y) ) dx

-Load these routines in program memory:

  01  LBL "T"
  02  +
  03  1/X
  04  RTN
  05  LBL "U"
  06  RTN
  07  LBL "V"
  08  X^2
  09  RTN
  10  LBL "W"
  11  LN1+X
  12  RTN

 "T"  ASTO 00   "U"  ASTO 04  "V"  ASTO 05   "W"  ASTO 17   SF 04
  1  STO 01  2  STO 02

-With   n = 2   2  STO 03  XEQ "DIN2"  >>>>  0.191218819
-With   n = 4   4  STO 03        R/S           >>>>  0.191218775

-If flag F04 is clear, "DIN2" performs the same double integral as "DIN"
 ( G is not taken into account )

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall