#!/usr/bin/env python # -*- coding: utf-8 -*- # (c) G. Beyerle 2016-2017 # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # You should have received a copy of the GNU General Public License # along with this program. If not, see . # This script requires python 2; for python3 the call to # raw_input() needs to be replaced by input() # Disclaimer: I'm a SymPy newbie, please bear with me. # Changes # 170507 - correct typo # Notes & Comments # # Consider five boosts, # instead of Cos() and sqrt(1-Cos()**2) # use half-angle parameterization # S = sin( alpha) # C = cos( alpha) # with # T = tan( alpha/2) # we get # S = 2 * T / (1 + T**2); # C = (1 - T**2) / (1 + T**2); # # Boosts # t = c/g * (sinh(g*tau/c)) # x = c^2/g * (cosh(g*tau/c)-1) # v = c * tanh(g*tau/c) # # with c = 1 and g = 1 we obtain # # t = (sinh(tau)) = beta*gamma # x = (cosh(tau)-1) = (gamma-1) # v = beta # # t_half = sinh(tau/2) # = sqrt((cosh(tau) - 1)/2) # = sqrt((gamma - 1)/2) # x_half = (cosh(tau/2) - 1) # = sqrt((cosh(tau) + 1)/2) - 1 # = sqrt((gamma + 1)/2) - 1 # v_half = tanh(tau/2) # = sqrt((gamma-1)/(gamma+1)) # # Simplifications / assumptions : # => 2-d geometry, i.e. all boosts are in one plane # => first boost is in [Ca, Sa] direction # Ca = C12*C23 - S12*S23; # Sa = -S12*C23 - C12*S23; # => second boost is in [Cb, Sb] direction # Cb = C23; # Sb = -S23; # => third boost is in x direction [1,0] # => 3rd boost in +x direction # i.e. e1 = [Ca,Sa] # e2 = [Cb,Sb] # e3 = [1,0] # e4 = [Cb,-Sb] # e5 = [Ca,-Sa] # and Ca/Sa/Cb/Sb are expressed in terms of C12/S12/C23/S23 # Ca = C12*C23 - S12*S23; # Sa = -S12*C23 - C12*S23; # Cb = C23; # Sb = S23; # => all angles are arbitrary # => time reversal symmetry, # i.e. symmetry with respect to 3rd boost, # i.e. C14 = C25, C23 = C34, C13 = C35 # => 3rd boost in x direction : # i.e. e1 = [Ca,Sa] # e2 = [Cb,Sb] # e3 = [1,0] # e4 = [Cb,-Sb] # e5 = [Ca,-Sa] # and Ca/Sa/Cb/Sb are expressed in terms of C12/S12/C23/S23 # Ca = C12*C23 - S12*S23; # Sa = -S12*C23 - C12*S23; # Cb = C23; # Sb = S23; # # time reversal symmetry implies that vec(B->E) = a * vec(C->D) # thus # (B->E)_x = a * (C->D)_x # (B->E)_y = a * (C->D)_y # thus # (B->E)_x / (B->E)_y = (C->D)_x / (C->D)_y # thus # (B->E)_x * (C->D)_y - (C->D)_x * (B->E)_y = 0 # # Imports import sys from sympy import * from sympy import sin,cos,tan from sympy import sinh,cosh,tanh from sympy import symbols from sympy import MatrixSymbol, Matrix # # helper functions # def press_enter() : raw_input(" ... press ENTER to continue ...") # print("*** debug *** Press ENTER to continue ...") # In python3 use function input() instead # w = g * tau / c # g : acceleration # tau : proper time # c : velocity of light # half-angle parameterization # cos(zeta) = (1-tan(zeta/2)^2) / (1+tan(zeta/2)^2) # sin(zeta) = (2*tan(zeta/2)) / (1+tan(zeta/2)^2) # T12 = tan(zeta12/2) # T23 = tan(zeta23/2) # C13 = (1 - T13**2) / (1 + T13**2); # C23 = (1 - T23**2) / (1 + T23**2); # S13 = 2 * T13 / (1 + T13**2); # S23 = 2 * T23 / (1 + T23**2); # gamma : Lorentz factor # ga = 1 / sqrt(1- v**2/c**2) # with speed of light c and speed v ga = Symbol('ga') # normalized speed # be = v/c be = sqrt(ga**2 - 1) / ga; T12 = symbols( 'T12') T23 = symbols( 'T23') # half-angle parameterization S12 = 2 * T12 / (1 + T12**2); S23 = 2 * T23 / (1 + T23**2); C12 = (1 - T12**2) / (1 + T12**2); C23 = (1 - T23**2) / (1 + T23**2); Ca = C12*C23 - S12*S23; Sa = -S12*C23 - C12*S23; Cb = C23; Sb = -S23; # 3rd boost in frame [3] is taken to be along x-axis # e3 = [1,0] # thus, 2nd boost in frame [2] is along # e2 = [cos(zeta2),sin(zeta2)] # = [Cb,Sb] # and, similarly, # e4 = [cos(zeta2),-sin(zeta2)] # = [Cb,-Sb] # # i.e. e1 = [Ca,Sa] # e2 = [Cb,Sb] # e3 = [1,0] # e4 = [Cb,-Sb] # e5 = [Ca,-Sa] print ('\nWarning: Script execution time may be minutes to hours ...\n') press_enter() print ('\nWe simplify the problem and consider a 1+2 dimensional,\n' 'instead of 1+3 dimensional space-time. This is possible\n' 'since the spatial motion is restricted to the xy-plane.\n\n' 'First we define the boost directions.\n\n' 'The 3rd boost is assumed to be along x-axis, i.e.:\n' ' e3 = [1,0]\n\n' 'The second and fourth boost directions are written as\n' ' e2 = [cos(zetaB), sin(zetaB)] = [cos(zeta2), sin(zeta2)] \n' ' e4 = [Cb,-Sb] = [cos(zeta2),-sin(zeta2)] \n' 'where zeta2 denotes the angle between boost \n' 'directions e2 and e3. \n\n' 'Similarly the first and fifth boost are \n' ' e1 = [Ca, Sa] = [cos(zeta1), sin(zeta1)] \n' ' e5 = [Ca,-Sa] = [cos(zeta1),-sin(zeta1)] \n' 'where \n' ' sin(zetaA) = -sin(zeta12)*cos(zeta23) - cos(zeta12)*sin(zeta23) \n' ' cos(zetaA) = cos(zeta12)*cos(zeta23) - sin(zeta12)*sin(zeta23) \n' ' sin(zetaB) = -sin(zeta12) \n' ' cos(zetaB) = cos(zeta12) \n' 'Here \n' ' zeta12 (zeta23) denotes the angle between boost unit \n' 'vector e1 and e2 (e2 and e3). For symmetry reasons zeta12 (zeta23) \n' 'is also the angle between boost vectors e4 and e5 (e3 and e4). \n') print 'With the half-angle parametrization this translates to' print ' e1 = ', [Ca,Sa] print ' e2 = ', [Cb,Sb] print ' e3 = ', [1,0] print ' e4 = ', [Cb,-Sb] print ' e5 = ', [Ca,-Sa] print ' ' press_enter() # 1st boost in [Ca, Sa] direction LT2To1 = Matrix([ [ ga, ga*(+be)*Ca, ga*(+be)*Sa ], [ ga*(+be)*Ca, 1+(ga-1)*Ca**2, (ga-1)*Sa*Ca ], [ ga*(+be)*Sa, (ga-1)*Sa*Ca, 1+(ga-1)*Sa**2 ] ]); # inverse 1st boost in [Ca, Sa] direction LT1To2 = Matrix([ [ ga, ga*(-be)*Ca, ga*(-be)*Sa ], [ ga*(-be)*Ca, 1+(ga-1)*Ca**2, (ga-1)*Sa*Ca ], [ ga*(-be)*Sa, (ga-1)*Sa*Ca, 1+(ga-1)*Sa**2 ] ]); # 2nd boost in [Cb, Sb] direction LT3To2 = Matrix([ [ ga, ga*(+be)*Cb, ga*(+be)*Sb ], [ ga*(+be)*Cb, 1+(ga-1)*Cb**2, (ga-1)*Sb*Cb ], [ ga*(+be)*Sb, (ga-1)*Sb*Cb, 1+(ga-1)*Sb**2 ] ]); # inverse 2nd boost in [Cb, Sb] direction LT2To3 = Matrix([ [ ga, ga*(-be)*Cb, ga*(-be)*Sb ], [ ga*(-be)*Cb, 1+(ga-1)*Cb**2, (ga-1)*Sb*Cb ], [ ga*(-be)*Sb, (ga-1)*Sb*Cb, 1+(ga-1)*Sb**2 ] ]); # 3rd boost in [1, 0] direction LT4To3 = Matrix([ [ ga, ga*(+be), 0 ], [ ga*(+be), ga, 0 ], [ 0, 0, 1 ] ]); # boost in [1, 0] direction LT3To4 = Matrix([ [ ga, ga*(-be), 0 ], [ ga*(-be), ga, 0 ], [ 0, 0, 1 ] ]); # boost in [Cb, -Sb] direction LT5To4 = Matrix([ [ ga, ga*(+be)*Cb, ga*(+be)*(-Sb) ], [ ga*(+be)*Cb, 1+(ga-1)*Cb**2, (ga-1)*(-Sb)*Cb ], [ ga*(+be)*(-Sb), (ga-1)*(-Sb)*Cb, 1+(ga-1)*Sb**2 ] ]); # boost in [Cb, -Sb] direction LT4To5 = Matrix([ [ ga, ga*(-be)*Cb, ga*(-be)*(-Sb) ], [ ga*(-be)*Cb, 1+(ga-1)*Cb**2, (ga-1)*(-Sb)*Cb ], [ ga*(-be)*(-Sb), (ga-1)*(-Sb)*Cb, 1+(ga-1)*Sb**2 ] ]); # boost in [Ca, -Sa] direction LT6To5 = Matrix([ [ ga, ga*(+be)*Ca, ga*(+be)*(-Sa) ], [ ga*(+be)*Ca, 1+(ga-1)*Ca**2, (ga-1)*(-Sa)*Ca ], [ ga*(+be)*(-Sa), (ga-1)*(-Sa)*Ca, 1+(ga-1)*Sa**2 ] ]); # boost in [Ca, -Sa] direction LT5To6 = Matrix([ [ ga, ga*(-be)*Ca, ga*(-be)*(-Sa) ], [ ga*(-be)*Ca, 1+(ga-1)*Ca**2, (ga-1)*(-Sa)*Ca ], [ ga*(-be)*(-Sa), (ga-1)*(-Sa)*Ca, 1+(ga-1)*Sa**2 ] ]); print ('\nIn 1+2 dimensional space-time the Lorentz \n' 'transformation matrices reduce to 3x3 matrices.\n\n' 'Transformation from frame [2] to [1] ... \n') print ' LT2To1 = ', LT2To1, '\n' press_enter() print '\nTransformation from frame [3] to [2] ... \n' print ' LT3To2 = ', LT3To2, '\n' press_enter() print '\nTransformation from frame [4] to [3] ... \n' print ' LT4To3 = ', LT4To3, '\n' press_enter() print '\nTransformation from frame [5] to [4] ... \n' print ' LT5To4 = ', LT5To4, '\n' press_enter() print '\nTransformation from frame [6] to [5] ... \n' print ' LT6To5 = ', LT6To5, '\n' press_enter() print ('\nTo transform the four-velocity at event F from \n' 'frame [6] to frame [1], i.e. \n\n' ' vlF_Fr1 = LT6To1 * vlF_Fr6\n\n' 'the corresponding Lorentz transformation matrix.\n') # eq = LT5To6 * LT4To5 * LT3To4 * LT2To3 * LT1To2; vlF_Fr6 = Matrix([ [1], [0], [0] ]); print ' vlF_Fr5 = LT6To5 * vlF_Fr6 ' vlF_Fr5 = factor( expand( LT6To5 * vlF_Fr6)); print ' vlF_Fr4 = LT5To4 * vlF_Fr5 ' vlF_Fr4 = factor( expand( LT5To4 * vlF_Fr5)); print ' vlF_Fr3 = LT4To3 * vlF_Fr4 ' vlF_Fr3 = factor( expand( LT4To3 * vlF_Fr4)); print ' vlF_Fr2 = LT3To2 * vlF_Fr3 ' vlF_Fr2 = factor( expand( LT3To2 * vlF_Fr3)); print ' vlF_Fr1 = LT2To1 * vlF_Fr2 ' vlF_Fr1 = factor( expand( LT2To1 * vlF_Fr2)); print ('\nNow evaluate\n' ' vlF_Fr1[0] = 1\n' 'or\n' ' vlF_Fr1[0] - 1 = 0\n\n' 'The result is\n') print ' vlF_Fr1[0] - 1 = ', factor( expand( vlF_Fr1[0] - 1)), ' = 0\n' press_enter() print '\nAssuming ga =|= 1 we find \n' # remove non-zero factors tmpRes = factor((vlF_Fr1[0] - 1)*(T12**2 + 1)**2*(T23**2 + 1)**2/(ga-1)) eq1 = (T12**2*T23**2 + T12**2 - 4*T12*T23*ga - 4*T12*T23 - 2*T23**2*ga - T23**2 + 4*ga**2 + 2*ga - 1); # make sure I copied/pasted this correctly assert eq1**2 == tmpRes print eq1, ' = 0' print ('\n(This is equation 29 in reference {1}.) \n\n' 'Now solve this for T23 ... \n') press_enter() eq1Pol = Poly( eq1, T23) eq1PlCf = eq1Pol.all_coeffs(); #print ' eq1PolCoef = ', eq1PolCoef #sys.exit(); # eq1 is a second-order polynomial in T23 #res = solve( collect( eq1, T23), T23 ); # solution of second order polynomial # a*x**2 + b*x + c = 0 # is # xa = ((-b + sqrt(-4*a*c + b**2))/(2*a)); # xb = ((-b - sqrt(-4*a*c + b**2))/(2*a)); sqrtArg = factor( expand( -4*eq1PlCf[0]*eq1PlCf[2] + eq1PlCf[1]**2)); T23a = ((-eq1PlCf[1] + sqrt( sqrtArg)) / (2*eq1PlCf[0])); T23b = ((-eq1PlCf[1] - sqrt( sqrtArg)) / (2*eq1PlCf[0])); # checks assert factor(expand(eq1PlCf[0]*T23a**2 + eq1PlCf[1]*T23a + eq1PlCf[2])) == 0 assert factor(expand(eq1PlCf[0]*T23b**2 + eq1PlCf[1]*T23b + eq1PlCf[2])) == 0 print ('\nSince this polynomial is quadratic in T23 we get \n' 'two solutions \n') print ' T23a = ', T23a print '\nand\n' print ' T23b = ', T23b, '\n' press_enter() # # calculate four-position of event F in frame [1] # evA_Fr1 = Matrix([ [0], [0], [0] ]); # new position := old position + boost distance evB_Fr1 = ( evA_Fr1 + Matrix([ [ga*be], [(ga-1)*(+Ca)], [(ga-1)*(+Sa)] ]) ); print ('\nDetermine sequentially evB_Fr2, evC_Fr3, ..., evF_Fr6 \n\n' 'Transform event B from frame [1] to frame [2].\n\n' ' evB_Fr2 = LT1To2 * evB_Fr1 \n') evB_Fr2 = factor( expand( LT1To2 * evB_Fr1 )) print ' evB_Fr2 = ', evB_Fr2 print ('\nMove from event B to event C in frame [2].\n\n' ' evC_Fr2 = ( evB_Fr2 \n' ' + Matrix([ [ga*be], \n' ' [(ga-1)*(+Cb)], \n' ' [(ga-1)*(+Sb)] ]) ); \n' ) evC_Fr2 = ( evB_Fr2 + Matrix([ [ga*be], [(ga-1)*(+Cb)], [(ga-1)*(+Sb)] ]) ); print ' evC_Fr2 = ', evC_Fr2, '\n' press_enter() print ('\nTransform event C from frame [2] to frame [3].\n\n' ' evC_Fr3 = LT2To3 * evC_Fr2 \n') evC_Fr3 = factor( expand( LT2To3 * evC_Fr2 )) print ' evC_Fr3 = ', evC_Fr3 print ('\nMove from event C to event D in frame [3].\n\n' ' evD_Fr3 = ( evC_Fr3 \n' ' + Matrix([ [ga*be], \n' ' [(ga-1)], \n' ' [0] ]) ); \n' ) evD_Fr3 = ( evC_Fr3 + Matrix([ [ga*be], [(ga-1)], [0] ]) ); print ' evD_Fr3 = ', evD_Fr3, '\n' press_enter() print ('\nTransform event D from frame [3] to frame [4].\n\n' ' evD_Fr4 = LT3To4 * evD_Fr3 \n') evD_Fr4 = factor( expand( LT3To4 * evD_Fr3 )) print ' evD_Fr4 = ', evD_Fr4 print ('\nMove from event D to event E in frame [4].\n\n' ' evE_Fr4 = ( evD_Fr4 \n' ' + Matrix([ [ga*be], \n' ' [(ga-1)*(+Cb)],\n' ' [(ga-1)*(-Sb)] ]) ); \n' ) evE_Fr4 = ( evD_Fr4 + Matrix([ [ga*be], [(ga-1)*(+Cb)], [(ga-1)*(-Sb)] ]) ); print ' evE_Fr4 = ', evE_Fr4, '\n' press_enter() print ('\nTransform event E from frame [4] to frame [5].\n\n' ' evE_Fr5 = LT4To5 * evE_Fr4 \n') evE_Fr5 = factor( expand( LT4To5 * evE_Fr4 )) print ' evE_Fr5 = ', evE_Fr5 print ('\nMove from event E to event F in frame [5].\n\n' ' evF_Fr5 = ( evE_Fr5 \n' ' + Matrix([ [ga*be], \n' ' [(ga-1)*(+Ca)], \n' ' [(ga-1)*(-Sa)] ]) ); \n' ) evF_Fr5 = ( evE_Fr5 + Matrix([ [ga*be], [(ga-1)*(+Ca)], [(ga-1)*(-Sa)] ]) ); print ' evF_Fr5 = ', evF_Fr5, '\n' press_enter() print ('\nTransform event F from frame [5] to frame [6].\n\n' ' evF_Fr6 = LT5To6 * evF_Fr5 \n') evF_Fr6 = factor( expand( LT5To6 * evF_Fr5 )) print '\nThe final result is\n\n' print ' evF_Fr6 = ', evF_Fr6 #print ' ' #print ' evF_Fr6[1]-evF_Fr6[2] = ', factor( expand( evF_Fr6[1]-evF_Fr6[2])) print ('\n\nNow things get a bit messy ... \n\n' 'We use the condition\n\n' ' evF_Fr6[1] = 0 = evF_Fr6[2] \n\n' 'Start with evF_Fr6[2] = 0 , insert T23 = T23(T12,ga) \n' 'and assume that \n' ' (T12**2 + 1) =|= 0 \n' 'and \n' ' (T23**2 + 1) =|= 0 \n') tmpExpr = factor( evF_Fr6[2] * ((T12**2 + 1)**3 *(T23**2 + 1)**3) / (4*(ga - 1))) print ('\na) \n' 'Insert first solution T23a = T23a(T12,ga) of quadratic \n' 'polynomial and collect all terms containing square \n' 'root appearing in the solution T23 = T23(T12,ga) : \n\n') expr1 = collect( factor( expand( tmpExpr.subs(T23, T23a))), sqrt( sqrtArg)); print ' 0 = ', expr1 print ('\n\nThis expression has the form \n' ' (-8)*(T12**2 + 1)*(ga + 1)**2 \n' ' * (expr1 + expr2 * sqrt(expr3)) \n' ' / (T12**2 - 2*ga - 1)**6 = 0 \n' 'We drop the non-zero factors \n' ' (-8)*(T12**2 + 1)*(ga + 1)**2 \n' 'and \n' ' (T12**2 - 2*ga - 1)**6 \n' 'move the second term to the RHS, square both sides, \n' 'move everything back to the LHS. \n') # Don't know how to do this programmatically in SymPy ... expr2 = ((14*T12**13*ga + 10*T12**13 + 8*T12**11*ga**3 + 8*T12**11*ga**2 + 76*T12**11*ga + 52*T12**11 - 224*T12**9*ga**4 - 1320*T12**9*ga**3 - 2904*T12**9*ga**2 - 2502*T12**9*ga - 754*T12**9 - 1184*T12**7*ga**5 - 5760*T12**7*ga**4 - 10224*T12**7*ga**3 - 8432*T12**7*ga**2 - 2904*T12**7*ga - 200*T12**7 - 2112*T12**5*ga**6 - 8704*T12**5*ga**5 - 13600*T12**5*ga**4 - 9808*T12**5*ga**3 - 3568*T12**5*ga**2 - 862*T12**5*ga - 154*T12**5 - 1536*T12**3*ga**7 - 5376*T12**3*ga**6 - 7456*T12**3*ga**5 - 5568*T12**3*ga**4 - 2456*T12**3*ga**3 - 536*T12**3*ga**2 + 12*T12**3*ga + 20*T12**3 - 256*T12*ga**8 - 1024*T12*ga**7 - 1728*T12*ga**6 - 1472*T12*ga**5 - 576*T12*ga**4 - 8*T12*ga**3 + 72*T12*ga**2 + 22*T12*ga + 2*T12)**2 - (-T12**4 + 8*T12**2*ga + 6*T12**2 + 8*ga**3 + 8*ga**2 - 1) *(T12**12*ga + 3*T12**12 - 12*T12**10*ga**2 - 28*T12**10*ga - 8*T12**10 - 216*T12**8*ga**3 - 788*T12**8*ga**2 - 899*T12**8*ga - 317*T12**8 - 640*T12**6*ga**4 - 2112*T12**6*ga**3 - 2568*T12**6*ga**2 - 1216*T12**6*ga - 120*T12**6 - 688*T12**4*ga**5 - 2256*T12**4*ga**4 - 2512*T12**4*ga**3 - 1272*T12**4*ga**2 - 389*T12**4*ga - 71*T12**4 - 320*T12**2*ga**6 - 832*T12**2*ga**5 - 960*T12**2*ga**4 - 640*T12**2*ga**3 - 236*T12**2*ga**2 - 36*T12**2*ga - 64*ga**6 - 144*ga**5 - 112*ga**4 - 24*ga**3 + 12*ga**2 + 7*ga + 1)**2) print '\n 0 = ', expr2 print '\nEvaluate this and obtain\n' expr3 = factor(expand( expr2)) print ' 0 = ', expr3 print ('\n\nwhich has the following two solutions :\n' '(The expression T12**2 - 2*ga - 1 is not zero, since\n' 'it appeared in the denominator before and was therefore\n' 'assumed not to evaluate to zero.) \n') expr3a = collect( T12**8 + 4*T12**6*ga + 8*T12**6 - 8*T12**4*ga**3 - 36*T12**4*ga**2 - 52*T12**4*ga - 18*T12**4 - 32*T12**2*ga**4 - 112*T12**2*ga**3 - 104*T12**2*ga**2 - 20*T12**2*ga + 8*T12**2 - 32*ga**5 - 64*ga**4 - 40*ga**3 - 4*ga**2 + 4*ga + 1, T12) expr3b = collect( T12**8*ga**2 + 6*T12**8*ga + 9*T12**8 - 24*T12**6*ga**3 + 60*T12**6*ga**2 + 48*T12**6*ga - 20*T12**6 - 48*T12**4*ga**4 + 88*T12**4*ga**3 + 110*T12**4*ga**2 - 52*T12**4*ga - 2*T12**4 - 32*T12**2*ga**5 + 64*T12**2*ga**4 + 56*T12**2*ga**3 - 20*T12**2*ga**2 - 4*T12**2 + 16*ga**4 + 8*ga**3 - 7*ga**2 - 2*ga + 1, T12) print ' solution (1) : ', expr3a, ' = 0' print 'and' print ' solution (2) : ', expr3b, ' = 0\n' press_enter() print ('\nb) \n' 'Insert second solution T23b = T23b(T12,ga) of quadratic \n' 'polynomial and collect all terms containing square \n' 'root appearing in the solution T23 = T23(T12,ga) : \n') expr4 = collect( factor( expand( tmpExpr.subs(T23, T23b))), sqrt( sqrtArg)); print ' 0 = ', expr4 print ('\nThis expression has the form \n' ' (-8)*(T12**2 + 1)*(ga + 1)**2 \n' ' * (expr1 + expr2 * sqrt(expr3)) \n' ' / (T12**2 - 2*ga - 1)**6 = 0 \n' 'We drop the non-zero factors \n' ' (-8)*(T12**2 + 1)*(ga + 1)**2 \n' 'and \n' ' (T12**2 - 2*ga - 1)**6 \n' 'move the second term to the RHS, square both sides, \n' 'move everything back to the LHS. \n') expr5 = ((14*T12**13*ga + 10*T12**13 + 8*T12**11*ga**3 + 8*T12**11*ga**2 + 76*T12**11*ga + 52*T12**11 - 224*T12**9*ga**4 - 1320*T12**9*ga**3 - 2904*T12**9*ga**2 - 2502*T12**9*ga - 754*T12**9 - 1184*T12**7*ga**5 - 5760*T12**7*ga**4 - 10224*T12**7*ga**3 - 8432*T12**7*ga**2 - 2904*T12**7*ga - 200*T12**7 - 2112*T12**5*ga**6 - 8704*T12**5*ga**5 - 13600*T12**5*ga**4 - 9808*T12**5*ga**3 - 3568*T12**5*ga**2 - 862*T12**5*ga - 154*T12**5 - 1536*T12**3*ga**7 - 5376*T12**3*ga**6 - 7456*T12**3*ga**5 - 5568*T12**3*ga**4 - 2456*T12**3*ga**3 - 536*T12**3*ga**2 + 12*T12**3*ga + 20*T12**3 - 256*T12*ga**8 - 1024*T12*ga**7 - 1728*T12*ga**6 - 1472*T12*ga**5 - 576*T12*ga**4 - 8*T12*ga**3 + 72*T12*ga**2 + 22*T12*ga + 2*T12)**2 -(-T12**4 + 8*T12**2*ga + 6*T12**2 + 8*ga**3 + 8*ga**2 - 1) *(-T12**12*ga - 3*T12**12 + 12*T12**10*ga**2 + 28*T12**10*ga + 8*T12**10 + 216*T12**8*ga**3 + 788*T12**8*ga**2 + 899*T12**8*ga + 317*T12**8 + 640*T12**6*ga**4 + 2112*T12**6*ga**3 + 2568*T12**6*ga**2 + 1216*T12**6*ga + 120*T12**6 + 688*T12**4*ga**5 + 2256*T12**4*ga**4 + 2512*T12**4*ga**3 + 1272*T12**4*ga**2 + 389*T12**4*ga + 71*T12**4 + 320*T12**2*ga**6 + 832*T12**2*ga**5 + 960*T12**2*ga**4 + 640*T12**2*ga**3 + 236*T12**2*ga**2 + 36*T12**2*ga + 64*ga**6 + 144*ga**5 + 112*ga**4 + 24*ga**3 - 12*ga**2 - 7*ga - 1)**2) print '\n 0 = ', expr5 print ('\nThis is the same expression as the one obtained for T23a \n' 'and therefore does not yield new solutions. \n') # check tmpRes = expand( expr2 - expr5) #print ' tmpRes = ', tmpRes assert tmpRes == 0 print ('\nNow consider the constraint evF_Fr6[1] = 0 , insert \n' 'T23a = T23a(T12,ga) and again assume that \n' ' (T12**2 + 1) =|= 0 \n' 'and \n' ' (T23**2 + 1) =|= 0 \n') tmpExpr = factor( evF_Fr6[1] * ((T12**2 + 1)**3 *(T23**2 + 1)**3) / ( -(ga - 1))) #print ' 0 = ', tmpExpr press_enter() print ('\nc)\n' 'Insert first solution T23a = T23a(T12,ga) of quadratic \n' 'polynomial and collect all terms containing square \n' 'root appearing in the solution T23 = T23(T12,ga) :') expr6 = collect( factor( expand( tmpExpr.subs(T23, T23a))), sqrt( sqrtArg)); print ' 0 = ', expr6 print ('\nThis expression has the form \n' ' (-32)*(T12**2 + 1)*(ga - 1)*(ga + 1)**2 \n' ' * (expr1 + expr2 * sqrt(expr3)) \n' ' / (T12**2 - 2*ga - 1)**6 = 0 \n' 'We drop the non-zero factors \n' ' (-32)*(T12**2 + 1)*(ga - 1)*(ga + 1)**2 \n' 'and \n' ' (T12**2 - 2*ga - 1)**6 \n' 'move the second term to the RHS, square both sides, \n' 'move everything back to the LHS. \n') expr7 = ((T12**14 + 8*T12**12*ga + 13*T12**12 - 8*T12**10*ga**3 - 64*T12**10*ga**2 - 64*T12**10*ga + T12**10 - 256*T12**8*ga**4 - 1464*T12**8*ga**3 - 3248*T12**8*ga**2 - 3032*T12**8*ga - 987*T12**8 - 960*T12**6*ga**5 - 4400*T12**6*ga**4 - 6928*T12**6*ga**3 - 4768*T12**6*ga**2 - 1408*T12**6*ga - 133*T12**6 - 1024*T12**4*ga**6 - 3264*T12**4*ga**5 - 3952*T12**4*ga**4 - 1904*T12**4*ga**3 + 64*T12**4*ga**2 + 344*T12**4*ga + 79*T12**4 - 128*T12**2*ga**7 - 384*T12**2*ga**6 - 64*T12**2*ga**5 + 688*T12**2*ga**4 + 792*T12**2*ga**3 + 352*T12**2*ga**2 + 64*T12**2*ga + 3*T12**2 + 128*ga**7 + 384*ga**6 + 448*ga**5 + 240*ga**4 + 40*ga**3 - 16*ga**2 - 8*ga - 1)**2 - (-T12**4 + 8*T12**2*ga + 6*T12**2 + 8*ga**3 + 8*ga**2 - 1) *(2*T12**11 - 24*T12**9*ga**2 - 72*T12**9*ga - 38*T12**9 - 208*T12**7*ga**3 - 768*T12**7*ga**2 - 992*T12**7*ga - 412*T12**7 - 416*T12**5*ga**4 - 1296*T12**5*ga**3 - 1360*T12**5*ga**2 - 592*T12**5*ga - 92*T12**5 - 192*T12**3*ga**5 - 480*T12**3*ga**4 - 368*T12**3*ga**3 + 96*T12**3*ga + 26*T12**3 + 64*T12*ga**5 + 192*T12*ga**4 + 208*T12*ga**3 + 104*T12*ga**2 + 24*T12*ga + 2*T12)**2) print '\n 0 = ', expr7 print '\nEvaluate this and obtain' expr8 = factor(expand( expr7)) print ' 0 = ', expr8 print ('\nwhich has the following two solutions \n' '(The factor T12**2 - 2*ga - 1 is non-zero.) \n') expr8a = collect(T12**8 + 24*T12**6*ga + 28*T12**6 - 8*T12**4*ga**3 + 56*T12**4*ga**2 + 48*T12**4*ga - 10*T12**4 + 16*T12**2*ga**3 + 16*T12**2*ga**2 - 8*T12**2*ga - 4*T12**2 - 8*ga**3 - 8*ga**2 + 1, T12); expr8b = collect(T12**8 + 4*T12**6*ga + 8*T12**6 - 8*T12**4*ga**3 - 36*T12**4*ga**2 - 52*T12**4*ga - 18*T12**4 - 32*T12**2*ga**4 - 112*T12**2*ga**3 - 104*T12**2*ga**2 - 20*T12**2*ga + 8*T12**2 - 32*ga**5 - 64*ga**4 - 40*ga**3 - 4*ga**2 + 4*ga + 1, T12); print ' solution (1) : ', expr8a, ' = 0' print 'and' print ' solution (2) : ', expr8b, ' = 0\n' press_enter() print ('\nFinally, we determine the Thomas-Wigner rotation angle thetaTW.\n' 'By construction the Lorentz transformation matrix connecting\n' 'frame [1] and frame [6] is a rotation in the xy-plane,\n' ' [1, 0, 0, 0]\n' ' LT1To6 = [0, R11, R12, 0]\n' ' [0, R21, R22, 0]\n' ' [0, 0, 0, 1]\n' 'Thus, \n' ' thetaTW = atan(R21,R11)\n' 'LT1To6 \n'); LT1To3 = factor( expand( LT2To3 * LT1To2)); print '\n LT1To3 = ', LT1To3, '\n' LT1To4 = factor( expand( LT3To4 * LT1To3)); print '\n LT1To4 = ', LT1To4, '\n' LT4To6 = factor( expand( LT5To6 * LT4To5)); print '\n LT4To6 = ', LT4To6, '\n' LT1To6 = factor( expand( LT4To6 * LT1To4)); print '\n LT1To6 = ', LT1To6, '\n' print ' LT1To6[1,1] + 1 = R11 + 1 = ', factor( expand( LT1To6[1,1] + 1)), '\n\n' print ' LT1To6[2,1] = R21 + 1 = ', LT1To6[2,1], '\n' print '\n ... q.e.d.' sys.exit(); # # end of file #