# hp41programs

Twin Einstein's Twin Paradox for the HP-41

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.

-Line 03:  this constant = 365.25*86400*9.80665/(299792458*4)

The "standard" gravity 9.80665 doesn't seem quite "standardized"
For instance, I've seen here and there the value 9.80616
In this case, replace line 02 by  0.2580609239

 01  LBL "ETP" 02  X<>Y 03  .2580738189 04  *             05  STO M   06  *             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

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.
-Line 14:  this constant is 4 times the value appearing line 02 of the "ETP" listing above.

 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 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 - paragraph 5 )
"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 )

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 )