lorentz.R
In lorentz: The Lorentz Transform in Relativistic Physics

### R code from vignette source 'lorentz.Rnw'

###################################################
### code chunk number 1: lorentz.Rnw:87-88
###################################################
library("lorentz")


###################################################
### code chunk number 2: floater
###################################################
                          # NB: speed of light = 1 by default
u <- as.3vel(c(0.6,0,0))  # coerce to a three-velocity
u
as.4vel(u)                # four-velocity is better for calculations
(B <- boost(u))           # transformation matrix


###################################################
### code chunk number 3: flighter
###################################################
B %*% c(1,0,0,0)


###################################################
### code chunk number 4: flayter
###################################################
B %*% c(4,6,-8,9)


###################################################
### code chunk number 5: flooter
###################################################
B %*% c(1,1,0,0)
B %*% c(1,-1,0,0)


###################################################
### code chunk number 6: lorentz.Rnw:289-291
###################################################
(B <- boost(as.3vel(c(0.8,0,0))))  # 0.8c left to right
solve(B) %*% c(1,0,0,0)            # active transform


###################################################
### code chunk number 7: lorentz.Rnw:300-304
###################################################
u <- as.3vel(c(0.3,-0.4,+0.8))
v <- as.3vel(c(0.4,+0.2,-0.1))
L <- boost(u) %*% boost(v)
L


###################################################
### code chunk number 8: lorentz.Rnw:314-316
###################################################
(P <- pureboost(L))  # pure boost
P - t(P)   # check for symmetry


###################################################
### code chunk number 9: lorentz.Rnw:321-325
###################################################
(U <- orthog(L))                  # rotation matrix
U[2:4,2:4]                        # inspect the spatial components
round(crossprod(U) - diag(4),10)  # check for orthogonality
## zero to within numerical uncertainty


###################################################
### code chunk number 10: lorentz.Rnw:335-337
###################################################
sol(299792458)
sol()


###################################################
### code chunk number 11: lorentz.Rnw:345-347
###################################################
u <- as.3vel(c(100,200,300))
as.4vel(u)


###################################################
### code chunk number 12: lorentz.Rnw:354-355
###################################################
gam(u)


###################################################
### code chunk number 13: lorentz.Rnw:360-361
###################################################
gam(u)-1


###################################################
### code chunk number 14: lorentz.Rnw:367-368
###################################################
gamm1(u)


###################################################
### code chunk number 15: lorentz.Rnw:373-374
###################################################
boost(u)


###################################################
### code chunk number 16: lorentz.Rnw:385-387
###################################################
sol(299792458)
disp  <- c(1,1,0,0)


###################################################
### code chunk number 17: lorentz.Rnw:395-396
###################################################
ptm(to_natural=TRUE) %*% disp


###################################################
### code chunk number 18: lorentz.Rnw:406-407
###################################################
ptm(to_natural=TRUE,change_time=FALSE) %*% disp


###################################################
### code chunk number 19: lorentz.Rnw:415-418
###################################################
sol(1)
B1 <- boost((2:4)/10) %*% boost(c(-5,1,3)/10)
orthog(B1)[2:4,2:4]


###################################################
### code chunk number 20: lorentz.Rnw:425-428
###################################################
sol(10)
B2 <- boost(2:4) %*% boost(c(-5,1,3))  # exactly the same as B1 above
orthog(B2)[2:4,2:4]


###################################################
### code chunk number 21: setcinfinite
###################################################
sol(Inf)


###################################################
### code chunk number 22: parallelograminfinitec
###################################################
u <- as.3vel(1:3)
v <- as.3vel(c(-6,8,3))
u+v


###################################################
### code chunk number 23: boostinfc
###################################################
boost(u)


###################################################
### code chunk number 24: gamm1google
###################################################
gamm1(1e100)


###################################################
### code chunk number 25: rboostinfc
###################################################
set.seed(0)
options(digits=3)
(B <- rboost(1))  # random boost, speed 1


