Special Relativity

Material for JUAS 2021 lectures

John Jowett

DateString[]

Setup

Constants and units

Recent versions of Mathematica have a system for handling physical units which is based on string names for physical quantities.   The built-in documentation describes this.

Earlier versions used a package in which units were symbols such as Meter.   
As a long-time user, I built a package on top of that system which I still use and find easier and more convenient.   The names of the symbols for units are systematic and predictable (but not the same as conventional abbreviations which are not logical enough).

If the package can be found on your $Path, then it will be loaded with

Needs["Accelerator`ConstantsUnits`"]

To make this notebook self-contained, there is  a copy of the code in the following subsection

Copy of Accelerator`ConstantsUnits package

To make this notebook self contained (and not rely on the package being installed for the user), this will be evaluated on initialisation.  No need to look at it.

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)

Notation`AutoLoadNotationPalette=False
BeginPackage["Accelerator`ConstantsUnits`",{"Notation`","Units`","PhysicalConstants`"}]
If[$VersionNumber>8,Print["This package should still work in Mathematica Version 9.\nIt implements a system for managing physical constants and units that is different and independent of the new scheme implemented in Version 9.\nIt may be adapted for compatibility at a later date."]]
ConstantsUnits::usage ="ConstantsUnits.m is a package that provides notations and values for constants occurring in accelerator physics.";
PeV::usage=TeV::usage=GeV::usage=MeV::usage= keV::usage="keV, MeV, GeV, TeV, PeV are units of energy and evaluate to (Kilo ElectronVolt), ... etc.";
lumiUnit::usage="lumiUnit as the conventional unit of colliding beam luminosity, \!\(\*SuperscriptBox[\(cm\), \(-2\)]\)\!\(\*SuperscriptBox[\(s\), \(-1\)]\).";
c::usage="c is a symbol for the speed of light in vacuum; it evaluates numerically with N[ ].";
ℏ::usage="ℏ is a symbol for the reduced Planck constant; it evaluates numerically with N[ ].";
e::usage="e is a symbol for the electron charge; it evaluates numerically with N[ ].";
Notation[(Subscript[r, e]) ⟺ rElectron];
rElectron::usage="\!\(\*SubscriptBox[\(r\), \(e\)]\) or rElectron is a symbol for the classical radius of the electron; it evaluates numerically with N[ ].";
Notation[ParsedBoxWrapper[SubscriptBox["α", "s"]]⟺ParsedBoxWrapper["alphas"]]
alphas::usage="\!\(\*SubscriptBox[\(α\), \(s\)]\) or alphas is a symbol for the fine structure constant; it evaluates numerically with N[ ].";
Notation[(Subscript[λ, e]) ⟺ lambdaElectron];
lambdaElectron::usage="\!\(\*SubscriptBox[\(λ\), \(e\)]\) or lambdaElectron is a symbol for the Compton wavelength of the electron; it evaluates numerically with N[ ].";
Notation[(Subscript[G, e]) ⟺ GElectron];
GElectron::usage="\!\(\*SubscriptBox[\(G\), \(e\)]\) or GElectron is a symbol for the anomalous magnetic moment of the electron; it evaluates numerically with N[ ].";
Notation[(Subscript[μ, 0]) ⟺ mu0];
mu0::usage="\!\(\*SubscriptBox[\(μ\), \(0\)]\) or mu0 is a symbol for the permeability of the vacuum; it evaluates numerically with N[ ].";
Notation[(Subscript[ϵ, 0]) ⟺ epsilon0];
epsilon0::usage="epsilon0 is a symbol for the permittivity of the vacuum; it evaluates numerically with N[ ].";
Notation[(Subscript[m, e]) ⟺ mElectron];
mElectron::usage="\!\(\*SubscriptBox[\(m\), \(e\)]\) or mElectron is a symbol for the mass of the electron; it evaluates numerically with N[ ].";
Notation[(Subscript[m, p]) ⟺ mProton];
mProton::usage="\!\(\*SubscriptBox[\(m\), \(p\)]\) or mProton is a symbol for the mass of the proton; it evaluates numerically with N[ ].";
Notation[(Subscript[m, n]) ⟺ mNeutron];
mNeutron::usage="\!\(\*SubscriptBox[\(m\), \(n\)]\) or mNeutron is a symbol for the mass of the neutron; it evaluates numerically with N[ ].";
Notation[(Subscript[m, d]) ⟺ mDeuteron];
mDeuteron::usage="\!\(\*SubscriptBox[\(m\), \(d\)]\) or mDeuteron is a symbol for the mass of the deuteron; it evaluates numerically with N[ ].";
Notation[(Subscript[m, υ]) ⟺ mMuon];
mMuon::usage="\!\(\*SubscriptBox[\(m\), \(υ\)]\) or mMuon is a symbol for the mass of the muon; it evaluates numerically with N[ ].";
Notation[(Subscript[m, Z]) ⟺ mZ];
mZ::usage="\!\(\*SubscriptBox[\(m\), \(Z\)]\) or mZ is a symbol for the mass of the Z-particle; it evaluates numerically with N[ ].";
Notation[(Subscript[m, W]) ⟺ mW];
mW::usage="\!\(\*SubscriptBox[\(m\), \(\(W\)\(\\\ \)\)]\)or mW is a symbol for the mass of the W-particle; it evaluates numerically with N[ ].";
Notation[(Subscript[a, ∞]) ⟺ aBohr];
aBohr::usage="\!\(\*SubscriptBox[\(a\), \(∞\)]\) or aBohr is a symbol for the Bohr radius; it evaluates numerically with N[ ].";
Notation[(Subscript[r, p]) ⟺ rProton];
rProton::usage="\!\(\*SubscriptBox[\(r\), \(p\)]\) or rProton is a symbol for the classical radius of the proton; it evaluates numerically with N[ ].";
Notation[(Subscript[Z, 0]) ⟺ Z0];
Z0::usage="\!\(\*SubscriptBox[\(Z\), \(0\)]\) or Z0 is a symbol for the impedance of free space; it evaluates numerically with N[ ].";
Notation[(Subscript[σ, T]) ⟺ sigmaThomson];
sigmaThomson::usage="\!\(\*SubscriptBox[\(σ\), \(T\)]\) or sigmaThomson is a symbol for the Thomson cross section; it evaluates
numerically with N[ ].";
Notation[(Subscript[σ, SB]) ⟺ sigmaStefanBoltzmann];
sigmaStefanBoltzmann::usage="\!\(\*SubscriptBox[\(σ\), \(SB\)]\) or sigmaStefanBoltzmann is a symbol for the Stefan-Boltzmann constant; it evaluates numerically with N[ ].";
Notation[(Subscript[μ, B]) ⟺ muBohr];
muBohr::usage="\!\(\*SubscriptBox[\(μ\), \(B\)]\) or muBohr is the Bohr magneton; it evaluates numerically with N[ ].";
Notation[(Subscript[μ, N]) ⟺ muNuclear]
muNuclear::usage="\!\(\*SubscriptBox[\(μ\), \(N\)]\) or muNuclear is the nuclear magneton; it evaluates numerically with N[ ].";
ToFundamentalSI::usage="toFundamentalSI[f] converts the quantity f to fundamental SI units.  Generally this helps to simplify the expressions of the units.";
ClassicalProtonRadius::usage="ClassicalProtonRadius is the value of the classical radius of the proton.";
DimensionCheck::usage="DimensionCheck[f] returns an expression showing the dimensional consistency of f, based on the internal database of symbols and units maintained by the ConstantsUnits package.\nThe alternative form DimensionCheck[f,symbolsUnits] uses a list of {symbol,unit} pairs given as second argument instead of the internal database.";
NumericalValueUnits::usage="NumericalValueUnits[f] interprets f as the product of a numerical value fn and symbolic units fu and returns {fn,fu}.";
PhysicalUnits::usage="PhysicalUnits[f] interprets f as the product of a numerical value fn and symbolic units fu and returns fu.";
NumericalValue::usage="NumericalValue[f] interprets f as the product of a numerical value fn and symbolic units fu and returns fn.";
IntroduceSymbol::usage="IntroduceSymbol[symbol,usagetext,units] introduces a symbol, registers its units for purposes of dimensional analysis and defines a usage message for it.";
SymbolsUnits::usage="SymbolsUnits[] returns a list of all symbols and units registered in the internal database.";
PhysicalUnitsPlot::usage="PhysicalUnitsPlot[f,{x,xmin,xmax}] works like Plot but removes the physical units of f and x before plotting.  Further, it creates labels for the axes from the symbolic forms of f and x and indicates their units.  It has all options of Plot except AxesLabel.";
nuclearMassGeV::usage="nuclearMassGeV[{Z,A}] returns the mass of a nucleus in GeV."
ConstantsUnits::badarg = "You called `1` with argument `2`!"
Begin["`Private`"]
minGodelPrimeIndex = 10;
ClassicalProtonRadius=Convert[ElectronCharge^2/((4 π VacuumPermittivity) ProtonMass SpeedOfLight^2),Meter]
ToFundamentalSI[f_]:=PowerExpand[SI[f]//.Units`Private`$ToFundamental,{Kilogram,Meter,Second,Candela}]
DimensionCheck[f_,symbolsWithUnits_List]:=Expand[f//.MapIndexed[(#1[[1]]->Prime[minGodelPrimeIndex+First[#2]] #1[[2]]&),Union[symbolsWithUnits]]];
DimensionCheck[f_]:=DimensionCheck[f,Accelerator`ConstantsUnits`Private`$SymbolsUnits]
PeV=Peta ElectronVolt;TeV= Tera ElectronVolt;GeV=Giga ElectronVolt;
MeV=Mega ElectronVolt;keV=Kilo ElectronVolt;eV=ElectronVolt;
lumiUnit = 1/((Centi*Meter)^2*Second);
N[c]=SpeedOfLight;
N[e]=ElectronCharge;
N[ℏ]=PlanckConstantReduced;
N[mElectron]=ElectronMass;
N[mProton]=ProtonMass;N[mNeutron]=NeutronMass;N[mDeuteron]=DeuteronMass;
N[mMuon]=MuonMass;
N[alphas]=FineStructureConstant;
N[rElectron]=ClassicalElectronRadius;
N[rProton]=ClassicalProtonRadius;
N[lambdaElectron]=ElectronComptonWavelength;
N[lambdabarElectron]=N[ElectronComptonWavelength/(2 π)];N[aBohr]=BohrRadius;N[sigmaThomson]=ThomsonCrossSection;N[sigmaStefanBoltzmann]=StefanConstant;
N[GElectron]=If[$VersionNumber>=4.1,(-ElectronGFactor-2)/2,ElectronGFactor-1]
N[mu0]=Convert[VacuumPermeability,Henry/Meter]
N[epsilon0]=Convert[VacuumPermittivity,Farad/Meter]
N[mZ]=Convert[N[91187.6`5.999999999999999  (GeV/1000)/(SpeedOfLight^2)],Kilogram]
N[mW]=Convert[N[80403.`5.000000000000002  (GeV/1000)/(SpeedOfLight^2)],Kilogram]
N[muBohr]=Convert[N[(ElectronCharge PlanckConstantReduced)/(2 ElectronMass)],Joule/Tesla]
N[muNuclear]=Convert[N[(ElectronCharge PlanckConstantReduced)/(2 ProtonMass)],Joule/Tesla]
N[Z0]=Convert[N[VacuumPermeability SpeedOfLight],Ohm]
nuclearMassGeV[{z_,a_}]:=Convert[(IsotopeData[{z,a},"AtomicMass"]/Quantity[1,"AtomicMassUnit"] AMU-z ElectronMass )SpeedOfLight^2//N,GeV]
NumericalValueUnits[Times[a_/;NumericQ[a],b__]]:={a,Times@b};
NumericalValueUnits[a_/;NumericQ[a]]:={a,1};
NumericalValueUnits[a_]:={1,a}
NumericalValue[f_]:=First[NumericalValueUnits[f]]
PhysicalUnits[f_]:=Last[NumericalValueUnits[f]]
SetAttributes[PhysicalUnitsPlot,HoldAll];
Options[PhysicalUnitsPlot]=Options[Plot];
PhysicalUnitsPlot[f_,{var_,varmin_,varmax_},opts___?OptionQ]/;PhysicalUnits[ToFundamentalSI[varmin]]==PhysicalUnits[ToFundamentalSI[varmax]]:=
	Module[{vunits,funits,nvar},
		vunits=PhysicalUnits[N[varmin]];funits=PhysicalUnits[N[f/.var->varmin]];
		fnew=f/.var->nvar vunits;
		If[FreeQ[funits,Plus],Plot[Evaluate[fnew/funits],{nvar,ToFundamentalSI[varmin/vunits],ToFundamentalSI[varmax/vunits]},AxesLabel->{StandardForm[HoldForm[var] [vunits]],StandardForm[HoldForm[f][ funits]]},opts],Print["Units of ordinate appear to be inconsistent."]]]
$SymbolsUnits=({#1,PhysicalUnits[N[#1]]}&)/@ToExpression[Names["Accelerator`ConstantsUnits`*"]]
$SymbolsUnits=Select[$SymbolsUnits,First[#1]=!=Last[#1]&& !(NumericQ[First[#1]/Last[#1]])&]
IntroduceSymbol[sym_Symbol,usagetext_String,unit_]:=($SymbolsUnits=Append[Select[$SymbolsUnits,(First[#1]=!=sym)&],
{sym,unit}];
sym::"usage"=ToString[sym]<>" "<>usagetext<>" (units: "<>ToString[unit]<>")")
SymbolsUnits[]:=$SymbolsUnits
(* Sin/: Sin[x_]^2 := 1 - Cos[x]^2 *)
(* Protect[ Evaluate[protected] ] *)
End[ ]
Protect[Evaluate[$Context <> "*"]]
EndPackage[ ]
If[$Notebooks,NotebookPut[Notebook[{Cell[BoxData[{GridBox[{{ButtonBox[SubscriptBox["r","e"]],ButtonBox["c"],ButtonBox["e"],ButtonBox["ℏ"],ButtonBox[SubscriptBox["α","s"]],ButtonBox[SubscriptBox["λ","e"]],"\[Placeholder]",ButtonBox[SubscriptBox["Z","0"]]},{ButtonBox[SubscriptBox["r","p"]],ButtonBox[SubscriptBox["μ","0"]],ButtonBox[SubscriptBox["ϵ","0"]],ButtonBox[SubscriptBox["G","e"]],ButtonBox[SubscriptBox["a","∞"]],ButtonBox[SubscriptBox["μ","B"]],ButtonBox[SubscriptBox["μ","N"]],ButtonBox[SubscriptBox["σ","T"]]},{ButtonBox[SubscriptBox["m","e"]],ButtonBox[SubscriptBox["m","p"]],ButtonBox[SubscriptBox["m","n"]],ButtonBox[SubscriptBox["m","d"]],ButtonBox[SubscriptBox["m","μ"]],ButtonBox[SubscriptBox["m","Z"]],ButtonBox[SubscriptBox["m","W"]],ButtonBox[SubscriptBox["σ","SB"]]}},GridBoxSpacings->{"Columns"->{Offset[0.27999999999999997`],{Offset[0.`]},Offset[0.27999999999999997`]},"ColumnsIndexed"->{},"Rows"->{Offset[0.2`],{Offset[0.`]},Offset[0.2`]},"RowsIndexed"->{}},GridDefaultElement:>ButtonBox["\\[Placeholder]"]],"\n",GridBox[{{GridBox[{{ButtonBox[RowBox[{"Convert","[",RowBox[{"\[SelectionPlaceholder]",",","\[Placeholder]"}],"]"}]]," ",ButtonBox[RowBox[{"SI","[","\[SelectionPlaceholder]","]"}]]}}]},{ButtonBox[RowBox[{"ToFundamentalSI","[","\[SelectionPlaceholder]","]"}]]},{ButtonBox[RowBox[{"DimensionCheck","[","\[SelectionPlaceholder]","]"}]]},{ButtonBox[RowBox[{RowBox[{RowBox[{"Notation","[",RowBox[{TagBox["\[SelectionPlaceholder]",NotationBoxTag,Editable->True]," ","⟺"," ",StyleBox[TagBox["\[Placeholder]",NotationBoxTag,Editable->True],FontColor->]}],"]"}],";"}],"\n",RowBox[{RowBox[{"IntroduceSymbol","[",RowBox[{StyleBox["\[Placeholder]",FontColor->],",","\"is the \[Placeholder].\"",",","\[Placeholder]"}],"]"}],";"}]}]]},{ButtonBox[RowBox[{"SymbolsUnits","[","]"}]]}},GridBoxSpacings->{"Columns"->{Offset[0.27999999999999997`],{Offset[0.`]},Offset[0.27999999999999997`]},"ColumnsIndexed"->{},"Rows"->{Offset[0.2`],{Offset[0.`]},Offset[0.2`]},"RowsIndexed"->{}},GridDefaultElement:>ButtonBox["\\[Placeholder]"]]}],NotebookDefault,CellFrame->{{1,1},{0,1}},CellBracketOptions->{"Color"->},Evaluatable->True,CellGroupingRules->"InputGrouping",TextClipboardType->"InputText",PageBreakAbove->True,PageBreakWithin->False,GroupPageBreakWithin->False,CellLabelMargins->{{12,Inherited},{Inherited,Inherited}},DefaultFormatType->DefaultInputFormatType,InputAutoReplacements->{"->"->"->",":>"->":>","<="->"≤",">="->"≥","!="->"≠","=="->"==",ParentList},CellChangeTimes->{{3.4034414238841486`*^9,3.403441425962247`*^9},3.403441473149143`*^9,{3.4035298714365263`*^9,3.4035299105454006`*^9},{3.415450682999962`*^9,3.415450730545924`*^9},{3.4154512436298227`*^9,3.4154512564733257`*^9}},ShowAutoStyles->True,"TwoByteSyntaxCharacterAutoReplacement"->True,HyphenationOptions->{"HyphenationCharacter"->""},LineSpacing->{1.25`,0,2},AutoItalicWords->{},LanguageCategory->"Mathematica",ScriptMinSize->9,NumberMarks->True,SingleLetterItalics->False,LinebreakAdjustments->{0.85`,2,10,0,1},CounterIncrements->"Input",MenuCommandKey->"9",StyleMenuListing->None,FontFamily->"Courier",FontWeight->"Bold",Background->,CellTags->{"b:0.1","b:0.3.1","b:0.3.2","b:0.4.2","b:0.6.0","b:0.6.1","b:0.6.3","b:0.6.5","b:0.6.6","b:0.6.7","b:0.6.8","b:0.6.9","b:0.6.10","b:0.6.11","b:0.7"},$CellContext`StripStyleOnPaste->True,$CellContext`MenuPosition->1500]},ClosingAutoSave->True,Editable->False,WindowToolbars->{},WindowSize->{Fit,Fit},WindowMargins->{{Automatic,47},{Automatic,218}},WindowFrame->"Palette",WindowElements->{},WindowFrameElements->{"CloseBox","MinimizeBox"},WindowClickSelect->False,ScrollingOptions->{"PagewiseScrolling"->True,"VerticalScrollRange"->Fit},PrivateNotebookOptions->{"VersionedStylesheet"->{"Default.nb"[8.`]->False}},ShowCellBracket->False,CellMargins->0,Deployed->True,CellOpen->True,ShowCellLabel->False,ShowCellTags->False,ImageMargins->0,Magnification->1.`,FrontEndVersion->"9.0 for Microsoft Windows (64-bit) (November 20, 2012)",StyleDefinitions->"Default.nb"]]]

