824 lines
24 KiB
C++
824 lines
24 KiB
C++
#include "stdafx.h"
|
|
|
|
#ifdef MAC
|
|
#define EOL "\n"
|
|
#else
|
|
#define EOL "\r\n"
|
|
#endif
|
|
|
|
char *example_script[8] = {
|
|
|
|
// 0. Smiley face
|
|
|
|
"# Next, click the Run Script button"EOL
|
|
EOL
|
|
"f = quote(test("EOL
|
|
" t < 0, 5 (cos(2 t), sin(2 t)),"EOL
|
|
" t < pi / 4, (-2, 2),"EOL
|
|
" t < pi / 2, (2, 2),"EOL
|
|
" 3 (cos(2 t), sin(2 t))"EOL
|
|
"))"EOL
|
|
EOL
|
|
"draw(f)"EOL,
|
|
|
|
// 1. GammaMatrixAlgebra
|
|
|
|
"# This script does a few of the exercises from Feynman's book \"Quantum Electrodynamics.\""EOL
|
|
"# To run this script, click on the 'Run Script' button (in the main window)."EOL
|
|
""EOL
|
|
"# Define the spacetime metric (for multiplying spacetime vectors)."EOL
|
|
""EOL
|
|
"metric = ((-1, 0, 0, 0),"EOL
|
|
" ( 0,-1, 0, 0),"EOL
|
|
" ( 0, 0,-1, 0),"EOL
|
|
" ( 0, 0, 0, 1))"EOL
|
|
""EOL
|
|
"# Define I, the identity matrix."EOL
|
|
""EOL
|
|
"I = ((1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1))"EOL
|
|
""EOL
|
|
"# Define the gamma matrices."EOL
|
|
""EOL
|
|
"gammax = (( 0, 0, 0, 1),"EOL
|
|
" ( 0, 0, 1, 0),"EOL
|
|
" ( 0,-1, 0, 0),"EOL
|
|
" (-1, 0, 0, 0))"EOL
|
|
""EOL
|
|
"gammay = (( 0, 0, 0,-i),"EOL
|
|
" ( 0, 0, i, 0),"EOL
|
|
" ( 0, i, 0, 0),"EOL
|
|
" (-i, 0, 0, 0))"EOL
|
|
""EOL
|
|
"gammaz = (( 0, 0, 1, 0),"EOL
|
|
" ( 0, 0, 0,-1),"EOL
|
|
" (-1, 0, 0, 0),"EOL
|
|
" ( 0, 1, 0, 0))"EOL
|
|
""EOL
|
|
"gammat = (( 1, 0, 0, 0),"EOL
|
|
" ( 0, 1, 0, 0),"EOL
|
|
" ( 0, 0,-1, 0),"EOL
|
|
" ( 0, 0, 0,-1))"EOL
|
|
""EOL
|
|
"# Define the gamma vector."EOL
|
|
"#"EOL
|
|
"# The gamma vector has gamma matrices for its components. We express it here"EOL
|
|
"# as a rank 3 tensor. We set up the tensor so that the vector component index"EOL
|
|
"# is the last (rightmost) index. With this configuration we can left-multiply"EOL
|
|
"# with a Feynman slash matrix using the dot function."EOL
|
|
"#"EOL
|
|
"# For example, in component notation, this is how we want to multiply with a"EOL
|
|
"# Feynman slash matrix:"EOL
|
|
"#"EOL
|
|
"# aslash[a,b] gamma[b,c,d]"EOL
|
|
"#"EOL
|
|
"# (summation over the repeated index b)"EOL
|
|
"#"EOL
|
|
"# The summation over b is exactly what the dot function does so we can do the"EOL
|
|
"# above multiply with dot(aslash,gamma)."EOL
|
|
"#"EOL
|
|
"# In the following outer products, placing the basis vector operands on the"EOL
|
|
"# right-hand side results in the desired index ordering."EOL
|
|
""EOL
|
|
"gamma = outer(gammax,(1,0,0,0)) +"EOL
|
|
" outer(gammay,(0,1,0,0)) +"EOL
|
|
" outer(gammaz,(0,0,1,0)) +"EOL
|
|
" outer(gammat,(0,0,0,1))"EOL
|
|
""EOL
|
|
"# DOT is for multiplying gamma vectors. This is a special multiply because we"EOL
|
|
"# have to dot the individual vector components (the gamma matrices) then we"EOL
|
|
"# have to sum over all the results. In component notation, this is how we want"EOL
|
|
"# to do the multiply:"EOL
|
|
"#"EOL
|
|
"# T[a,c] = A[a,b,d] B[b,c,d]"EOL
|
|
"#"EOL
|
|
"# To do this, we start with an outer product which results in the following"EOL
|
|
"# rank 6 tensor:"EOL
|
|
"#"EOL
|
|
"# T[a,b,d,b,c,d]"EOL
|
|
"#"EOL
|
|
"# Next we sum over b (indices 2 and 4) to get the following:"EOL
|
|
"#"EOL
|
|
"# T[a,d,c,d]"EOL
|
|
"#"EOL
|
|
"# Then we sum over d (indices 2 and 4 again) to get"EOL
|
|
"#"EOL
|
|
"# T[a,c]"EOL
|
|
"#"EOL
|
|
"# One final note, dot(B,metric) applies the spacetime metric to the rightmost"EOL
|
|
"# index of B, the vector index."EOL
|
|
""EOL
|
|
"DOT(A,B) = contract(contract(outer(A,dot(B,metric)),2,4),2,4)"EOL
|
|
""EOL
|
|
"# Define arbitrary spacetime vectors a, b and c."EOL
|
|
""EOL
|
|
"a = (ax,ay,az,at)"EOL
|
|
"b = (bx,by,bz,bt)"EOL
|
|
"c = (cx,cy,cz,ct)"EOL
|
|
""EOL
|
|
"# Define generic Feynman slash matrices."EOL
|
|
"# Note: The order of dot operands here is different from the book. This is"EOL
|
|
"# because we defined gamma to have its vector index on the right. Therefore"EOL
|
|
"# we have to right-multiply with the spacetime vector so that dot contracts"EOL
|
|
"# over the correct indices. In component notation we have"EOL
|
|
"#"EOL
|
|
"# aslash[u,v] = gamma[u,v,w] a[w]"EOL
|
|
"#"EOL
|
|
"# where summation is over the repeated index w."EOL
|
|
""EOL
|
|
"aslash = dot(gamma,metric,a)"EOL
|
|
"bslash = dot(gamma,metric,b)"EOL
|
|
"cslash = dot(gamma,metric,c)"EOL
|
|
""EOL
|
|
"# The Feynman slash matrices are 4x4 matrices. For example, aslash looks like"EOL
|
|
"# this:"EOL
|
|
"#"EOL
|
|
"# at 0 -az -ax + i ay"EOL
|
|
"#"EOL
|
|
"# 0 at -ax - i ay az"EOL
|
|
"#"EOL
|
|
"# az ax - i ay -at 0"EOL
|
|
"#"EOL
|
|
"# ax + i ay -az 0 -at"EOL
|
|
""EOL
|
|
"# Now we are ready to try the exercises. We want to show that each of the"EOL
|
|
"# following identities is true."EOL
|
|
""EOL
|
|
"\"Checking the following identities:\""EOL
|
|
""EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
"#"EOL
|
|
"# aslash = at gammat - ax gammax - ay gammay - az gammaz"EOL
|
|
"#"EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
""EOL
|
|
"display(quote(aslash = at gammat - ax gammax - ay gammay - az gammaz))"EOL
|
|
""EOL
|
|
"A = aslash"EOL
|
|
""EOL
|
|
"B = at gammat - ax gammax - ay gammay - az gammaz"EOL
|
|
""EOL
|
|
"check(A = B) # if A = B then continue, else stop"EOL
|
|
""EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
"#"EOL
|
|
"# aslash bslash = -bslash aslash + 2 a b"EOL
|
|
"#"EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
""EOL
|
|
"display(quote(aslash bslash = -bslash aslash + 2 a b))"EOL
|
|
""EOL
|
|
"A = dot(aslash,bslash)"EOL
|
|
""EOL
|
|
"B = -dot(bslash,aslash) + 2 dot(a,metric,b) I"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
"#"EOL
|
|
"# gamma gamma = 4"EOL
|
|
"#"EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
""EOL
|
|
"display(quote(gamma gamma = 4))"EOL
|
|
""EOL
|
|
"A = DOT(gamma,gamma)"EOL
|
|
""EOL
|
|
"B = 4 I"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
"#"EOL
|
|
"# gamma aslash gamma = -2 aslash"EOL
|
|
"#"EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
""EOL
|
|
"display(quote(gamma aslash gamma = -2 aslash))"EOL
|
|
""EOL
|
|
"A = DOT(gamma,dot(aslash,gamma))"EOL
|
|
""EOL
|
|
"B = -2 aslash"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
"#"EOL
|
|
"# gamma aslash bslash gamma = 4 a b"EOL
|
|
"#"EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
""EOL
|
|
"display(quote(gamma aslash bslash gamma = 4 a b))"EOL
|
|
""EOL
|
|
"A = DOT(gamma,dot(aslash,bslash,gamma))"EOL
|
|
""EOL
|
|
"B = 4 dot(a,metric,b) I"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
"#"EOL
|
|
"# gamma aslash bslash cslash gamma = -2 cslash bslash aslash"EOL
|
|
"#"EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
""EOL
|
|
"display(quote(gamma aslash bslash cslash gamma = -2 cslash bslash aslash))"EOL
|
|
""EOL
|
|
"A = DOT(gamma,dot(aslash,bslash,cslash,gamma))"EOL
|
|
""EOL
|
|
"B = -2 dot(cslash,bslash,aslash)"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
"#"EOL
|
|
"# If we get here then everything worked."EOL
|
|
"#"EOL
|
|
"#------------------------------------------------------------------------------"EOL
|
|
""EOL
|
|
"\"OK\""EOL,
|
|
|
|
// 2. VectorCalculus
|
|
|
|
"# This script tests 10 vector calculus identities."EOL
|
|
""EOL
|
|
"# Define the cross product, div, grad, curl and laplacian for"EOL
|
|
"# rectangular coordinates."EOL
|
|
""EOL
|
|
"cross(u,v) = ("EOL
|
|
" u[2] v[3] - u[3] v[2],"EOL
|
|
" u[3] v[1] - u[1] v[3],"EOL
|
|
" u[1] v[2] - u[2] v[1]"EOL
|
|
")"EOL
|
|
""EOL
|
|
"div(v) = contract(d(v,(x,y,z)),1,2)"EOL
|
|
""EOL
|
|
"grad(v) = d(v,(x,y,z))"EOL
|
|
""EOL
|
|
"curl(f) = ("EOL
|
|
" d(f[3],y) - d(f[2],z),"EOL
|
|
" d(f[1],z) - d(f[3],x),"EOL
|
|
" d(f[2],x) - d(f[1],y)"EOL
|
|
")"EOL
|
|
""EOL
|
|
"laplacian(f) = d(d(f,x),x) + d(d(f,y),y) + d(d(f,z),z)"EOL
|
|
""EOL
|
|
"# Note: Functions can be left undefined, such as FX(), FY(), etc."EOL
|
|
"# These \"generic\" functions, when evaluated by the derivative function d(),"EOL
|
|
"# are considered to be dependent on all variables."EOL
|
|
"# Basically what this means is that d() does no evaluation at all."EOL
|
|
"# For example, d(FX(),x) returns the expression d(FX(),x)."EOL
|
|
""EOL
|
|
"# Define generic vector functions F and G."EOL
|
|
""EOL
|
|
"F = (FX(),FY(),FZ())"EOL
|
|
"G = (GX(),GY(),GZ())"EOL
|
|
""EOL
|
|
"# Now check the 10 identities."EOL
|
|
""EOL
|
|
"\"Checking the following identities:\""EOL
|
|
""EOL
|
|
"\"1. div(curl F) = 0\""EOL
|
|
""EOL
|
|
"A = div(curl(F))"EOL
|
|
""EOL
|
|
"check(A = 0)"EOL
|
|
""EOL
|
|
"\"2. curl(grad f) = 0\""EOL
|
|
""EOL
|
|
"A = curl(grad(f())) # Note the use of generic scalar function f() here."EOL
|
|
""EOL
|
|
"check(A = 0)"EOL
|
|
""EOL
|
|
"\"3. div(grad f) = laplacian f\""EOL
|
|
""EOL
|
|
"A = div(grad(f()))"EOL
|
|
""EOL
|
|
"B = laplacian(f())"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"\"4. curl(curl F) = grad(div F) - laplacian F\""EOL
|
|
""EOL
|
|
"A = curl(curl(F))"EOL
|
|
""EOL
|
|
"B = grad(div(F)) - laplacian(F)"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"\"5. grad(fg) = f grad g + g grad(f)\""EOL
|
|
""EOL
|
|
"A = grad(f() g())"EOL
|
|
""EOL
|
|
"B = f() grad(g()) + g() grad(f())"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"\"6. grad(F . G) = (G . grad)F + (F . grad)G + G x curl F + F x curl G\""EOL
|
|
""EOL
|
|
"A = grad(dot(F,G))"EOL
|
|
""EOL
|
|
"B = dot(grad(F),G)+dot(grad(G),F)+cross(G,curl(F))+cross(F,curl(G))"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"# Note: It turns out that (G . grad)F actually means (grad F) . G"EOL
|
|
""EOL
|
|
"\"7. div(fF) = f div F + grad f . F\""EOL
|
|
""EOL
|
|
"A = div(f() F)"EOL
|
|
""EOL
|
|
"B = f() div(F) + dot(grad(f()),F)"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"\"8. div(F x G) = G . curl F - F . curl G\""EOL
|
|
""EOL
|
|
"A = div(cross(F,G))"EOL
|
|
""EOL
|
|
"B = dot(G,curl(F)) - dot(F,curl(G))"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"\"9. curl(fF) = f curl F + grad f x F\""EOL
|
|
""EOL
|
|
"A = curl(f() F)"EOL
|
|
""EOL
|
|
"B = f() curl(F) + cross(grad(f()),F)"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"\"10. curl(F x G) = F div G - G div F + (G . grad)F - (F . grad)G\""EOL
|
|
""EOL
|
|
"A = curl(cross(F,G))"EOL
|
|
""EOL
|
|
"B = F div(G) - G div(F) + dot(grad(F),G) - dot(grad(G),F)"EOL
|
|
""EOL
|
|
"check(A = B)"EOL
|
|
""EOL
|
|
"# If we get here then everything worked."EOL
|
|
""EOL
|
|
"\"OK\""EOL,
|
|
|
|
// 3. RotationMatrix
|
|
|
|
"# This script demonstrates that the rotation matrix is orthogonal."EOL
|
|
"# We also demonstrate that the rotation matrix leaves the length of a vector"EOL
|
|
"# unchanged."EOL
|
|
"#"EOL
|
|
"# For an orthogonal matrix, inverse and transpose are equivalent."EOL
|
|
"#"EOL
|
|
"# -1 T"EOL
|
|
"# R = R"EOL
|
|
"#"EOL
|
|
"# First we define the rotation matrix R then we convert sine and cosine"EOL
|
|
"# functions to their exponential forms."EOL
|
|
"#"EOL
|
|
"# The exponential forms lead to the required trigonometric simplifications."EOL
|
|
""EOL
|
|
"# These are the components of the rotation matrix."EOL
|
|
""EOL
|
|
"R11 = cos(phi2) cos(phi1) - cos(theta) sin(phi1) sin(phi2)"EOL
|
|
"R12 = -cos(phi2) sin(phi1) - cos(theta) cos(phi1) sin(phi2)"EOL
|
|
"R13 = sin(phi2) sin(theta)"EOL
|
|
""EOL
|
|
"R21 = sin(phi2) cos(phi1) + cos(theta) sin(phi1) cos(phi2)"EOL
|
|
"R22 = -sin(phi2) sin(phi1) + cos(theta) cos(phi1) cos(phi2)"EOL
|
|
"R23 = -cos(phi2) sin(theta)"EOL
|
|
""EOL
|
|
"R31 = sin(theta) sin(phi1)"EOL
|
|
"R32 = sin(theta) cos(phi1)"EOL
|
|
"R33 = cos(theta)"EOL
|
|
""EOL
|
|
"# R is the rotation matrix."EOL
|
|
""EOL
|
|
"R = ((R11,R12,R13),(R21,R22,R23),(R31,R32,R33))"EOL
|
|
""EOL
|
|
"# Print R"EOL
|
|
""EOL
|
|
"R"EOL
|
|
""EOL
|
|
"# Convert R to exponential form."EOL
|
|
""EOL
|
|
"R = circexp(R)"EOL
|
|
""EOL
|
|
"# Subtract transpose from inverse. The difference should be zero."EOL
|
|
""EOL
|
|
"check(inv(R) = transpose(R,1,2))"EOL
|
|
""EOL
|
|
"# Show that the rotation matrix leaves the length of a vector unchanged."EOL
|
|
"# The length difference between U and RU should be zero."EOL
|
|
""EOL
|
|
"U = (U1,U2,U3)"EOL
|
|
""EOL
|
|
"check(U^2 = dot(R,U)^2)"EOL
|
|
""EOL
|
|
"\"OK\""EOL,
|
|
|
|
// 4. QuantumHarmonicOscillator
|
|
|
|
"# \"Harmonic oscillator\" is the generic term for a system that has potential"EOL
|
|
"# energy V proportional to x squared."EOL
|
|
"#"EOL
|
|
"# For total energy E, kinetic energy K and potential energy V we have"EOL
|
|
"#"EOL
|
|
"# E = K + V"EOL
|
|
"#"EOL
|
|
"# For QHO the equivalent Schroedinger equation is:"EOL
|
|
"#"EOL
|
|
"# (2n + 1) psi = -d^2 psi / dx^2 + x^2 psi"EOL
|
|
"#"EOL
|
|
"# This differential equation can only be solved for integer values of n."EOL
|
|
"# The fact that n must be an integer is in agreement with physical reality:"EOL
|
|
"# the total energy of a QHO system is quantized."EOL
|
|
"#"EOL
|
|
"# The solution to the above equation is"EOL
|
|
"#"EOL
|
|
"# psi(n,x) = exp(-1/2 x^2) Hn(x)"EOL
|
|
"#"EOL
|
|
"# Hn(x) is the nth Hermite polynomial in x."EOL
|
|
""EOL
|
|
"# Define the wave function."EOL
|
|
" "EOL
|
|
"psi(n) = exp(-1/2 x^2) hermite(x,n)"EOL
|
|
""EOL
|
|
"# We want to show that psi does indeed solve the Schroedinger equation."EOL
|
|
"# We start by defining the energy functions E, K and V."EOL
|
|
""EOL
|
|
"E(n) = (2 n + 1) psi(n)"EOL
|
|
""EOL
|
|
"K(n) = -d(d(psi(n),x),x)"EOL
|
|
""EOL
|
|
"V(n) = x^2 psi(n)"EOL
|
|
""EOL
|
|
"# Display a few wave functions."EOL
|
|
""EOL
|
|
"psi8 = condense(psi(8))"EOL
|
|
"psi9 = condense(psi(9))"EOL
|
|
"psi10 = condense(psi(10))"EOL
|
|
""EOL
|
|
"display(psi8)"EOL
|
|
"display(psi9)"EOL
|
|
"display(psi10)"EOL
|
|
""EOL
|
|
"\"Checking E = K + V for n = 0 to 10\""EOL
|
|
""EOL
|
|
"for(n, 0, 10, check(E(n) = K(n) + V(n)))"EOL
|
|
""EOL
|
|
"\"OK\""EOL,
|
|
|
|
// 5. HydrogenicWavefunctions
|
|
|
|
"# This script is a \"math manipulative\" for hydrogenic wavefunctions."EOL
|
|
"#"EOL
|
|
"# First we generate hydrogenic wavefunctions using all three quantum numbers."EOL
|
|
"#"EOL
|
|
"# Then we check that the wavefunctions solve the Schroedinger equation for the"EOL
|
|
"# hydrogen atom:"EOL
|
|
"#"EOL
|
|
"# psi(n,L,m)/n^2 = del^2 psi(n,L,m) + 2 psi(n,L,m)/r"EOL
|
|
"#"EOL
|
|
"# We use upper case L to avoid confusion with the numeral 1."EOL
|
|
""EOL
|
|
"# Define the Laplacian operator for spherical coordinates (del-squared)."EOL
|
|
""EOL
|
|
"delsq(x) = 1/r^2 d(r^2 d(x,r),r) +"EOL
|
|
" 1/(r^2 sin(theta)) d(sin(theta) d(x,theta),theta) +"EOL
|
|
" 1/(r sin(theta))^2 d(d(x,phi),phi)"EOL
|
|
""EOL
|
|
"# psi is the product of radial wavefunction R and spherical harmonic Y"EOL
|
|
""EOL
|
|
"R(n,L) = r^L exp(-r/n) laguerre(2r/n,n-L-1,2L+1)"EOL
|
|
""EOL
|
|
"Y(L,m) = legendre(cos(theta),L,abs(m)) exp(i m phi)"EOL
|
|
""EOL
|
|
"psi(n,L,m) = R(n,L) Y(L,m)"EOL
|
|
""EOL
|
|
"# Define the energy functions E, K and V (from the Schroedinger equation)."EOL
|
|
""EOL
|
|
"E(n,L,m) = psi(n,L,m) / n^2"EOL
|
|
""EOL
|
|
"K(n,L,m) = delsq(psi(n,L,m))"EOL
|
|
""EOL
|
|
"V(n,L,m) = 2 psi(n,L,m) / r"EOL
|
|
""EOL
|
|
"\"Checking E = K + V for n = 1..5, L = 0..n-1, m = -L..L\""EOL
|
|
""EOL
|
|
"F(n,L,m) = simplify(E(n,L,m) - K(n,L,m) - V(n,L,m))"EOL
|
|
""EOL
|
|
"for(n,1,5,"EOL
|
|
" for(L,0,n-1,"EOL
|
|
" for(m,-L,L,"EOL
|
|
" check(F(n,L,m) = 0)"EOL
|
|
")))"EOL
|
|
""EOL
|
|
"# If we get here then everything worked, print OK."EOL
|
|
""EOL
|
|
"\"OK\""EOL
|
|
""EOL
|
|
"# Print a few wavefunctions."EOL
|
|
""EOL
|
|
"psi320 = condense(psi(3,2,0))"EOL
|
|
"psi321 = condense(psi(3,2,1))"EOL
|
|
"psi322 = condense(psi(3,2,2))"EOL
|
|
""EOL
|
|
"display(psi320)"EOL
|
|
"display(psi321)"EOL
|
|
"display(psi322)"EOL,
|
|
|
|
// 6. StaticSphericalMetric
|
|
|
|
"# This script calculates the Einstein tensor for a static spherically symmetric"EOL
|
|
"# metric."EOL
|
|
"#"EOL
|
|
"# Cf. \"A first course in general relativity,\" Bernard F. Schutz, p. 255."EOL
|
|
"#"EOL
|
|
"# The book tells us exactly what the Einstein tensor components should be."EOL
|
|
"# If we get the right answer then we can be reasonably sure that the script is"EOL
|
|
"# correct. Once that is known then we can use the functions defined here in"EOL
|
|
"# other scripts."EOL
|
|
"#"EOL
|
|
"# This is the line element for the metric (Equation 10.7)"EOL
|
|
"#"EOL
|
|
"# 2 2 Phi 2 2 Lambda 2 2 2"EOL
|
|
"# ds = -e dt + e dr + r d Omega"EOL
|
|
"#"EOL
|
|
"# where"EOL
|
|
"#"EOL
|
|
"# 2 2 2 2 2 2"EOL
|
|
"# r d Omega = r (d theta + sin theta d phi )"EOL
|
|
"#"EOL
|
|
"# Note: Phi and Lambda are both functions of r."EOL
|
|
""EOL
|
|
"# Given the line element we can write the metric tensor by inspection:"EOL
|
|
""EOL
|
|
"gdd = ((-exp(2 Phi(r)), 0, 0, 0),"EOL
|
|
" ( 0, exp(2 Lambda(r)), 0, 0),"EOL
|
|
" ( 0, 0, r^2, 0),"EOL
|
|
" ( 0, 0, 0, r^2 sin(theta)^2))"EOL
|
|
""EOL
|
|
"# Note: \"dd\" stands for two \"down\" indices, \"uu\" stands for two \"up\" indices."EOL
|
|
""EOL
|
|
"# X is our coordinate system. We need this for computing gradients."EOL
|
|
""EOL
|
|
"X = (t,r,theta,phi)"EOL
|
|
""EOL
|
|
"# Step 1: Calculate guu."EOL
|
|
""EOL
|
|
"guu = inv(gdd)"EOL
|
|
""EOL
|
|
"# Step 2: Calculate the connection coefficients. Cf. Gravitation, p. 210."EOL
|
|
"#"EOL
|
|
"# Gamma = 1/2 (g + g - g )"EOL
|
|
"# abc ab,c ac,b bc,a"EOL
|
|
"#"EOL
|
|
"# Note: The comma means gradient which increases the rank of gdd by 1."EOL
|
|
""EOL
|
|
"gddd = d(gdd,X)"EOL
|
|
""EOL
|
|
"# Note: We transpose indices so they match up with Gamma, i.e., we put them in"EOL
|
|
"# alphabetical order."EOL
|
|
""EOL
|
|
"GAMDDD = 1/2 (gddd + # indices are already in correct order"EOL
|
|
"transpose(gddd,2,3) - # transpose c and b"EOL
|
|
"transpose(transpose(gddd,2,3),1,2)) # transpose c and a, then b and a"EOL
|
|
""EOL
|
|
"# Raise first index."EOL
|
|
"#"EOL
|
|
"# a au"EOL
|
|
"# Gamma = g Gamma"EOL
|
|
"# bc ubc"EOL
|
|
"#"EOL
|
|
"# Note: Sum over index u means contraction."EOL
|
|
""EOL
|
|
"GAMUDD = contract(outer(guu,GAMDDD),2,3)"EOL
|
|
""EOL
|
|
"# Step 3. Calculate the Riemann tensor. Cf. Gravitation, p. 219."EOL
|
|
"#"EOL
|
|
"# a is alpha"EOL
|
|
"# b is beta"EOL
|
|
"# c is gamma"EOL
|
|
"# d is delta"EOL
|
|
"# u is mu"EOL
|
|
"#"EOL
|
|
"# a a a a u a u"EOL
|
|
"# R = Gamma - Gamma + Gamma Gamma - Gamma Gamma"EOL
|
|
"# bcd bd,c bc,d uc bd ud bc"EOL
|
|
"#"EOL
|
|
"# Do the gradient once and save in a temporary variable."EOL
|
|
""EOL
|
|
"tmp1 = d(GAMUDD,X)"EOL
|
|
""EOL
|
|
"# The Gamma Gamma product is a rank 6 tensor with dim 4 per rank."EOL
|
|
"# That works out to 4 to the 6th or 4,096 elements."EOL
|
|
"# Of course, we'll do the outer product and contract over u just once and save"EOL
|
|
"# the result in a second temporary variable."EOL
|
|
""EOL
|
|
"tmp2 = contract(outer(GAMUDD,GAMUDD),2,4)"EOL
|
|
""EOL
|
|
"# Now put it all together. Do the transpositions so the indices get matched up"EOL
|
|
"# with R on the left, i.e., put them in alphabetical order."EOL
|
|
""EOL
|
|
"RUDDD = transpose(tmp1,3,4) - # transpose d and c"EOL
|
|
" tmp1 + # already in correct order"EOL
|
|
" transpose(tmp2,2,3) - # transpose c and b"EOL
|
|
" transpose(transpose(tmp2,2,3),3,4) # transpose d and b, then d and c"EOL
|
|
""EOL
|
|
"# Step 4: Calculate the Ricci tensor. Cf. Gravitation, p. 343."EOL
|
|
"#"EOL
|
|
"# a"EOL
|
|
"# R = R"EOL
|
|
"# uv uav"EOL
|
|
"#"EOL
|
|
"# Contract over \"a\" (1st and 3rd indices)."EOL
|
|
""EOL
|
|
"RDD = contract(RUDDD,1,3)"EOL
|
|
""EOL
|
|
"# Step 5: Calculate the Ricci scalar. Cf. Gravitation, p. 343."EOL
|
|
"#"EOL
|
|
"# uv"EOL
|
|
"# R = g R"EOL
|
|
"# vu ...the book has uv, does it give the same result?"EOL
|
|
"# Yes because the metric tensor is symmetric so it's ok to"EOL
|
|
"# transpose."EOL
|
|
"# I prefer vu because it looks like raising an index."EOL
|
|
""EOL
|
|
"R = contract(contract(outer(guu,RDD),2,3),1,2)"EOL
|
|
""EOL
|
|
"# Step 6: Finally, calculate the Einstein tensor. Cf. Gravitation, p. 343."EOL
|
|
"#"EOL
|
|
"# G = R - 1/2 g R"EOL
|
|
"# uv uv uv"EOL
|
|
""EOL
|
|
"GDD = RDD - 1/2 gdd R"EOL
|
|
""EOL
|
|
"# Next we compare this result with Schutz' book. Schutz p. 255 gives the"EOL
|
|
"# following Einstein tensor components (all other components are zero):"EOL
|
|
"#"EOL
|
|
"# 1 d"EOL
|
|
"# G = ---- exp(2 Phi) ---- [r (1 - exp(-2 Lambda))]"EOL
|
|
"# tt 2 dr"EOL
|
|
"# r"EOL
|
|
"#"EOL
|
|
"# 1 2"EOL
|
|
"# G = - ---- exp(2 Lambda) (1 - exp(-2 Lambda)) + --- Phi'"EOL
|
|
"# rr 2 r"EOL
|
|
"# r"EOL
|
|
"#"EOL
|
|
"# 2 2"EOL
|
|
"# G = r exp(-2 Lambda) [Phi'' + (Phi') + Phi'/r"EOL
|
|
"# theta theta"EOL
|
|
"#"EOL
|
|
"# - Phi' Lambda' - Lamda'/r]"EOL
|
|
"#"EOL
|
|
"# 2"EOL
|
|
"# G = sin theta G"EOL
|
|
"# phi phi theta theta"EOL
|
|
""EOL
|
|
"Gtt = 1/r^2 exp(2 Phi(r)) d(r (1 - exp(-2 Lambda(r))),r)"EOL
|
|
""EOL
|
|
"Grr = -1/r^2 exp(2 Lambda(r)) (1 - exp(-2 Lambda(r))) + 2/r d(Phi(r),r)"EOL
|
|
""EOL
|
|
"Gthetatheta = r^2 exp(-2 Lambda(r)) ("EOL
|
|
" d(d(Phi(r),r),r) +"EOL
|
|
" d(Phi(r),r)^2 +"EOL
|
|
" d(Phi(r),r) / r -"EOL
|
|
" d(Phi(r),r) d(Lambda(r),r) -"EOL
|
|
" d(Lambda(r),r) / r)"EOL
|
|
""EOL
|
|
"Gphiphi = sin(theta)^2 Gthetatheta"EOL
|
|
""EOL
|
|
"# Put together the expected tensor:"EOL
|
|
""EOL
|
|
"expect = ((Gtt, 0, 0, 0),"EOL
|
|
" ( 0, Grr, 0, 0),"EOL
|
|
" ( 0, 0, Gthetatheta, 0),"EOL
|
|
" ( 0, 0, 0, Gphiphi))"EOL
|
|
""EOL
|
|
"# Check that GDD is correct."EOL
|
|
""EOL
|
|
"check(GDD = expect)"EOL
|
|
""EOL
|
|
"# Display the non-zero components of GDD."EOL
|
|
""EOL
|
|
"display(Gtt)"EOL
|
|
"display(Grr)"EOL
|
|
"display(Gthetatheta)"EOL
|
|
"display(Gphiphi)"EOL
|
|
""EOL
|
|
"\"OK\""EOL,
|
|
|
|
// 7. FreeParticleDiracEquation
|
|
|
|
"# This script demonstrates the free particle Dirac equation and a few of its"EOL
|
|
"# solutions."EOL
|
|
"#"EOL
|
|
"# \"Free particle\" means that there is no force pushing the particle around. In"EOL
|
|
"# other words, the potential energy field V is zero everywhere. What remains is"EOL
|
|
"# just the total energy T and kinetic energy K."EOL
|
|
"#"EOL
|
|
"# T = K"EOL
|
|
"#"EOL
|
|
"# The equivalent Dirac equation is"EOL
|
|
"#"EOL
|
|
"# m Psi = i delslash Psi"EOL
|
|
"#"EOL
|
|
"# where m is the particle's rest mass, Psi is a 4-vector and delslash is a"EOL
|
|
"# differential operator involving gamma matrices."EOL
|
|
"#"EOL
|
|
"# In terms of physical units, the above equation has m in the wrong place. The"EOL
|
|
"# mass m should be a divisor in the kinetic energy term. However, the equation"EOL
|
|
"# is more cool-looking as shown and anyway that's how Feynman wrote it. We just"EOL
|
|
"# have to remember to divide the above equation by m to get physical units of"EOL
|
|
"# energy."EOL
|
|
""EOL
|
|
"# Verify a few solutions (Psi) to the above free-particle Dirac equation."EOL
|
|
""EOL
|
|
"# Define the spacetime metric."EOL
|
|
""EOL
|
|
"metric = ((-1,0,0,0),(0,-1,0,0),(0,0,-1,0),(0,0,0,1))"EOL
|
|
""EOL
|
|
"# Define the gamma matrices."EOL
|
|
""EOL
|
|
"gammax = (( 0, 0, 0, 1),( 0, 0, 1, 0),( 0,-1, 0, 0),(-1, 0, 0, 0))"EOL
|
|
"gammay = (( 0, 0, 0,-i),( 0, 0, i, 0),( 0, i, 0, 0),(-i, 0, 0, 0))"EOL
|
|
"gammaz = (( 0, 0, 1, 0),( 0, 0, 0,-1),(-1, 0, 0, 0),( 0, 1, 0, 0))"EOL
|
|
"gammat = (( 1, 0, 0, 0),( 0, 1, 0, 0),( 0, 0,-1, 0),( 0, 0, 0,-1))"EOL
|
|
""EOL
|
|
"# Define the delslash operator."EOL
|
|
""EOL
|
|
"delslash(f) ="EOL
|
|
" dot(gammax,d(f,x))+"EOL
|
|
" dot(gammay,d(f,y))+"EOL
|
|
" dot(gammaz,d(f,z))+"EOL
|
|
" dot(gammat,d(f,t))"EOL
|
|
""EOL
|
|
"# Define energy E, momentum p and coordinate X."EOL
|
|
""EOL
|
|
"E = sqrt(px^2 + py^2 + pz^2 + m^2)"EOL
|
|
"p = (px,py,pz,E)"EOL
|
|
"X = (x,y,z,t)"EOL
|
|
""EOL
|
|
"# Verify that p.p = m^2"EOL
|
|
""EOL
|
|
"check(dot(p,metric,p) = m^2) # continue if equal, else stop"EOL
|
|
""EOL
|
|
"# Define the solutions."EOL
|
|
""EOL
|
|
"PsiA = (E+m,0,pz,px+i py) exp(-i dot(p,metric,X))"EOL
|
|
"PsiB = (0,E+m,px-i py,-pz) exp(-i dot(p,metric,X))"EOL
|
|
"PsiC = (pz,px+i py,E+m,0) exp(i dot(p,metric,X))"EOL
|
|
"PsiD = (px-i py,-pz,0,E+m) exp(i dot(p,metric,X))"EOL
|
|
""EOL
|
|
"# Verify the solutions."EOL
|
|
""EOL
|
|
"check(m PsiA = i delslash(PsiA))"EOL
|
|
"check(m PsiB = i delslash(PsiB))"EOL
|
|
"check(m PsiC = i delslash(PsiC))"EOL
|
|
"check(m PsiD = i delslash(PsiD))"EOL
|
|
""EOL
|
|
"# Try a linear combination of all solutions."EOL
|
|
""EOL
|
|
"Psi = A PsiA + B PsiB + C PsiC + D PsiD"EOL
|
|
""EOL
|
|
"check(m Psi = i delslash(Psi))"EOL
|
|
""EOL
|
|
"# For PsiA and PsiB it turns out that"EOL
|
|
"#"EOL
|
|
"# i delslash Psi = pslash Psi"EOL
|
|
"#"EOL
|
|
"# So another form of the free particle Dirac equation is"EOL
|
|
"#"EOL
|
|
"# m Psi = pslash Psi"EOL
|
|
"#"EOL
|
|
"# Verify solutions to the above equation."EOL
|
|
""EOL
|
|
"# Define the gamma tensor. See Gamma Matrix Algebra for details."EOL
|
|
""EOL
|
|
"gamma ="EOL
|
|
" outer(gammax,(1,0,0,0)) +"EOL
|
|
" outer(gammay,(0,1,0,0)) +"EOL
|
|
" outer(gammaz,(0,0,1,0)) +"EOL
|
|
" outer(gammat,(0,0,0,1))"EOL
|
|
""EOL
|
|
"# Dot gamma with p to get pslash."EOL
|
|
""EOL
|
|
"pslash = dot(gamma,metric,p)"EOL
|
|
""EOL
|
|
"# Verify the solutions again."EOL
|
|
""EOL
|
|
"check(m PsiA = dot(pslash,PsiA))"EOL
|
|
"check(m PsiB = dot(pslash,PsiB))"EOL
|
|
"check(m PsiC = -dot(pslash,PsiC))"EOL
|
|
"check(m PsiD = -dot(pslash,PsiD))"EOL
|
|
""EOL
|
|
"# Display pslash on the computer screen."EOL
|
|
""EOL
|
|
"pslash = subst(quote(E),E,pslash) # subst letter E for the energy expression"EOL
|
|
""EOL
|
|
"display(pslash)"EOL
|
|
""EOL
|
|
"# If we get here then everything worked, print OK."EOL
|
|
""EOL
|
|
"\"OK\""EOL
|
|
};
|