###################################################
### code chunk number 26: orthoginfc
###################################################
orthog(B)
pureboost(B)


###################################################
### code chunk number 27: lorentz.Rnw:494-499
###################################################
sol(10)
u <- as.3vel(c(5,-6,4))
(U <- as.4vel(u))
B <- boost(U)
B %*% as.vector(U)


###################################################
### code chunk number 28: lorentz.Rnw:511-514
###################################################
u <- 1:7  # speed in the x-direction [c=10]
jj <- cbind(gam(u),gam(u)*u,0,0)
(U <- as.4vel(jj))


###################################################
### code chunk number 29: lorentz.Rnw:519-520
###################################################
(B <- boost(as.3vel(c(6,0,0))))  # 60% speed of light


###################################################
### code chunk number 30: lorentz.Rnw:530-531
###################################################
U %*% t(B)


###################################################
### code chunk number 31: lorentz.Rnw:536-537
###################################################
is.consistent.4vel(U %*% t(B))


###################################################
### code chunk number 32: lorentz.Rnw:544-545
###################################################
tcrossprod(U,B)


###################################################
### code chunk number 33: lorentz.Rnw:562-563
###################################################
solve(B)


###################################################
### code chunk number 34: lorentz.Rnw:568-569
###################################################
tcrossprod(U,solve(B))


###################################################
### code chunk number 35: lorentz.Rnw:576-577
###################################################
is.consistent.4vel(tcrossprod(U,solve(B)))


###################################################
### code chunk number 36: lorentz.Rnw:586-590
###################################################
sol(100)
B1 <- boost(r3vel(1)) %*% boost(r3vel(1))
B2 <- boost(r3vel(1)) %*% boost(r3vel(1)) 
(U <- r4vel(5))


###################################################
### code chunk number 37: lorentz.Rnw:596-600
###################################################
U %*% t(B1) %*% t(B2)
U %*% t(B2 %*% B1)    # note order of operations
tcrossprod(U, B2 %*% B1)
U %>% tcrossprod(B2 %*% B1)


###################################################
### code chunk number 38: lorentz.Rnw:610-615
###################################################
U %*% B  # Young Frankstein: Do Not Use This Brain!
## The above idiom is incorrect.  See
## https://www.youtube.com/watch?v=m7-bMBuVmHo&t=1s
## (in particular @1:08) for a technical explanation of why 
## this is a Very Bad Idea (tm).


###################################################
### code chunk number 39: lorentz.Rnw:630-633
###################################################
sol(1)        # revert to natural units 
D <- dust(1)  # Dust is the simplest nontrivial SET, with 
D             # only one nonzero component


###################################################
### code chunk number 40: lorentz.Rnw:640-642
###################################################
B <- boost(as.3vel(c(0.0,0.8,0.0)))
transform_uu(D,B)


###################################################
### code chunk number 41: lorentz.Rnw:653-656
###################################################
pg <- photongas(3)
pg
transform_uu(pg,B)


###################################################
### code chunk number 42: lorentz.Rnw:665-667
###################################################
raise(transform_dd(lower(pg),lower(B)))
raise(transform_dd(lower(pg),lower(B))) - transform_uu(pg,B) #zero to numerical precision


###################################################
### code chunk number 43: lorentz.Rnw:679-685
###################################################
B1 <- boost(as.3vel(c(0.5,-0.4,0.6)))
B2 <- boost(as.3vel(c(0.1,-0.1,0.3)))
pf <- perfectfluid(4,1)
pf
pf %>% transform_uu(B1) %>% transform_uu(B2)
pf %>% transform_uu(B2 %*% B1)  # should match


###################################################
### code chunk number 44: lorentz.Rnw:691-692
###################################################
lower(pf) %>% transform_dd(lower(B1) %*% lower(B2)) %>% raise()


###################################################
### code chunk number 45: lorentz.Rnw:704-707
###################################################
sol(10)
pf_rest <- perfectfluid(1,4)
pf_rest


