The Museum of HP Calculators


Special Relativity for the HP-41

This program is Copyright © 2007 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

 1°)  The Lorentz Transformation
 2°)  Composition of Speeds & Relative Velocity
 3°)  Doppler Effect & Aberration

-See also "Einstein's Twin Paradox for the HP-41"
 

1°)  The Lorentz Transformation

                           y'
        y                  |       z'
        |                   |      /
        |       z          |    /
        |      /            |  /
        |    /              |/-------------------------- x'                ( Ox ) // ( O'x' )
        |  /               O'                                                          ( Oy ) // ( O'y' )
        |/----------------------------------- x                        ( Oz ) // ( O'z' )
      O

-We assume that O' has a constant velocity ß = ( ßx , ßy , ßz ) with respect to ( Ox , Oy , Oz ).
 and that the event ( 0 , 0 , 0 , 0 ) has the same coordinates in both reference-frames.

-The following program performs the conversion  ( x , y , z , ct )  < >  ( x' , y' , z' , ct' )  for a given event
 

Formulae:     With  r(x,y,z)  and  r'(x',y',z')  ,  g = ( 1 - b2 ) -1/2  ,   ß = || ß || = ( ßx2 +  ßy2z2 ) 1/2

   r  =  r' + ß.g [ g/(g+1) ß.r' + c t' ]
  c t = g ( c t' + ß.r' )
                                                                           where  ß.r and ß.r' are the dot products
   r'  =  r + ß.g [ g/(g+1) ß.r - c t ]
  c t' = g ( c t' - ß.r )
 