Needs["Accelerator`ElectronRings`"]

Some examples of handling units in the Accelerator`ConstantsUnits` package

The package generates a handy input palette that looks like this (does not work in cloud evaluation, just close it):

In this package, units are represented by capitalised symbols, in American spelling

?Meter

?Fortnight

?ToFundamentalSI

ToFundamentalSI[5.5 Centi Meter]

?GeV

?Convert

The function Convert changes the units in which a quantity is expressed

Convert[1 GeV, Joule]

?lumiUnit

Luminosity of a collider in conventional units multiplied by running time, converted to conventional units of integrated luminosity

1.7 (10^32)lumiUnit

Convert[(1.7 (10^32)lumiUnit)( 1.5 Month), (((Pico Barn))^(-1))]

SI[1.7 (10^32)lumiUnit 1.5 Month]

Some symbols represent physical constants and do not have values

(((e^2)/(4π (Subscript[ϵ, 0])))/( ℏ c))

but evaluate numerically with the N function (like π) when you want them to:

N[(((e^2)/(4π (Subscript[ϵ, 0])))/( ℏ c))]

To shake down the units use

?ToFundamentalSI

ToFundamentalSI[N[(((e^2)/(4π (Subscript[ϵ, 0])))/( ℏ c))]]

In this case the result is dimensionless, and the inverse of the fine-structure constant