###################################################
### code chunk number 46: lorentz.Rnw:715-718
###################################################
u <- as.3vel(3:5)
pf_moving <- perfectfluid(1,4,u)
pf_moving


###################################################
### code chunk number 47: lorentz.Rnw:725-726
###################################################
transform_uu(perfectfluid(1,4,u),boost(u))


###################################################
### code chunk number 48: lorentz.Rnw:736-738
###################################################
sol(1)
(A <- as.photon(as.3vel(cbind(0.9,1:5/40,5:1/40))))


###################################################
### code chunk number 49: lorentz.Rnw:746-747
###################################################
inner4(A)


###################################################
### code chunk number 50: lorentz.Rnw:754-755
###################################################
tcrossprod(A,boost(as.3vel(c(0.7,0,0))))


###################################################
### code chunk number 51: lorentz.Rnw:762-763
###################################################
tcrossprod(A,boost(as.3vel(c(-0.7,0,0))))


###################################################
### code chunk number 52: lorentz.Rnw:768-769
###################################################
tcrossprod(A,solve(boost(as.3vel(c(0.7,0,0)))))


###################################################
### code chunk number 53: lorentz.Rnw:785-787
###################################################
m <- c(1,1,1)
B <- boost(as.3vel(c(0.5,0,0)))


###################################################
### code chunk number 54: lorentz.Rnw:797-799
###################################################
A
(A <- as.4mom(A %*% t(B)))


###################################################
### code chunk number 55: lorentz.Rnw:806-807
###################################################
(A <- reflect(A,m))


###################################################
### code chunk number 56: lorentz.Rnw:816-817
###################################################
(A <- as.4mom(A %*% solve(t(B))))


###################################################
### code chunk number 57: lorentz.Rnw:826-828
###################################################
A <- as.photon(as.3vel(cbind(0.9,1:5/40,5:1/40)))
A %>% tcrossprod(B) %>% reflect(m) %>% tcrossprod(solve(B)) %>% as.4mom


###################################################
### code chunk number 58: disco_ball
###################################################
disco <- matrix(rnorm(3000),ncol=3) %>% sweep(1, sqrt(rowSums(.^2)),`/`)
head(disco)


###################################################
### code chunk number 59: lorentz.Rnw:843-845
###################################################
p <- as.photon(c(1,0,0))
reflect(p,head(disco))


###################################################
### code chunk number 60: disco_reflect_percentage
###################################################
 table(reflect(p,disco)[,2]>0) # should be TRUE with probability sqrt(0.5)


###################################################
### code chunk number 61: relativistic_disco
###################################################
B <- boost(as.3vel(c(0.5,0,0)))
p %>% tcrossprod(B) %>% reflect(head(disco)) %>% tcrossprod(solve(B))


###################################################
### code chunk number 62: lorentz.Rnw:875-886
###################################################
sol(1)
light_start <- as.photon(as.3vel(cbind(0.9,1:5/40,5:1/40)))
m <- c(1,0,0)     # mirror normal to x-axis
B1 <- boost(as.3vel(c(-0.5, 0.1, 0.0)))
B2 <- boost(as.3vel(c( 0.2, 0.0, 0.0)))
B3 <- boost(as.3vel(c( 0.0, 0.0, 0.6)))
B <- B1 %*% B2 %*% B3   # matrix multiplication is associative!
light <- light_start %*% t(B)
light <- reflect(light,m)
light <- as.4mom(light %*% solve(t(B)))
light


###################################################
### code chunk number 63: lorentz.Rnw:894-895
###################################################
light_start %>% tcrossprod(B) %>% reflect(m) %>% tcrossprod(solve(B)) %>% as.4mom


###################################################
### code chunk number 64: lorentz.Rnw:899-907
###################################################
sol(Inf)
light_start <- as.photon(as.3vel(cbind(0.9,1:5/40,5:1/40)))
B1 <- boost(as.3vel(c(-0.5, 0.1, 0.0)))
B2 <- boost(as.3vel(c( 0.2, 0.0, 0.0)))
B3 <- boost(as.3vel(c( 0.0, 0.0, 0.6)))
B <- B1 %*% B2 %*% B3
light_start
light_start %>% tcrossprod(B) %>% reflect(m) %>% tcrossprod(solve(B)) %>% as.4mom