Data Registers:              R00 = Dilatation Factor = g = ( 1 - ß2 ) -1/2          ( Registers R01 thru R03 are to be initialized before executing "LRTF" )

                                        R01 = ßx         R04 = x    ( or x' if SF 00 )             R08 = the invariant    c2 t2 - x2 - y2 - z2 = c2 t'2 - x'2 - y'2 - z'2
                                        R02 = ßy         R05 = y    ( or y' if SF 00 )
                                        R03 = ßz          R06 = z    ( or z' if SF 00 )
                                                                 R07 = ct  ( or ct' if SF 00 )

Flags:   F00        CF 00 for the transformation  ( r , ct )   >>>  ( r' , ct' )
                            SF 00 for the transformation  ( r' , ct' )  >>>  ( r , ct )
Subroutines: /
 

01  LBL "LRTF"
02  STO 04
03  X^2
04  X<>Y
05  STO 05
06  X^2
07  +
08  X<>Y
09  STO 06
10  X^2
11  +
12  X<>Y
13  STO 07
14  X^2
15  X<>Y
16  -
17  STO 08
18  1
19  RCL 01
20  X^2
21  RCL 02
22  X^2
23  RCL 03
24  X^2
25  +
26  +
27  -
28  SQRT
29  1/X
30  STO 00
31  RCL 07
32  RCL 01
33  RCL 04
34  *
35  RCL 02
36  RCL 05
37  *
38  +
39  RCL 03
40  RCL 06
41  *
42  +
43  FS? 00
44  ST+ Y
45  FC? 00
46  ST- Y
47  RCL 00
48  ST* Z
49  ST* Y
50  1
51  +
52  /
53  RCL 07
54  FC? 00
55  CHS
56  +
57  RCL 00
58  *
59  RCL 03
60  X<>Y
61  *
62  RCL 02
63  LASTX
64  *
65  RCL 01
66  ST* L
67  CLX
68  RCL 06
69  ST+ Z
70  CLX
71  RCL 05
72  ST+ Y
73  CLX
74  RCL 04
75  ST+ L
76  X<> L
77  END

( 98 bytes / SIZE 009 )
 
 
 STACK(CF00)        INPUTS      OUTPUTS
           T            c.t           c.t'
           Z             z            z'
           Y             y            y'
           X             x            x'
 STACK(SF 00)        INPUTS      OUTPUTS
           T            c.t'           c.t
           Z             z'            z
           Y             y'            y
           X             x'            x

Example1:    Velocity  ß = ( 0.4 ; 0.5 ; 0.6 )   Event = ( x , y , z , c.t ) = ( 1 , 2 , 3 , 4 )

  0.4  STO 01
  0.5  STO 02
  0.6  STO 06

   CF 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "LRTF"  >>>>   x =  -0.5324
                               RDN   y  =   0.0846
                               RDN   z  =   0.7015
                               RDN  c.t =   1.6681

-We also have the dilatation factor g = R00 = 2.0851  and the invariant   c2 t2 - x2 - y2 - z2 = R08 = 2
 

Example2:    With the same  ß = ( 0.4 ; 0.5 ; 0.6 )   Event = ( x' , y' , z' , c.t' ) = ( 1 , 2 , 3 , 4 )

    SF 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "LRTF"  >>>>   x =  6.1401
                               RDN   y  =  8.4251
                               RDN   z  = 10.7102
                               RDN  c.t = 15.0130
 

Notes:

-If R08 > 0  it is a time-like event.
-If R08 < 0  it is a space-like event.
 

2°)  Composition of Speeds & Relative Velocity
 

-We assume that the units are chosen so that  c = speed of light = 1

                           y'
        y                  |       z'
        |                   |      /
        |       z          |    /          . P = Particle
        |      /            |  /
        |    /              |/-------------------------- x'
        |  /               O'
        |/----------------------------------- x
      O

- O' has a constant velocity u = ( ux , uy , uz ) with respect to ( Ox , Oy , Oz )

Problem1:   The speed of a particle P is v = ( vx , vy , vz ) with respect to ( Ox' , Oy' , Oz' ). Calculate its speed w with respect to ( Ox , Oy , Oz ) ?

Answer:   w = ( wx , wy , wz )  with

   wx = [  vx + g.ux ( g u.v /(g+1) + 1 ] / [ g ( 1 + u.v ) ]
   wy = [ vy + g.uy ( g u.v /(g+1) + 1 ] / [ g ( 1 + u.v ) ]             u.v = dot product     g = ( 1 - ux2 - uy2 - uz2 ) -1/2
   wz = [  vz + g.uz ( g u.v /(g+1) + 1 ] / [ g ( 1 + u.v ) ]

Problem2:   The speed of this particle P is v = ( vx , vy , vz ) with respect to ( Ox , Oy , Oz ). Calculate its relative velocity w with respect to ( Ox' , Oy' , Oz' ) ?

Answer:   w = ( wx , wy , wz )  with

   wx = [  vx + g.ux ( g u.v /(g+1) - 1 ] / [ g ( 1 - u.v ) ]
   wy = [ vy + g.uy ( g u.v /(g+1) - 1 ] / [ g ( 1 - u.v ) ]
   wz = [  vz + g.uz ( g u.v /(g+1) - 1 ] / [ g ( 1 - u.v ) ]

-Since these formulae only differ in a few changes of sign, we can use a unique program and a flag.
 

Data Registers:              R00 = Dilatation Factor = g = ( 1 - u2 ) -1/2        ( Registers R01 thru R06 are to be initialized before executing "SPEED" )

                                        R01 = ux        R04 = vx         R07 = wx
                                        R02 = uy       R05 = vy         R08 = wy
                                        R03 = uz        R06 = vz          R09 = wz

Flag:   F01        CF 01  for the composition of speeds
                          SF 01  for the relative velocity
Subroutines: /
 
 

01  LBL "SPEED"
02  1
03  RCL 01
04  X^2
05  RCL 02
06  X^2
07  RCL 03
08  X^2
09  +
10  +
11  -
12  SQRT
13  1/X
14  STO 00
15  RCL 01
16  RCL 04
17  *
18  RCL 02
19  RCL 05
20  *
21  +
22  RCL 03
23  RCL 06
24  *
25  +
26  STO 07
27  *
28  X<>Y
29  ST* 07
30  FC? 01
31  ST+ 07
32  FS? 01
33  ST- 07
34  1
35  +
36  /
37  *
38  FS? 01
39  CHS
40  +
41  STO 08
42  RCL 03
43  *
44  RCL 06
45  FS? 01
46  CHS
47  +
48  RCL 02
49  RCL 08
50  *
51  RCL 05
52  FS? 01
53  CHS
54  +
55  RCL 01
56  RCL 08
57  *
58  RCL 04
59  FS? 01
60  CHS
61  +
62  RCL 07
63  ST/ T
64  ST/ Z
65  /
66  STO 07
67  X^2
68  X<>Y
69  STO 08
70  X^2
71  +
72  X<>Y
73  STO 09
74  X^2
75  +
76  SQRT
77  RCL 09
78  RCL 08
79  RCL 07
80  END

( 101 bytes / SIZE 010 )
 
 
      STACK        INPUTS      OUTPUTS
           T             /            w
           Z             /            wz
           Y             /            wy
           X             /            wx

  where  w is the norm of the 3-vector w = ( wx , wy , wz )

-if CF 01  w = the composition of  u  and v
-if SF 01  w = the relative speed  ( with respect to the 2nd reference-frame )

Example:    u = ( 0.4 , 0.5 , 0.6 )   v =  ( 0.27 , 0.37 , 0.47 )

    0.4  STO 01       0.27  STO 04
    0.5  STO 02       0.37  STO 05
    0.6  STO 03       0.47  STO 06

 CF 01  XEQ "SPEED"  >>>>   wx =  0.434880406           ( execution time = 3 seconds )
                                      RDN    wy =  0.553496668
                                      RDN    wz =  0.672112930
                                      RDN    w  =  0.973249875

 SF 01  XEQ "SPEED"  >>>>    wx =  -0.270737319
                                      RDN    wy =  -0.301747643
                                      RDN    wz =  -0.332757967
                                      RDN    w  =   0.524478981

Notes:

-The composition of speeds is not commutative in the general case.
-However, the norms  w  are identical
-There will be a DATA ERROR line 12 or 13 if  u = || u ||  is equal to or greater than 1
-However, the program also works if  v = || v || = 1 ( a photon )
                                               and if  v = || v ||  is greater than 1. In this case, the particle P is a ( hypothetic ) tachyon!

-The formulas are much simpler  if  uy = uz = vy =  vz = 0
 then  wy = wz = 0  and for the composition of speeds:
 wx = ( ux + vx ) / ( 1 + ux vx )

-If you only want to compute the magnitude  w , the routine hereafter is shorter than "SPEED"
-It uses the formulae:

  w2 = [ ( u + v )2 - ( u x v )2 ] / ( 1 + u.v )2    for the composition of speeds           where   u.v  is the dot product
  w2 = [ ( u - v )2 - ( u x v )2 ] / ( 1 - u.v )2     for the relative velocities                     and   u x v  is the cross-product
 

Data Registers:              R00 = w                                   ( Registers R01 thru R06 are to be initialized before executing "MAGW" )

                                        R01 = ux        R04 = vx
                                        R02 = uy       R05 = vy
                                        R03 = uz        R06 = vz

Flag:   F01        CF 01  for the composition of speeds
                          SF 01  for the relative speed

Subroutines:  "DOT" & "CROSS" ( cf "Dot-product & Cross-product for the HP-41" )
 
 

01  LBL "MAGW"
02  1
03  RCL 06
04  RCL 05
05  RCL 04
06  XEQ "CROSS"
07  X^2
08  X<>Y
09  X^2
10  +
11  X<>Y
12  X^2
13  +
14  RCL 01
15  RCL 04
16  FS? 01
17  CHS
18  +
19  X^2
20  RCL 02
21  RCL 05
22  FS? 01
23  CHS
24  +
25  X^2
26  +
27  RCL 03
28  RCL 06
29  FS? 01
30  CHS
31  +
32  X^2
33  +
34  X<>Y
35  -
36  SQRT
37  STO 00
38  SIGN
39  4.006
40  XEQ "DOT"
41  FS? 01
42  CHS
43  1
44  +
45  ST/ 00
46  RCL 00
47  END

( 75 bytes / SIZE 007 )
 
 
      STACK        INPUTS      OUTPUTS
           X             /            w

  where  w is the norm of the 3-vector w = ( wx , wy , wz )

-if CF 01  w = the composition of  u and v
-if SF 01  w = the relative speed  ( with respect to the 2nd reference-frame )

Example:    u = ( 0.4 , 0.5 , 0.6 )   v =  ( 0.27 , 0.37 , 0.47 )

    0.4  STO 01       0.27  STO 04
    0.5  STO 02       0.37  STO 05
    0.6  STO 03       0.47  STO 06

 CF 01  XEQ "MAGW"  >>>>   w  =  0.973249875

 SF 01  XEQ "MAGW"  >>>>    w  =  0.524478980
 
 

3°)  Doppler Effect & Aberration
 

-The frequency and the direction of a photon are changed by the relative motion of the source and the observer.
-The following program performs these conversions.

    S = Source                                                        S = Source

     S\\---> ß                                                      S/ /---> ß
        \ \                                           y                 /
         \                                         |               /   /
     cS \   \ c0                                   |         cS/    / c                                                 ß  = velocity of the source , we assume ß // (Ox)
          .     .                                     |        .      .                                                      cS = velocity of a photon,  seen in the source reference-frame
          . µS   .  µ0                             |     . µS   . µ0                                                   c0  = velocity of the photon,  seen in the observer reference-frame
   ----------------------------------|---------------------------------- x
                                                     O = Observer                                                  µS = the angle ( ß , -cS ) = ( Ox , -cS )
                                                                                                                             µ0  = the angle ( ß , -c0 ) = ( Ox , -c0 ) is measured by the observer
 

                     fS = frequency in the source reference-frame                                     Of course  || cS || = || c0 || = 1
                     f0  = frequency measured by the observer
 

Formulae:

  f0  =  fS ( 1 - ß2 ) -1/2  ( 1 - ß cos µS )                cos µ0 = ( cos µS - ß ) ( 1 - ß cos µS ) -1
  fS =   f0 ( 1 - ß2 ) -1/2  ( 1 + ß cos µ0 )               cos µS = ( cos µ0 + ß ) ( 1 + ß cos µ0 ) -1

-We assume that µS and µ0  are between 0 and 180°  and  ß  is non-negative
-The source is moving away if µ0 < 90°
 

Data Registers:  /
Flag:     F00         CF 00 for the conversion  ( f0 , µ0 )  >>>  ( fS , µS )
                              SF 00 ------------------  ( fS , µS )  >>>  ( f0 , µ0 )
Subroutines: /
 

01  LBL "DAB"
02  RCL Z
03  FS? 00
04  CHS
05  STO T
06  X^2
07  SIGN
08  ST- L
09  X<> L
10  CHS
11  SQRT
12  /
13  RDN
14  COS
15  ST+ Z
16  *
17  1
18  +
19  ST* Z
20  /
21  ACOS
22  X<>Y
23  END

( 38 bytes / SIZE 000 )
 
 
 STACK(CF 00)        INPUTS      OUTPUTS
           Z             ß            /
           Y             µ0            µS
           X             f0            fS
 STACK(SF00)        INPUTS      OUTPUTS
           Z             ß            /
           Y             µS            µ0
           X             fS            f0

Example1:      ß = 0.7     µ0 = 41°     f0 = 10

   CF 00

   0.7  ENTER^
   41  ENTER^
   10  XEQ "DAB"   yields     fS = 21.4004   X<>Y   µS =  17.8522°

-Likewise,   ß = 0.7     µ0 = 150°     f0 = 10   give   fS = 5.5141  and   µS =  114.9367°

Example2:      ß = 0.7     µS = 41°     fS = 10

   SF 00

   0.7  ENTER^
   41  ENTER^
   10  XEQ "DAB"  gives in 1.8 second   f0 = 6.6052   X<>Y   µ0 =  83.3397°

-Likewise,   ß = 0.7     µS = 150°     fS = 10   produce   f0 = 22.4915  and   µ0 =  167.1555°

Notes:

-The frequencies  f0 & fS are equal  iff  ( ß = 0  or  ß cos2 µ0 + 2 cos µ0 + ß = 0 )
                                                     i-e  ( ß = 0  or  ß cos2 µS - 2 cos µS + ß = 0 )

      for  ß = 0.7  the solution is  µ0 = 114.1023165°  i-e  µS =  65.89768346°

-"DAB" works in all angular modes.
-However, the program listed above should not be used  if  µ  is very near 0° or 180°  because ACOS doesn't give enough accurate results in that cases.
-For example, with  ß = 0.7 and  µ0 = 0.01°  it yields  µS = 0.003969568048°  whereas the correct value is µS = 0.004200840261°

-The following variant uses R-P instead of ACOS , it employs the formulae:

   Tan µ0 = Sin µS ( cos µS - ß ) -1 ( 1 - ß2 ) 1/2
   Tan µS = Sin µ0 ( cos µ0 + ß ) -1 ( 1 - ß2 ) 1/2
 

01  LBL "DAB"
02  X<> Z
03  FS? 00
04  CHS
05  STO T
06  X^2
07  SIGN
08  ST- L
09  X<> L
10  CHS
11  SQRT
12  ST/ Z
13  X<>Y
14  SIN
15  ST* Y
16  X<> L
17  COS
18  STO L
19  R^
20  ST* Y
21  ST+ L
22  CLX
23  1
24  ST+ Y
25  RDN
26  ST* Z
27  X<> L
28  R-P
29  X<> Z
30  END

( 53 bytes / SIZE 000 )
 
 
 STACK(CF 00)        INPUTS      OUTPUTS
           Z             ß            /
           Y             µ0            µS
           X             f0            fS
 STACK(SF00)        INPUTS      OUTPUTS
           Z             ß            /
           Y             µS            µ0
           X             fS            f0

-Examples 1 and 2 produce the same results, but this time,  ß = 0.7 and  µ0 = 0.01°  yield  µS = 0.004200840259°
-Moreover, this program also works if  -180° < µ < 0°

-6 bytes may be saved if you replace lines 06 to 11 by  ASIN  COS  but the execution time becomes 3.4 seconds instead of 2.5 seconds.
 

References:

[1]  Robin M. Green - "Spherical Astronomy" - Cambridge University Press - ISBN 0-521-31779-7
[2]  Stamatia Mavridès - "L'Univers relativiste" - Masson  ISBN 2-225-36080-7  ( in French )
 
 

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