1/%

In some US labs they used to sometimes use units such that c=1, or nearly:

Convert[N[c],Foot/(Nano Second)]

Other settings

SetDirectory[NotebookDirectory[]];
SetOptions[#,Axes->True,ImageSize->Medium,BaseStyle->{FontFamily->"Tahoma",FontSize->16}]&/@{Plot,LogLinearPlot,LogPlot,RegionPlot,ContourPlot,Plot3D,ContourPlot3D};
SetOptions[#,AxesStyle->Thick]&/@{Plot,LogLinearPlot,LogPlot,ContourPlot,RegionPlot};
SetOptions[#,FrameLabel->{"\!\(\*StyleBox[\"x\",FontSlant->\"Italic\"]\)","\!\(\*StyleBox[\"ct\",FontSlant->\"Italic\"]\)"},
AxesLabel->None]&/@{Plot,RegionPlot,ContourPlot};
SetOptions[#,AxesLabel->{"\!\(\*StyleBox[\"x\",FontSlant->\"Italic\"]\)","\!\(\*StyleBox[\"y\",FontSlant->\"Italic\"]\)","\!\(\*StyleBox[\"ct\",FontSlant->\"Italic\"]\)"}]&/@{Plot3D,ContourPlot3D};

Lorentz transformation

Reference frames

triad[p_String]:={Arrow[{{0,0,0},{1,0,0}}],Arrow[{{0,0,0},{0,1,0}}],Arrow[{{0,0,0},{0,0,1}}],
 Text["x"<>p,{1.1,0,0}],Text["y"<>p,{0,1.1,0}],Text["z"<>p,{0,0,1.1}]}

Graphics3D[{
{Black,Text["S",{.1,0,.8}],triad[""]},
Translate[{Red,Text["S'",{.1,0,.7}],
{Darker[Red],Arrow[{{.2,0,.5},{.4,0,.5}}],Text["v",{.44,0,.5}]},
triad["'"]},{.3,.041 ,0}]
},AxesLabel->{"x","y","z"}, Boxed->False,BaseStyle->Large,
ViewPoint->{1,-10,4}]

Manipulate[Graphics3D[{Text["t="<>ToString[t],{1,0,.9}],
{Black,Text["S",{.1,0,.8}],triad[""]},
Translate[{Red,Text["S'",{.1,0,.7}],
{Darker[Red],Arrow[{{.2,0,.5},{.4,0,.5}}],Text["v",{.44,0,.5}]},
triad["'"]},{t,.041 ,0}]
},AxesLabel->{"x","y","z"}, PlotRange->{{-1,2},{0,1},{0,1}},
Boxed->False,BaseStyle->Large,ImageSize->Large,
ViewPoint->{1,-10,4}],{{t,0},-1,1}]

Video for slides

Export["frames.mp4",Table[Graphics3D[{Text["t="<>ToString[t],{1,0,.9}],
{Black,Text["S",{.1,0,.8}],triad[""]},
Translate[{Red,Text["S'",{.1,0,.7}],
{Darker[Red],Arrow[{{.2,0,.5},{.4,0,.5}}],Text["v",{.44,0,.5}]},
triad["'"]},{t,.041 ,0}]
},AxesLabel->{"x","y","z"}, PlotRange->{{-1,2},{0,1},{0,1}},
Boxed->False,BaseStyle->Large,ImageSize->Large,
ViewPoint->{1,-10,4}],{t,-1,1,.05}]]

SystemOpen["frames.mp4"]

Galilean Transformation

Define β=v/c.  Add p to variable names instead of primes.  The Galilean transformation is

galileo={tp->ct,xp->(x-β ct),yp->y,zp->z}

Derivation of Lorentz Transformation

boost={
ctp->a00 ct + a01 x + a02 y + a03 z,
xp->a10 ct + a11 x + a12 y + a13 z,
yp->a20 ct + a21 x + a22 y + a23 z,
zp->a30 ct + a31 x + a32 y + a33 z}

Isotropy of space: events on positive and negative y- and z-axes occur at equal times in S'

ctp/.boost/.{x->0,y->y0,z->0}

isoy=Solve[(ctp/.boost/.{x->0,y->y0,z->0})==(ctp/.boost/.{x->0,y->-y0,z->0}),a02]

isoz=Solve[(ctp/.boost/.{x->0,y->0,z->z0})==(ctp/.boost/.{x->0,y->0,z->-z0}),a03]

Origin in S is mapped to

{ctp,xp,yp,zp}/.boost/.{x->0,y->0,z->0}

which must stay on the xp - axis

isoO=Solve[({yp,zp}/.boost/.{x->0,y->0,z->0})=={0,0},{a20,a30}]

Gather all simplifications so far:

iso=Flatten[{isoy,isoz,isoO}]

boost/.iso//TableForm

Cylindrical symmetry around the x-axis requires that each of the three axes remain parallel

cylyz={a33->a22,a12->0,a13->0,a21->0,a23->0,a31->0,a32->0}

TableForm[boost1=boost/.Flatten[{iso,cylyz}]]

Now a22 can only depend on the magnitude of β.  Considering the inverse transformation by the principle of relativity

boostIy={y->a22 yp}

(y/.boostIy/.boost1)==y

soly=Last[Solve[(y/. boostIy/. boost1)==y,a22]]

Redefine boost with all these simplifications

TableForm[boost1/.soly]

Assuming that the general transformation boost reduces to galileo in the limit of small β, we see that in that limit, we should have

smallβ={a00->1,a01->0,a10->-β ,a11->1}

(boost1//.Flatten[{iso,cylyz,soly}]/.smallβ)

boost1//.Flatten[{iso,cylyz,soly}]

invarianceInterval=Expand[ (ct^2)-(x^2)-(y^2)-(z^2)-((ctp^2)-(xp^2)-(yp^2)-(zp^2)/.boost1) ]

Coefficient[invarianceInterval,{(ct^2),ct x,(x^2)}]

rela=Map[(#==0)&,Coefficient[invarianceInterval,{(ct^2),ct x,(x^2)}]]

Solve[ rela,{a11,a01,a00}]//TableForm

Get a 4th equation from the motion of the origin in S' which is, by definition,

xp==-β c tp

Construct this in terms of the coordinates in S

movO=(xp==-β ctp/.boost/.{x->0,y->0,z->0})

Collect the 4 equations for 4 unknowns

Flatten[{rela,movO}]

TableForm[PowerExpand[Solve[Flatten[{rela,movO}],{a00,a01,a10,a11}],Assumptions->{β>0}]]

From

smallβ

we see that we should choose the solution

solxt=PowerExpand[Solve[Flatten[{rela,movO}],{a00,a01,a10,a11}],Assumptions->{β>0}][[3]]

Define the Lorentz transformation as a function of β as a set of rules:

TableForm[
Lorentz[β_]=boost//.Flatten[{iso,soly,cylyz,solxt}]
]

Simplify@Together@PowerExpand[{xp+ctp,xp-ctp}/. Lorentz[β],Assumptions->{β>0}]

γdef=γ->(1/(Sqrt[1-(β^2)]))

Simplify[Together[(Sqrt[-1+(γ^2)])/. γdef],Assumptions->{c>0,β>0}]

Lorentz[2/3]

Coefficient[ctp/.Lorentz[β],{ct,x,y}]

mL[β_]=(Coefficient[#1,{ct,x,y,z}]&)/@({ctp,xp,yp,zp}/. Lorentz[β])

mL[β]/.(β^2)->((-1+(γ^2))/(γ^2))//Together//PowerExpand//MatrixForm

Solve[(1/(Sqrt[1-b2]))==γ,b2]

Lorentz[β]/.(β^2)->(-1+(γ^2))/(γ^2)//Together//PowerExpand

mL[β].mL[-β]//Simplify//MatrixForm

Useful relations and plots

Plot[(1/(Sqrt[1-(β^2)])),{β,0,1},PlotRange->{0,10},GridLines->Automatic,AxesLabel->{β,γ},PlotStyle->Thick]

LogPlot[(1/(Sqrt[1-(β^2)])),{β,0.1,1},PlotRange->{0,10000},GridLines->Automatic,AxesLabel->{β,γ}]

Solve[β==((k-1)/(k+1)),k]

Solve[β==(((K^2)-1)/((K^2)+1)),K]

Solve[γ==((1/(Sqrt[1-(β^2)]))),β]

βf[γ_]=(Sqrt[1-(1/(γ^2))])

Rapidity and composition of Lorentz transformations (“addition” of velocity)

rapiditydef = β -> Tanh[ζ]

MatrixForm[mLr[ζ_]=mL[β]/.rapiditydef//Simplify//PowerExpand]

PowerExpand[Simplify[γ/. γdef/. rapiditydef]]

PowerExpand[Simplify[β γ/. γdef/. rapiditydef]]

Plot[ArcTanh[β],{β,-1,1},AxesLabel->{β,ζ},PlotRange->{-10,10},GridLines->Automatic]

ArcCosh[1]

How does γ combine?

Cosh[ζ1+ζ2]//TrigExpand

PowerExpand[Expand[Simplify[γ1 γ2 (1+βf[γ1] βf[γ2])]],Assumptions->{γ1>1}]

PowerExpand[Expand[Simplify[γ1 γ2 (1+βf[γ1] βf[γ2])]],Assumptions->{γ1>1}]

γ1 γ2+(Sqrt[((γ1)^2)-1]) (Sqrt[((γ2)^2)-1])

LogLinearPlot[ArcCosh[γ],{γ,1,10000},AxesLabel->{"γ",ζ},PlotRange->{0,10},GridLines->Automatic]

mLr[ζ1].mLr[ζ2]//Simplify//MatrixForm

mL[β1].mL[β2]//Simplify//MatrixForm

Change notation for MathType in PowerPoint

%/.{β1->(Subscript[β, 1]),β2->(Subscript[β, 2])}//MatrixForm

(mL[β1].mL[β2])[[1,1]]//Simplify

((mL[β1].mL[β2])[[1,1]]/.{β1->(Subscript[β, 1]),β2->(Subscript[β, 2])})==(mL[β][[1,1]])//Simplify

Solve[(mL[β1].mL[β2])[[1,1]]==(mL[β][[1,1]]),β]

PowerExpand[Factor[Solve[PowerExpand[(mL[β1].mL[β2])[[2,1]]]==mL[β][[2,1]],β]]]

PowerExpand[Factor[Solve[PowerExpand[(mL[β1].mL[β2])[[1,2]]]==mL[β][[1,2]],β]]]

((β1+β2)/(1+β1 β2))/.{β1->(Subscript[β, 1]),β2->(Subscript[β, 2])}

Plot3D[(((Subscript[β, 1])+(Subscript[β, 2]))/(1+(Subscript[β, 1]) (Subscript[β, 2]))),{(Subscript[β, 1]),0,1},{(Subscript[β, 2]),0,1},AxesLabel->{(Subscript[β, 1]),(Subscript[β, 2]),β},BoxRatios->{1, 1, 1},AxesStyle->Thick,PlotStyle->Opacity[.8]]

Plot with both signs

Plot3D[(((Subscript[β, 1])+(Subscript[β, 2]))/(1+(Subscript[β, 1]) (Subscript[β, 2]))),{(Subscript[β, 1]),-1,1},{(Subscript[β, 2]),-1,1},AxesLabel->{(Subscript[β, 1]),(Subscript[β, 2]),β},BoxRatios->{1, 1, 1},AxesStyle->Thick,PlotStyle->Opacity[.4]]

Plot3D[γ1 γ2+(Sqrt[-1+((γ1)^2)]) (Sqrt[-1+((γ2)^2)]),{γ1,1,10},{γ2,1,10},AxesLabel->{(Subscript[γ, 1]),(Subscript[γ, 2]),γ},BoxRatios->{1, 1, 1},AxesStyle->Thick,PlotStyle->Opacity[.8]]

D[(((Subscript[β, 1])+(Subscript[β, 2]))/(1+(Subscript[β, 1]) (Subscript[β, 2]))),(Subscript[β, 2])]//Simplify

Spacetime

RegionPlot[(t^2)-(x^2)<0,{x,-4,4},{t,-4,4},PlotLabel->"World line and light cone",ImageSize->Medium]

ContourPlot[Evaluate[Table[(t^2)-(x^2)==d,{d,{-4,-1,0,1,4}}]],{x,-4,4},{t,-4,4},ContourLabels->True, PlotLabel->"Isometric curves in (1+1)D spacetime"]

ContourPlot[(t^2)-(x^2),{x,-4,4},{t,-4,4},ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&),PlotLabel->"Isometric contours in (1+1)D spacetime"]

ContourPlot[(t^2)+(x^2),{x,-4,4},{t,-4,4},ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&),PlotLabel->"Isometric contours in Euclidean space",PlotLabel->{"x","y"}]

Two dimensional reduction of Lorentz transformation:

mL2[β_]=mL[β][[{1,2},{1,2}]]

mL2[14/15]//N

mL2[.9].{1,1}

mL2[.9].{1,-1}

mL2[.9].{-1,-1}

square={{-1,1},{1,1},{1,-1},{-1,-1},{-1,1}}

mL2[.5]

lorentzText[β_]:=Text["β="<>ToString[β]<>"  γ="<>ToString[(1/(Sqrt[1-(β^2)]))],{0.2,1.8},BaseStyle->Large]

ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==0,(t^2)-(x^2)==-4},{x,-4,4},{t,-4,4},PlotLabel->"Transformation of light-cone",Epilog->{Translate[lorentzText[-.5],{0,2}],Line[square],Red,Line[mL2[-.5].#&/@( square)]}]

CloudDeploy[ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==0,(t^2)-(x^2)==-4},{x,-4,4},{t,-4,4},PlotLabel->"Transformation of light-cone",Epilog->{Line[square],Red,Line[mL2[-.5].#&/@( square)]}]]

Table[ContourPlot3D[(t^2)-(x^2)-(y^2)==d,{x,-3,3},{y,-3,3},{t,-3,3}],{d,{-2,0,2}}]

With[{d=2},ContourPlot3D[{(t^2)-(x^2)-(y^2)==-(d^2),(t^2)-(x^2)-(y^2)==0,(t^2)-(x^2)-(y^2)==(d^2)},{x,-5.5,5.5},{y,-5.5,5.5},{t,-5,5},Mesh->4,ContourStyle->(Opacity[#]&/@{.4,.2,.41}),ViewPoint->{-5,-15,4},ImageSize->Large]]

With[{d=2},ContourPlot3D[{(t^2)-(x^2)-(y^2)==-(d^2),(t^2)-(x^2)-(y^2)==0,(t^2)-(x^2)-(y^2)==(d^2)},{x,-1.2,-1.2+10},{y,-5,5},{t,-5,5},ContourStyle->Opacity[0.27],Mesh->4,ViewPoint->{-15,-2,4},ImageSize->Large]]

Manipulate[ContourPlot3D[(t^2)-(x^2)-(y^2)==d,{x,-4,4},{y,-4,4},{t,-4,4}],{{d,0},-2,2}]

Lorentz contraction

Construct snapshots of a unit rod that is stationary in the lab frame.  The left end of the rod is at x=0.

Pick a snapshot at a time such that its right end will mapped to the xp-axis

mL2[-.5].{1,t}=={xp,0}

Solve[0.5773502691896258` +1.1547005383792517` t==0,t]

rod=Table[{{0,t},{1,t}},{t,{-1,-0.5,0,1}}]

Map[mL2[-.5].#&,rod,{2}]

With[{β=0},ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},
Epilog->{Thickness[0.01],Line/@rod,}]]

With[{β=-0.5},ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},
Epilog->{lorentzText[β],Thickness[0.01],Line/@rod,Red,Line/@Map[mL2[β].#&,rod,{2}]}]]

With[{β=0.5},ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},
Epilog->{lorentzText[β],Thickness[0.01],Line/@rod,Red,Line/@Map[mL2[β].#&,rod,{2}]}]]

Manipulate[ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},
Epilog->{lorentzText[β],Thickness[0.01],Line/@rod,Red,Line/@Map[mL2[β].#&,rod,{2}]}],{{β,0},-1,1}]

Export["LorentzContraction.mp4",Table [ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},
Epilog->{lorentzText[β],Thickness[0.01],Line/@rod,Red,Line/@Map[mL2[β].#&,rod,{2}]}],{β,-.9,.9,0.025}]]

CloudDeploy[Manipulate[ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},
Epilog->{lorentzText[β],Thickness[0.01],Line/@rod,Red,Line/@Map[mL2[β].#&,rod,{2}]}],{{β,0},-1,1}],"LorentzContraction",Permissions->"Public"]

Time dilation

(1/(Sqrt[1-(0.5^2)]))

duration=Table[{{x,0},{x,1}},{x,-1,1,1}];

ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},Epilog->{Thickness[0.01],Line/@duration}]

With[{β=-.5},ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},Epilog->{lorentzText[β],Thickness[0.01],Line/@duration,Red,Line/@Map[mL2[β].#&,duration,{2}]}]]

Reverse the velocity to illustrate relativity of simultaneity.

With[{β=.5},ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},Epilog->{lorentzText[β],Thickness[0.01],Line/@duration,Red,Line/@Map[mL2[β].#&,duration,{2}]}]]

Manipulate[ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==2,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},
PlotLabel->"Time dilation",Epilog->{lorentzText[β],Thickness[0.01],Line/@duration,Red,Line/@Map[mL2[β].#&,duration,{2}]}],{{β,0},-1,1}]

Export["TimeDilation.mp4",Table[ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==2,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},
PlotLabel->"Time dilation",Epilog->{lorentzText[β],Thickness[0.01],Line/@duration,Red,Line/@Map[mL2[β].#&,duration,{2}]}],{β,-.9,.9,0.05}]]

CloudDeploy[Manipulate[ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==2,(t^2)-(x^2)==1,(t^2)-(x^2)==0,(t^2)-(x^2)==-1,(t^2)-(x^2)==-4},{x,-2,2},{t,-2,2},
PlotLabel->"Time dilation",Epilog->{lorentzText[β],Thickness[0.01],Line/@duration,Red,Line/@Map[mL2[-β].#&,duration,{2}]}],{{β,0},-1,1}],"Time dilation",Permissions->"Public"]

grid=Flatten[Table[{t,x},{t,-3,3,.5},{x,-3,3,.5}],1];

Table[{Point[p],Red,Arrow[{p,mL2[-.5].p}]},{p,grid}]//Graphics

ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==0,(t^2)-(x^2)==-4},{x,-4,4},{t,-4,4},Epilog->{,Table[{Point[p],Red,},{p,grid}]},ImageSize->Medium]

With[{β=0.5},ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==0,(t^2)-(x^2)==-4},{x,-4,4},{t,-4,4},Epilog->{Text["β="<>ToString[β],{1,3.5},BaseStyle->Large],Table[{Point[p],Red,Arrow[{p,mL2[-β].p}]},{p,grid}]},ImageSize->Medium]]

Manipulate[ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==0,(t^2)-(x^2)==-4},{x,-4,4},{t,-4,4},Epilog->{Translate[lorentzText[β],{0,2}],Table[{Point[p],Red,Arrow[{p,mL2[-β].p}]},{p,grid}]},ImageSize->Medium],{{β,0},-.99,0.99}]

Export["LorentzFlow.mp4",Table[ContourPlot[{(t^2)-(x^2)==4,(t^2)-(x^2)==0,(t^2)-(x^2)==-4},{x,-4,4},{t,-4,4},Epilog->{Translate[lorentzText[β],{0,2}],Table[{Point[p],Red,Arrow[{p,mL2[-β].p}]},{p,grid}]},ImageSize->Medium],{β,-.9,0.9,0.05}]]

Lorentz[0.9]

Lorentz[14/15]

{x,t}/.{t->((15 (t-((14 x)/15)))/(Sqrt[29])),x->((15 (-((14 t)/15)+x))/(Sqrt[29])),y->y,z->z}

{((15 (-((14 t)/15)+x))/(Sqrt[29])),((15 (t-((14 x)/15)))/(Sqrt[29]))}/.{x->1,t->-1}//N

x+t/.Lorentz[β]//Simplify

(t^2)-(x^2)/.Lorentz[14/15]//Expand

{(t^2)-(x^2)==4,(t^2)-(x^2)==0,(t^2)-(x^2)==-4}/.Lorentz[0.5]//Simplify

ContourPlot[Evaluate[{(t^2)-(x^2)==4,(t^2)-(x^2)==0,(t^2)-(x^2)==-4}/.Lorentz[0.9]//Simplify],{x,-4,4},{t,-4,4}]

Mass, momentum, energy

D[m v[t] (1/(Sqrt[1-((v[t])^2)])),t]

Special relativity in accelerators

Two-body decay

Solve[{M==E1+E2,(E1^2)-(m1^2)==(E2^2)-(m2^2)},{E1,E2}]

Convert[(Sqrt[2m (Eb+m)])/.Eb->300 GeV/.m->(Subscript[m, p])(c^2)//N,GeV]

(Sqrt[(((Eb-m))/((Eb+m)))])/.Eb->300 GeV/.m->(Subscript[m, p])(c^2)//N//ToFundamentalSI

Masses of electrons, protons, ions and nuclei

?nuclearMassGeV

Reference sources, like the following function, give access to the atomic mass

IsotopeData[{82,208},"AtomicMass"]

UnitConvert[Quantity[207.976652071`12.,"AtomicMassUnit"],HoldForm[(("Electronvolts")/(("SpeedOfLight")^2))]]

To get the mass of the bare nucleus, we need to subtract the mass of the 82 electrons in the atom

IsotopeData[{82,208},"AtomicMass"]-82 ParticleData["Electron","Mass"]

In our package, we have a function that does that and returns

mPb=nuclearMassGeV[{82,208}]

((82. 6.37 TeV)/(nuclearMassGeV[{82,208}]))//ToFundamentalSI

(Subscript[m, e])

N[(Subscript[m, e])]

Convert[N[(Subscript[m, e])(c^2)],MeV]

LHC proton beam energy and momentum vs. French TGV train

Convert[N[(Subscript[m, p])(c^2)],GeV]

First, some typical parameters of the TGV Euroduplex train (Wikipedia)

TGVtrain={v->320. Kilo Meter/Hour,m->380. Tonne}

Contrary to the claim, and the impression you get when one goes by, TGVs are really quite slow

βTGV=ToFundamentalSI[N[(v/c)/.TGVtrain]]

γTGV=ToFundamentalSI[N[γ/.γdef/.β->v/c/.TGVtrain]]

Subtract 1 from that result to see the difference

%-1

So we can use the non-relativistic formula for the kinetic energy

KTGV=Convert[(1/2)m (v^2)/.TGVtrain,Mega Joule]

Momentum of a TGV

PTGV=ToFundamentalSI[m v/.TGVtrain]

The energy of a proton, number of protons in one “bunch” and the number of bunches in the LHC beam (original design values):

LHCprotonbeam={Eb->7. TeV,Nb->1.15 (10^11),kb->2808}

Gamma factor of a proton in LHC

ToFundamentalSI[N[(Eb/((Subscript[m, p]) (c^2)))/.LHCprotonbeam]]

So the LHC beam is extremely relativistic.

Energy in the beam

Eb Nb kb/.LHCprotonbeam

ELHCbeam=Convert[Eb Nb kb/.LHCprotonbeam,Mega Joule]

Momentum of the beam (highly relativistic)

N[Eb Nb kb /c/.LHCprotonbeam]

PLHCbeam=ToFundamentalSI[N[((Eb Nb kb)/c)/.LHCprotonbeam]]

(KTGV/ELHCbeam)

So the kinetic energy of the TGV  and the LHC beam are comparable.

(PTGV/PLHCbeam)

But the momentum of the TGV is >10 million times larger.   
The momentum of the LHC beam is comparable to that of a child’s toy car.
The LHC beam dump gets very hot but the impact of the beam is feeble.

This is the difference in relativistic dynamics p=E/c, instead of p=(Sqrt[2m E]), where c is has a large value in everyday units.

LEIR energy

mPb54=Convert[N[nuclearMassGeV[{82,208}]+54(Subscript[m, e])(c^2)],MeV]

LEIR={TinjN->4.2 MeV,TextN->72 MeV,mb->mPb54,A->208,Q->54,Z->82}

mPb54MeV=Convert[mPb54,MeV]/MeV

208 72

((Sqrt[(((mPb54MeV+208 72))^2)-(((mPb54MeV))^2)])/54)

%/208

((Sqrt[(((mPb54MeV+208 4.2))^2)-(((mPb54MeV))^2)])/54)

%/208

mPb54MeV+208 72

Plot[((Sqrt[(((mPb54MeV+208 TN))^2)-(((mPb54MeV))^2)])/54),{TN,4.2,72 },PlotLabel->"(Subscript[p, p]) in LEIR ramp",AxesLabel->{"(Subscript[T, N])/MeV","(Subscript[p, p])/MeV"},GridLines->Automatic,PlotStyle->Thick,ImageSize->Medium,AspectRatio->1,PlotRange->{0,Automatic}]

Plot[((Sqrt[(((mPb54MeV+208 TN))^2)-(((mPb54MeV))^2)])/((mPb54MeV+208 TN))),{TN,4.2,72 },PlotLabel->"Beam velocity in LEIR ramp",AxesLabel->{"(Subscript[T, N])/MeV","v/c"},GridLines->Automatic,PlotRange->{0,Automatic},PlotStyle->Directive[{Thick,Red}],ImageSize->Medium,AspectRatio->1]