###################################################
### code chunk number 65: kickoff
###################################################
sol(1)
u <- as.3vel(c(-0.7,+0.2,-0.3))
v <- as.3vel(c(+0.3,+0.3,+0.4))
w <- as.3vel(c(+0.1,+0.3,+0.8))
x <- as.3vel(c(-0.2,-0.1,-0.9))
u


###################################################
### code chunk number 66: try
###################################################
u+v
v+u


###################################################
### code chunk number 67: lorentz.Rnw:979-980
###################################################
(u+v)-gyr(u,v,v+u)


###################################################
### code chunk number 68: funcid
###################################################
f <- gyrfun(u,v)
(u+v)-f(v+u)    # should be zero


###################################################
### code chunk number 69: vec9
###################################################
u9 <- r3vel(9)
u9


###################################################
### code chunk number 70: vecfun
###################################################
f <- gyrfun(u9,v)
f(x)


###################################################
### code chunk number 71: u9+x
###################################################
u9+x


###################################################
### code chunk number 72: nonass
###################################################
(u+v)+w
u+(v+w)


###################################################
### code chunk number 73: nonass1
###################################################
(u+(v+w)) - ((u+v)+gyr(u,v,w))
((u+v)+w) - (u+(v+gyr(v,u,w)))


###################################################
### code chunk number 74: viss
###################################################
u <- as.3vel(c(0.4,0,0))
v <- seq(as.3vel(c(0.4,-0.2,0)), as.3vel(c(-0.3,0.9,0)),len=20)
w <- as.3vel(c(0.8,-0.4,0))


###################################################
### code chunk number 75: comfail1_fig
###################################################
comm_fail1(u=u, v=v)


###################################################
### code chunk number 76: comfail2_fig
###################################################
comm_fail2(u=u, v=v)


###################################################
### code chunk number 77: assfail_fig
###################################################
ass_fail(u=u, v=v, w=w, bold=10)


###################################################
### code chunk number 78: defuvw
###################################################
 u <- as.3vel(c(+0.5,0.1,-0.2))
 v <- as.3vel(c(+0.4,0.3,-0.2))
 w <- as.3vel(c(-0.3,0.2,+0.2))


###################################################
### code chunk number 79: lorentz.Rnw:1113-1116
###################################################
jj1 <- u %>% add(v)
jj2 <- u+v
speed(jj1-jj2)


###################################################
### code chunk number 80: lorentz.Rnw:1121-1124
###################################################
jj1 <- u %>% add(v) %>% add(w)
jj2 <- (u+v)+w
speed(jj1-jj2)


###################################################
### code chunk number 81: lorentz.Rnw:1131-1134
###################################################
jj1 <- u %>% add(v %>% add(w))
jj2 <- u+(v+w)
speed(jj1-jj2)


###################################################
### code chunk number 82: funcnotation
###################################################
x <- as.3vel(c(0.7, 0.0, -0.7))
y <- as.3vel(c(0.1, 0.3, -0.6))
u <- as.3vel(c(0.0, 0.8, +0.1))   # x,y,u: single three-velocities
v <- r3vel(5,0.9)
w <- r3vel(5,0.8)                 # v,w: vector of three-velocities
f <- gyrfun(u,v)
g <- gyrfun(v,u)


###################################################
### code chunk number 83: testeq3to8
###################################################
max(speed((u+v) - f(v+u)))              # equation 3
max(abs(prod3(f(x),f(y)) - prod3(x,y))) # equation 4
max(speed(f(x+y) - (f(x)+f(y))))        # equation 5
max(speed(f(g(x)) - g(f(x))))           # equation 6
max(speed((u+(v+w)) - ((u+v)+f(w))))    # equation 7
max(speed(((u+v)+w) - (u+(v+g(w)))))    # equation 8