-
Notifications
You must be signed in to change notification settings - Fork 20
/
OGRe.m
3195 lines (3018 loc) · 154 KB
/
OGRe.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
(* ::Package:: *)
(*
OGRe: An (O)bject-Oriented (G)eneral (Re)lativity Package for Mathematica
By Barak Shoshany ([email protected]) (baraksh.com)
GitHub repository: https://github.com/bshoshany/OGRe
Copyright (c) 2021 Barak Shoshany. Licensed under the MIT license.
If you use this package in published research, please cite it as follows:
* Shoshany, B., (2021). OGRe: An Object-Oriented General Relativity Package for Mathematica. Journal of Open Source Software, 6(65), 3416, https://doi.org/10.21105/joss.03416
You can also use the following BibTeX entry:
@article{Shoshany2021_OGRe,
author = {Barak Shoshany},
doi = {10.21105/joss.03416},
journal = {Journal of Open Source Software},
number = {65},
pages = {3416},
publisher = {The Open Journal},
title = {OGRe: An Object-Oriented General Relativity Package for Mathematica},
url = {https://doi.org/10.21105/joss.03416},
volume = {6},
year = {2021}
}
*)
BeginPackage["OGRe`"];
(* Check if the package has already been loaded, in case Get was used instead of Needs. *)
If[
ValueQ[OGRe`Private`AlreadyLoaded],
(* Then *)
(* Unprotect and clear all symbols, so they can be redefined. Useful for debugging, or to reload the package after updating. *)
Unprotect["OGRe`*"];
Unprotect["OGRe`Private`*"];
(* Keep the tensor objects and settings created so far during the session, so they don't get deleted when the package is reloaded. *)
OGReTemp`TensorData = OGRe`Private`TensorData;
ClearAll["OGRe`*"];
ClearAll["OGRe`Private`*"];
OGRe`Private`TensorData = OGReTemp`TensorData;
Remove["OGReTemp`*"];
OGRe`Private`AlreadyLoaded = True,
(* Else *)
OGRe`Private`AlreadyLoaded = True;
(* Initialize the symbol TensorData, which is used to store the data for the tensor objects, as well as user settings. This is done only on first load. *)
OGRe`Private`TensorData = Association[];
];
(* A dirty trick to make the package's public modules globally visible without defining their usage messages in advance. I prefer to define each usage message at the same time as the module itself, so it can also serve as documentation for the code. *)
Null[{
TAddCoordTransformation,
TCalc,
TCalcChristoffel,
TCalcEinsteinTensor,
TCalcGeodesicFromChristoffel,
TCalcGeodesicFromLagrangian,
TCalcGeodesicWithTimeParameter,
TCalcLagrangian,
TCalcNormSquared,
TCalcRicciScalar,
TCalcRicciTensor,
TCalcRiemannTensor,
TChangeDefaultCoords,
TChangeDefaultIndices,
TChangeID,
TChangeSymbol,
TCheckForUpdates,
TCite,
TCovariantD,
TDelete,
TDocs,
TExport,
TExportAll,
TGetComponents,
TImport,
TImportAll,
TInfo,
TLineElement,
TList,
TMessage,
TNewCoordinates,
TNewMetric,
TNewTensor,
TPartialD,
TSetAllowOverwrite,
TSetAssumptions,
TSetAutoUpdates,
TSetCurveParameter,
TSetIndexLetters,
TSetParallelization,
TSetReservedSymbols,
TShow,
TSimplify,
TVolumeElementSquared
}];
Begin["`Private`"]; (* OGRe`Private` *)
(* DO NOT change the format of the next line. It is used by TCheckForUpdates to detect the version of this file. Changing it will break the automatic update mechanism. Only change the version number and date. *)
OGReVersion = "v1.7.0 (2021-09-17)";
(* The raw URL of this file on GitHub. *)
OGReURL = "https://raw.githubusercontent.com/bshoshany/OGRe/master/OGRe.m";
(* If the global options LocalSymbol has been previously set, then it will have the Head Association. Otherwise it will have the Head LocalSymbol, and we create it now for later use. *)
If[
Head[LocalSymbol["OGReGlobalOptions"]] =!= Association,
(* Then *)
LocalSymbol["OGReGlobalOptions"] = Association[];
];
(* Load the global options from the persistent storage. *)
OGReGlobalOptions = LocalSymbol["OGReGlobalOptions"];
(* Set the "OGReVersion" key to the current version of the package. *)
OGReGlobalOptions["OGReVersion"] = OGReVersion;
(* If the option "AutoUpdates" has not been set, or is not a boolean value, we set it now to the default value of True. *)
If[
!KeyExistsQ[OGReGlobalOptions, "AutoUpdates"] || !BooleanQ[OGReGlobalOptions["AutoUpdates"]],
(* Then *)
OGReGlobalOptions["AutoUpdates"] = True;
];
(* If the option "AllowOverwrite" has not been set, or is not a boolean value, we set it now to the default value of False. *)
If[
!KeyExistsQ[OGReGlobalOptions, "AllowOverwrite"] || !BooleanQ[OGReGlobalOptions["AllowOverwrite"]],
(* Then *)
OGReGlobalOptions["AllowOverwrite"] = False;
];
(* Save the global options to the persistent storage. *)
LocalSymbol["OGReGlobalOptions"] = OGReGlobalOptions;
(* ===================================================
Public modules (accessible to the user) start here.
=================================================== *)
(* Create a nicely-formatted usage message. *)
CreateUsageMessage[f_, msg_String : {}] := Evaluate[f::usage] = ToString[TextCell[Row[Flatten[{List @@ StringReplace[msg, {"`" ~~ (x : Shortest[__]) ~~ "`" :> Style[x, Bold]}]}]]], StandardForm];
CreateUsageMessage[TAddCoordTransformation, "TAddCoordTransformation[`sourceID` \[Rule] `targetID`, `rules`] adds a transformation from the coordinate system `sourceID` to the coordinate system `targetID`.
`rules` must be a list of transformation rules. For example, {x \[Rule] r Sin[\[Theta]] Cos[\[Phi]], y \[Rule] r Sin[\[Theta]] Sin[\[Phi]], z \[Rule] r Cos[\[Theta]]} is a transformation from Cartesian to spherical coordinates."];
TAddCoordTransformation::ErrorRulesForm = "The transformation rules must be a list of rules of the form x \[Rule] y.";
TAddCoordTransformation::ErrorDifferentCoords = "The source and target coordinate systems must be different.";
TAddCoordTransformation::ErrorNotSameDim = "The source and target coordinate systems must be of the same dimension.";
TAddCoordTransformation[sourceID_String -> targetID_String, rules_List] := TAddCoordTransformation[sourceID, targetID, rules];
TAddCoordTransformation[sourceID_String, targetID_String, rules_List] := Module[
{
allJacobians,
ChristoffelJacobian,
dim,
i,
inverseJacobian,
j,
jacobian,
k,
newCoordSymbols,
oldCoordSymbols
},
(* Check that the tensor object sourceID exists. *)
CheckIfTensorExists[sourceID];
(* Check that the rules are of the correct form. *)
If[
!AllTrue[rules, MatchQ[#, _->_] &],
(* Then *)
Message[TAddCoordTransformation::ErrorRulesForm];
Abort[];
];
(* Check that both tensor objects represents coordinate systems. *)
CheckIfCoordinates[sourceID];
CheckIfCoordinates[targetID];
(* Check that the source and target coordinate systems are different. *)
If[
sourceID === targetID,
(* Then *)
Message[TAddCoordTransformation::ErrorDifferentCoords];
Abort[];
];
(* Check that the source and target coordinate systems are of the same dimension. *)
If[
Length[TensorData[sourceID]["Components"][{{1}, sourceID}]] != Length[TensorData[targetID]["Components"][{{1}, targetID}]],
(* Then *)
Message[TAddCoordTransformation::ErrorNotSameDim];
Abort[];
];
(* Add the transformation to the CoordTransformations key of the source object, or create it if it doesn't already exist. *)
If[
KeyExistsQ[TensorData[sourceID], "CoordTransformations"] && AssociationQ[TensorData[sourceID]["CoordTransformations"]],
(* Then *)
ChangeTensorKey[sourceID, "CoordTransformations", Append[TensorData[sourceID]["CoordTransformations"], targetID -> rules]],
(* Else *)
ChangeTensorKey[sourceID, "CoordTransformations", Association[targetID -> rules]];
];
(* Calculate the Jacobian, inverse Jacobian, and the "Christoffel Jacobian", i.e. the extra second-derivative term in the coordinate transformation of the Christoffel symbols, and store them in the tensor object of the source coordinates, to be used later whenever a coordinate transformation is performed. *)
oldCoordSymbols = TensorData[sourceID]["Components"][{{1}, sourceID}];
newCoordSymbols = TensorData[targetID]["Components"][{{1}, targetID}];
dim = Length[oldCoordSymbols];
jacobian = TensorSimplify[Table[
D[oldCoordSymbols[[i]] /. rules, newCoordSymbols[[j]]],
{i, 1, dim},
{j, 1, dim}
]];
inverseJacobian = TensorSimplify[Inverse[jacobian]];
ChristoffelJacobian = TensorSimplify[Table[
D[
oldCoordSymbols[[i]] /. rules,
newCoordSymbols[[j]],
newCoordSymbols[[k]]
],
{i, 1, dim},
{j, 1, dim},
{k, 1, dim}
]];
allJacobians = Association[
"Jacobian" -> jacobian,
"InverseJacobian" -> inverseJacobian,
"ChristoffelJacobian" -> ChristoffelJacobian
];
(* Add all three Jacobians to the Jacobians key of the source object, or create it if it doesn't already exist. *)
If[
KeyExistsQ[TensorData[sourceID], "Jacobians"] && AssociationQ[TensorData[sourceID]["Jacobians"]],
(* Then *)
ChangeTensorKey[sourceID, "Jacobians", Append[TensorData[sourceID]["Jacobians"], targetID -> allJacobians]],
(* Else *)
ChangeTensorKey[sourceID, "Jacobians", Association[targetID -> allJacobians]];
];
Return[sourceID];
];
DefaultResultID = "Result";
DefaultSymbol = "\[DottedSquare]";
CreateUsageMessage[TCalc, "TCalc[`formula`] calculates a tensor `formula`, which may involve any number of tensors in the format `ID`[`indices`], where `ID` is a tensor object and `indices` is a string representing the order of indices, along with any combination of the following operations:
\[Bullet] Addition: For example, \"A\"[\"\[Mu]\[Nu]\"] + \"B\"[\"\[Mu]\[Nu]\"].
\[Bullet] Contraction: For example, \"A\"[\"\[Mu]\[Lambda]\"] . \"B\"[\"\[Lambda]\[Nu]\"].
\[Bullet] Multiplication by scalar: For example, 2 * \"A\"[\"\[Mu]\[Nu]\"].
TCalc[`targetID`[`targetIndices`], `formula`, `symbol`] calculates a tensor `formula` and stores the result in a new tensor object.
`targetID` specifies the ID of the tensor object in which to store the result. If omitted, the ID \"" <> DefaultResultID <> "\" will be used.
`targetIndices` specifies the order of indices of the resulting tensor. The indices must be a permutation of the free indices of `formula`. If omitted, the indices are assumed to be in the same order as they appear in `formula`.
`symbol` specifies the symbol to use for the resulting tensor. If omitted, the placeholder symbol " <> DefaultSymbol <> " will be used."];
TCalc::ErrorIndices = "The LHS index specification \"`1`\" and the RHS index specification \"`2`\" must be the same up to permutation.";
TCalc::ErrorResult = "Invalid tensor expression obtained: `1`. Please check that the tensor expression you entered contains only tensor references of the form \"ID\"[\"indices\"] combined using addition, contraction (dot product), or multiplication by scalar."
TCalc[RHSExpression_, symbol_String : DefaultSymbol] := TCalc[DefaultResultID[""], RHSExpression, symbol];
TCalc[LHSTensorID_String, RHSExpression_, symbol_String : DefaultSymbol] := TCalc[LHSTensorID[""], RHSExpression, symbol];
TCalc[LHSTensorID_String[LHSIndices_String], RHSExpression_, symbol_String : DefaultSymbol] := Module[
{
allVars,
components,
LHSVars,
newComponents,
newIndices,
result,
resultID,
resultIndices,
RHSVars,
rules,
useCoords,
useIndices
},
(* Check that the tensor LHSTensorID doesn't already exist, but only if it's not the default ID. *)
If[
LHSTensorID =!= DefaultResultID,
(* Then *)
CheckIfOverwriting[LHSTensorID];
];
(* Define the rules for computing tensor formulas. *)
rules = {
(* Trace *)
ID_String[indices_String /; !DuplicateFreeQ[Characters[indices]]] :>
TensorTrace[ID[indices]],
(* Tensor addition *)
firstID_String[firstIndices_String] + secondID_String[secondIndices_String] :>
AddTensors[TensorTrace[firstID[firstIndices]], TensorTrace[secondID[secondIndices]]],
(* Multiplication of tensor by scalar *)
scalar_ * ID_String[indices_String] :>
TensorByScalar[TensorTrace[ID[indices]], scalar],
(* Contraction of tensors *)
firstID_String[firstIndices_String] . secondID_String[secondIndices_String] :>
ContractTensors[TensorTrace[firstID[firstIndices]], TensorTrace[secondID[secondIndices]]],
(* Partial derivative *)
TPartialD[derivativeIndex_String] . tensorID_String[tensorIndices_String] :>
DivOrGrad[derivativeIndex, TensorTrace[tensorID[tensorIndices]]],
(* Covariant derivative *)
TCovariantD[derivativeIndex_String] . tensorID_String[tensorIndices_String] :>
CovariantDivOrGrad[derivativeIndex, TensorTrace[tensorID[tensorIndices]]]
};
(* Repeatedly replace tensor operations with their results until we reach a fixed point. *)
result = ReplaceRepeated[RHSExpression, rules];
(* Check that the result is valid, i.e. of the form "tensorID"["indices"]. *)
If[
!MatchQ[result, _String[_String]],
(* Then *)
Message[TCalc::ErrorResult, result];
(* Clear the temporary tensors that were created for the purpose of the calculation. *)
ClearTemp[];
Abort[];
];
resultID = result[[0]];
resultIndices = result[[1]];
(* Get the indices, coordinates, and components of the result. *)
useIndices = TensorData[resultID]["DefaultIndices"];
useCoords = TensorData[resultID]["DefaultCoords"];
components = TensorData[resultID]["Components"][{useIndices, useCoords}];
(* Simplify the components. *)
components = TensorSimplify[components];
If[
LHSIndices === "",
(* Then *)
(* Either a scalar, or no rearranging of indices is desired. Store the result directly in a new tensor object. *)
SetTensorID[LHSTensorID, Association[
"Components" -> Association[{useIndices, useCoords} -> components],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> useIndices,
"Metric" -> TensorData[resultID]["Metric"],
"Role" -> "Calculated",
"Symbol" -> symbol
]],
(* Else *)
(* Check that the LHS and RHS index specifications are the same up to permutation. *)
If[
Sort[Characters[LHSIndices]] != Sort[Characters[resultIndices]],
(* Then *)
Message[TCalc::ErrorIndices, LHSIndices, resultIndices];
(* Clear the temporary tensors that were created for the purpose of the calculation. *)
ClearTemp[];
Abort[];
];
(* Collect the variables to be used for rearranging the indices. Both LHSVars and RHSVars will be the same set of variables, but potentially in a different order. *)
allVars = Association[];
Scan[(allVars[#] = Unique["var"]) &, Characters[LHSIndices]];
LHSVars = allVars[#]& /@ Characters[LHSIndices];
RHSVars = allVars[#]& /@ Characters[resultIndices];
(* Rearrange the components and indices to allow for a LHS with a different index order than the RHS. *)
newComponents = Table[
components[[Sequence @@ RHSVars]],
Evaluate[Sequence @@ ({#, 1, Length[components]} & /@ LHSVars)]
];
newIndices = Table[
useIndices[[StringPosition[resultIndices, Characters[LHSIndices][[n]]][[1, 1]]]],
{n, 1, StringLength[LHSIndices]}
];
(* Store the result in a new tensor object. *)
SetTensorID[LHSTensorID, Association[
"Components" -> Association[{newIndices, useCoords} -> newComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> newIndices,
"Metric" -> TensorData[resultID]["Metric"],
"Role" -> "Calculated",
"Symbol" -> symbol
]];
];
(* Print the explicit formula we calculated, including the correct index placement. TODO: Uncomment this in a future version, once this feature works properly. *)
(* OGRePrint["Calculated: ", TensorData[resultID]["Symbol"]]; *)
(* Clear the temporary tensors that were created for the purpose of the calculation. *)
ClearTemp[];
Return[LHSTensorID];
];
CreateUsageMessage[TCalcChristoffel, "TCalcChristoffel[`metricID`] calculates the Christoffel symbols (the coefficients of the Levi-Civita connection) from the metric `metricID` and stores the result in a new tensor object with ID \"`metricID`Christoffel\". Note that the Christoffel symbols are not the components of a tensor, but this tensor object will know to transform according to the correct rules under change of coordinates."];
TCalcChristoffel::ErrorNotMetric = "The Christoffel symbols can only be calculated from a tensor object representing a metric.";
TCalcChristoffel[metricID_String] := Module[
{
christoffelID,
inverseMetricID = NewTempID[]
},
(* Check that metricID exists. *)
CheckIfTensorExists[metricID];
(* Check that metricID is indeed a metric. *)
If[
TensorData[metricID]["Role"] =!= "Metric",
(* Then *)
Message[TCalcChristoffel::ErrorNotMetric];
Abort[];
];
(* Create a temporary tensor for the inverse metric, with two upper indices as its default configuration, to force the Christoffel symbols to have the correct index configuration. We do this to increase performance, since if we don't, then we'll have to raise the first index later, which is a costly operation. *)
SetTensorID[inverseMetricID, Association[
"Components" -> TensorData[metricID]["Components"],
"DefaultCoords" -> TensorData[metricID]["DefaultCoords"],
"DefaultIndices" -> {1, 1},
"Metric" -> TensorData[metricID]["Metric"],
"Role" -> "Temporary",
"Symbol" -> Superscript[TensorData[metricID]["Symbol"], "\[Lambda]\[Sigma]"]
]];
(* Calculate the Christoffel symbols, and give the tensor object the correct ID and symbol. *)
christoffelID = TCalc[
(metricID <> "Christoffel")["\[Lambda]\[Mu]\[Nu]"],
1/2 inverseMetricID["\[Lambda]\[Sigma]"] . (
TPartialD["\[Mu]"] . metricID["\[Nu]\[Sigma]"] +
TPartialD["\[Nu]"] . metricID["\[Sigma]\[Mu]"] -
TPartialD["\[Sigma]"] . metricID["\[Mu]\[Nu]"]
),
"\[CapitalGamma]"
];
(* Set the role of the tensor to Christoffel, so that OGRe will know to transform it as a connection and not as a tensor. *)
ChangeTensorKey[christoffelID, "Role", "Christoffel"];
Return[christoffelID];
];
CreateUsageMessage[TCalcEinsteinTensor, "TCalcEinsteinTensor[`metricID`] calculates the Einstein tensor from the metric `metricID` and stores the result in a new tensor object with ID \"`metricID`Einstein\". If a tensor with ID \"`metricID`RicciTensor\" exists, it will be assumed to be the Ricci tensor of the metric, and will be used in the calculation. Otherwise, \"`metricID`RicciTensor\" will be created using TCalcRicciTensor[]."];
TCalcEinsteinTensor::ErrorNotMetric = "The Einstein tensor can only be calculated from a tensor object representing a metric.";
TCalcEinsteinTensor[metricID_String] := Module[
{
EinsteinID
},
(* Check that metricID exists. *)
CheckIfTensorExists[metricID];
(* Check that metricID is indeed a metric. *)
If[
TensorData[metricID]["Role"] =!= "Metric",
(* Then *)
Message[TCalcEinsteinTensor::ErrorNotMetric];
Abort[];
];
(* If the Ricci tensor was not already calculated, calculate it now. *)
If[
!KeyExistsQ[TensorData, metricID <> "RicciTensor"],
(* Then *)
TCalcRicciTensor[metricID];
];
(* Calculate the Einstein tensor, and give the tensor object the correct ID and symbol. *)
EinsteinID = TCalc[
(metricID <> "Einstein")["\[Mu]\[Nu]"],
(metricID <> "RicciTensor")["\[Mu]\[Nu]"] - 1/2 metricID["\[Mu]\[Nu]"] . (metricID <> "RicciTensor")["\[Rho]\[Rho]"],
"G"
];
(* Set the role of the tensor to Einstein for future reference. *)
ChangeTensorKey[EinsteinID, "Role", "Einstein"];
Return[EinsteinID];
];
CreateUsageMessage[TCalcGeodesicFromChristoffel, "TCalcGeodesicFromChristoffel[`metricID`, `coordinatesID`] calculates the geodesic equations obtained for each of the coordinates in `coordinatesID` using the Christoffel symbols of the metric `metricID` and stores the result in a new rank-1 tensor object with ID \"`metricID`GeodesicFromChristoffel\". Equating the components to zero will yield the full system of geodesic equations.
The result will be given in terms of the coordinate symbols as functions of the curve parameter and their derivatives with respect to the curve parameter. The curve parameter can be selected using TSetCurveParameter[].
If `coordinatesID` is not specified, the default coordinate system of the metric will be used."]
TCalcGeodesicFromChristoffel::ErrorNotMetric = "The geodesic equation vector can only be calculated from a tensor object representing a metric.";
TCalcGeodesicFromChristoffel[metricID_String, coordinatesID_String : "_UseDefault_"] := Module[
{
accelComponents,
accelID,
christComponents,
christID,
coordSymbols,
newID = metricID <> "GeodesicFromChristoffel",
parameter = Symbol[TensorData[Options]["CurveParameter"]],
tangentComponents,
tangentID,
tempMetricComponents,
tempMetricID,
useCoords
},
(* Check that metricID exists. *)
CheckIfTensorExists[metricID];
(* Check that metricID is indeed a metric. *)
If[
TensorData[metricID]["Role"] =!= "Metric",
(* Then *)
Message[TCalcGeodesicFromChristoffel::ErrorNotMetric];
Abort[];
];
(* If a specific coordinate system is not given, use the metric's default coordinate system. *)
If[
coordinatesID === "_UseDefault_",
(* Then *)
useCoords = TensorData[metricID]["DefaultCoords"],
(* Else *)
(* Check that the tensor object coordinatesID exists and represents a coordinate system. *)
CheckIfTensorExists[coordinatesID];
CheckIfCoordinates[coordinatesID];
useCoords = coordinatesID;
];
(* If the Christoffel symbols were not already calculated, calculate them now. This should be done first, otherwise the temporary tensors will be deleted prematurely when TCalc finishes. *)
If[
!KeyExistsQ[TensorData, metricID <> "Christoffel"],
(* Then *)
TCalcChristoffel[metricID];
];
(* Make sure we have the components of the Christoffel symbols in the correct coordinate system. *)
AddRepresentation[metricID <> "Christoffel", {1, -1, -1}, useCoords];
coordSymbols = TensorData[useCoords]["Components"][{{1}, useCoords}];
(* Define a temporary metric with any instance of the coordinate symbols replaced with coordinate functions in terms of the curve parameter. For example, with coordinates {x, y, z}, the expression z * f[x, y] will be replaced with z[c] * f[x[c], y[c]] where c is the curve parameter. *)
tempMetricID = NewTempID[];
tempMetricComponents = TensorData[metricID]["Components"][{{-1, -1}, useCoords}] /. (# -> #[parameter] & /@ coordSymbols);
SetTensorID[tempMetricID, Association[
"Components" -> Association[{{-1, -1}, useCoords} -> tempMetricComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {-1, -1},
"Metric" -> tempMetricID,
"Role" -> "Temporary",
"Symbol" -> TensorData[metricID]["Symbol"]
]];
(* Define a temporary vector representing the tangent vector to the curve. The vector consists of the derivatives of the coordinates with respect to the curve parameter. The vector will be associated to the metric tempMetricID. *)
tangentID = NewTempID[];
tangentComponents = coordSymbols /. ((# -> #'[parameter]) & /@ coordSymbols);
SetTensorID[tangentID, Association[
"Components" -> Association[{{1}, useCoords} -> tangentComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {1},
"Metric" -> tempMetricID,
"Role" -> "Temporary",
"Symbol" -> OverDot[TensorData[useCoords]["Symbol"]]
]];
(* Define a temporary vector representing the acceleration of the curve. The vector consists of the second derivatives of the coordinates with respect to the curve parameter. The vector will be associated to the metric tempMetricID. *)
accelID = NewTempID[];
accelComponents = coordSymbols /. ((# -> #''[parameter]) & /@ coordSymbols);
SetTensorID[accelID, Association[
"Components" -> Association[{{1}, useCoords} -> accelComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {1},
"Metric" -> tempMetricID,
"Role" -> "Temporary",
"Symbol" -> Overscript[TensorData[useCoords]["Symbol"], "\[DoubleDot]"]
]];
(* Copy the Christoffel symbols to a new temporary tensor associated with the new metric and replace the coordinate symbols with coordinate functions of the curve parameters as above. *)
christID = NewTempID[];
christComponents = TensorData[metricID <> "Christoffel"]["Components"][{{1, -1, -1}, useCoords}] /. (# -> #[parameter] & /@ coordSymbols);
SetTensorID[christID, Association[
"Components" -> Association[{{1, -1, -1}, useCoords} -> christComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {1, -1, -1},
"Metric" -> tempMetricID,
"Role" -> "Christoffel",
"Symbol" -> "\[CapitalGamma]"
]];
(* Calculate the geodesic equation vector and give it the symbol 0 since it is equal to the zero vector. *)
TCalc[newID, accelID["\[Sigma]"] + christID["\[Sigma]\[Mu]\[Nu]"] . tangentID["\[Mu]"] . tangentID["\[Nu]"], "0"];
(* Delete the copy of the Christoffel symbols tensor. We did not mark it as temporary, to ensure that it is correctly treated as Christoffel symbols and not as a proper tensor, therefore TCalc will not remove it automatically. *)
RemoveTensorID[christID];
(* Change the geodesic equation vector's associated metric back to the original metric. *)
ChangeTensorKey[newID, "Metric", metricID];
(* Set the geodesic equation vector's role. *)
ChangeTensorKey[newID, "Role", "GeodesicFromChristoffel"];
Return[newID];
];
CreateUsageMessage[TCalcGeodesicWithTimeParameter, "TCalcGeodesicWithTimeParameter[`metricID`, `coordinatesID`] calculates the geodesic equations obtained for each of the coordinates in `coordinatesID` using the Christoffel symbols of the metric `metricID` and stores the result in a new rank-1 tensor object with ID \"`metricID`GeodesicFromChristoffel\". Equating the components to zero will yield the full system of geodesic equations.
The result will be given in terms of the spatial coordinate symbols as functions of the time coordinate and their derivatives with respect to time. The first coordinate of `coordinatesID` will be assumed to be the time coordinate, even if its symbol is not t.
If `coordinatesID` is not specified, the default coordinate system of the metric will be used."]
TCalcGeodesicWithTimeParameter::ErrorNotMetric = "The geodesic equation vector can only be calculated from a tensor object representing a metric.";
TCalcGeodesicWithTimeParameter[metricID_String, coordinatesID_String : "_UseDefault_"] := Module[
{
accelComponents,
accelID,
christComponents,
christID,
christWith0Components,
christWith0ID,
coordSymbols,
coordSymbolsWithoutTime,
newID = metricID <> "GeodesicWithTimeParameter",
parameter,
tangentComponents,
tangentID,
tempMetricComponents,
tempMetricID,
useCoords
},
(* Check that metricID exists. *)
CheckIfTensorExists[metricID];
(* Check that metricID is indeed a metric. *)
If[
TensorData[metricID]["Role"] =!= "Metric",
(* Then *)
Message[TCalcGeodesicWithTimeParameter::ErrorNotMetric];
Abort[];
];
(* If a specific coordinate system is not given, use the metric's default coordinate system. *)
If[
coordinatesID === "_UseDefault_",
(* Then *)
useCoords = TensorData[metricID]["DefaultCoords"],
(* Else *)
(* Check that the tensor object coordinatesID exists and represents a coordinate system. *)
CheckIfTensorExists[coordinatesID];
CheckIfCoordinates[coordinatesID];
useCoords = coordinatesID;
];
(* If the Christoffel symbols were not already calculated, calculate them now. This should be done first, otherwise the temporary tensors will be deleted prematurely when TCalc finishes. *)
If[
!KeyExistsQ[TensorData, metricID <> "Christoffel"],
(* Then *)
TCalcChristoffel[metricID];
];
(* Make sure we have the components of the Christoffel symbols in the correct coordinate system. *)
AddRepresentation[metricID <> "Christoffel", {1, -1, -1}, useCoords];
coordSymbols = TensorData[useCoords]["Components"][{{1}, useCoords}];
(* Use the first coordinate as the curve parameter. *)
parameter = coordSymbols[[1]];
(* For the replacements below we only want to replace the non-time coordinates. *)
coordSymbolsWithoutTime = coordSymbols[[2 ;;]];
(* Define a temporary metric with any instance of the coordinate symbols, except t, replaced with coordinate functions in terms of t. For example, with coordinates {t, x, y, z}, the expression z * f[x, y] will be replaced with z[t] * f[t, x[t], y[t]]. *)
tempMetricID = NewTempID[];
tempMetricComponents = TensorData[metricID]["Components"][{{-1, -1}, useCoords}] /. (# -> #[parameter] & /@ coordSymbolsWithoutTime);
SetTensorID[tempMetricID, Association[
"Components" -> Association[{{-1, -1}, useCoords} -> tempMetricComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {-1, -1},
"Metric" -> tempMetricID,
"Role" -> "Temporary",
"Symbol" -> TensorData[metricID]["Symbol"]
]];
(* Define a temporary vector representing the tangent vector to the curve. The vector consists of the derivatives of the coordinates with respect to t, thus the first component will be equal to 1. The vector will be associated to the metric tempMetricID. *)
tangentID = NewTempID[];
tangentComponents = coordSymbols /. ((# -> #'[parameter]) & /@ coordSymbolsWithoutTime);
tangentComponents[[1]] = 1;
SetTensorID[tangentID, Association[
"Components" -> Association[{{1}, useCoords} -> tangentComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {1},
"Metric" -> tempMetricID,
"Role" -> "Temporary",
"Symbol" -> OverDot[TensorData[useCoords]["Symbol"]]
]];
(* Define a temporary vector representing the acceleration of the curve. The vector consists of the second derivatives of the coordinates with respect to t, thus the first component will be equal to 0. The vector will be associated to the metric tempMetricID. *)
accelID = NewTempID[];
accelComponents = coordSymbols /. ((# -> #''[parameter]) & /@ coordSymbolsWithoutTime);
accelComponents[[1]] = 0;
SetTensorID[accelID, Association[
"Components" -> Association[{{1}, useCoords} -> accelComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {1},
"Metric" -> tempMetricID,
"Role" -> "Temporary",
"Symbol" -> Overscript[TensorData[useCoords]["Symbol"], "\[DoubleDot]"]
]];
(* Copy the Christoffel symbols to a new temporary tensor associated with the new metric and replace the coordinate symbols, except t, with coordinate functions of t as above. *)
christID = NewTempID[];
christComponents = TensorData[metricID <> "Christoffel"]["Components"][{{1, -1, -1}, useCoords}] /. (# -> #[parameter] & /@ coordSymbolsWithoutTime);
SetTensorID[christID, Association[
"Components" -> Association[{{1, -1, -1}, useCoords} -> christComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {1, -1, -1},
"Metric" -> tempMetricID,
"Role" -> "Christoffel",
"Symbol" -> "\[CapitalGamma]"
]];
(* For the geodesic equations in terms of t, we also need the Christoffel symbols with 0 in the first index, which is a rank-2 tensor. The following is just a temporary solution until I implement the option to do something like christID["0AB"] and get the components automatically. *)
christWith0ID = NewTempID[];
christWith0Components = christComponents[[1, All, All]];
SetTensorID[christWith0ID, Association[
"Components" -> Association[{{-1, -1}, useCoords} -> christWith0Components],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {-1, -1},
"Metric" -> tempMetricID,
"Role" -> "Temporary",
"Symbol" -> "\[CapitalGamma]"
]];
(* Calculate the geodesic equation vector and give it the symbol 0 since it is equal to the zero vector. *)
TCalc[newID, accelID["\[Sigma]"] + (christID["\[Sigma]\[Mu]\[Nu]"] - christWith0ID["\[Mu]\[Nu]"] . tangentID["\[Sigma]"]) . tangentID["\[Mu]"] . tangentID["\[Nu]"], "0"];
(* Delete the copy of the Christoffel symbols tensor. We did not mark it as temporary, to ensure that it is correctly treated as Christoffel symbols and not as a proper tensor, therefore TCalc will not remove it automatically. *)
RemoveTensorID[christID];
(* Change the geodesic equation vector's associated metric back to the original metric. *)
ChangeTensorKey[newID, "Metric", metricID];
(* Set the geodesic equation vector's role. *)
ChangeTensorKey[newID, "Role", "GeodesicWithTimeParameter"];
Return[newID];
];
CreateUsageMessage[TCalcGeodesicFromLagrangian, "TCalcGeodesicFromLagrangian[`metricID`, `coordinatesID`] calculates the geodesic equations obtained for each of the coordinates in `coordinatesID` using the curve Lagrangian of the metric `metricID` and stores the result in a new rank-1 tensor object with ID \"`metricID`GeodesicFromLagrangian\". Equating the components to zero will yield the full system of geodesic equations.
Derivatives with respect to the curve parameter in the Euler-Lagrange equation will be left unevaluated using Inactive[], which can sometimes help solve the geodesic equations by inspection. Use Activate[] to evaluate the derivatives.
The result will be given in terms of the coordinate symbols as functions of the curve parameter and their derivatives with respect to the curve parameter. The curve parameter can be selected using TSetCurveParameter[].
If `coordinatesID` is not specified, the default coordinate system of the metric will be used."]
TCalcGeodesicFromLagrangian::ErrorNotMetric = "The geodesic equation vector can only be calculated from a tensor object representing a metric.";
TCalcGeodesicFromLagrangian[metricID_String, coordinatesID_String : "_UseDefault_"] := Module[
{
coordSymbols,
EulerLagrangeComponents,
LagrangianComponents,
newID = metricID <> "GeodesicFromLagrangian",
parameter = Symbol[TensorData[Options]["CurveParameter"]],
useCoords
},
(* Check that metricID exists. *)
CheckIfTensorExists[metricID];
(* Check that metricID is indeed a metric. *)
If[
TensorData[metricID]["Role"] =!= "Metric",
(* Then *)
Message[TCalcGeodesicFromLagrangian::ErrorNotMetric];
Abort[];
];
(* If a specific coordinate system is not given, use the metric's default coordinate system. *)
If[
coordinatesID === "_UseDefault_",
(* Then *)
useCoords = TensorData[metricID]["DefaultCoords"],
(* Else *)
(* Check that the tensor object coordinatesID exists and represents a coordinate system. *)
CheckIfTensorExists[coordinatesID];
CheckIfCoordinates[coordinatesID];
useCoords = coordinatesID;
];
(* If the Lagrangian was not already calculated, calculate it now. *)
If[
!KeyExistsQ[TensorData, metricID <> "Lagrangian"],
(* Then *)
TCalcLagrangian[metricID, useCoords];
];
(* Make sure we have the components of the Lagrangian in the correct coordinate system. We divide them by 2 since geodesics calculated in this way will inevitably get additional factors of 2 from taking the derivatives of squares. *)
LagrangianComponents = AddRepresentation[metricID <> "Lagrangian", {}, useCoords][[1]] / 2;
coordSymbols = TensorData[useCoords]["Components"][{{1}, useCoords}];
(* Calculate the geodesic equation vector from the Lagrangian by calculating the Euler-Lagrange equation for each variable. We leave the derivative with respect to the curve parameter implicit using Inactive. *)
EulerLagrangeComponents = Table[D[LagrangianComponents, coordSymbols[[n]][parameter]] - Inactive[D][D[LagrangianComponents, coordSymbols[[n]]'[parameter]], parameter], {n, 1, Length[coordSymbols]}];
(* Store the result in a new tensor object and give it the symbol 0 since it is equal to the zero vector. *)
SetTensorID[newID, Association[
"Components" -> Association[{{1}, useCoords} -> TensorSimplify[EulerLagrangeComponents]],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {1},
"Metric" -> metricID,
"Role" -> "GeodesicFromLagrangian",
"Symbol" -> "0"
]];
Return[newID];
];
CreateUsageMessage[TCalcLagrangian, "TCalcLagrangian[`metricID`, `coordinatesID`] calculates the curve Lagrangian of the metric `metricID`, defined as the norm-squared of the tangent to the curve, and stores the result in a new tensor object with ID \"`metricID`Lagrangian\".
Taking the square root of (the absolute value of) the Lagrangian yields the integrand of the curve length functional. Varying the Lagrangian using the Euler-Lagrange equations yields the geodesic equations.
The result will be given in terms of the coordinate symbols as functions of the curve parameter and their derivatives with respect to the curve parameter. The curve parameter can be selected using TSetCurveParameter[].
If `coordinatesID` is not specified, the default coordinate system of the metric will be used."]
TCalcLagrangian::ErrorNotMetric = "The curve Lagrangian can only be calculated from a tensor object representing a metric.";
TCalcLagrangian[metricID_String, coordinatesID_String : "_UseDefault_"] := Module[
{
coordSymbols,
newID = metricID <> "Lagrangian",
parameter = Symbol[TensorData[Options]["CurveParameter"]],
tangentComponents,
tangentID,
tempMetricComponents,
tempMetricID,
useCoords
},
(* Check that metricID exists. *)
CheckIfTensorExists[metricID];
(* Check that metricID is indeed a metric. *)
If[
TensorData[metricID]["Role"] =!= "Metric",
(* Then *)
Message[TCalcLagrangian::ErrorNotMetric];
Abort[];
];
(* If a specific coordinate system is not given, use the metric's default coordinate system. *)
If[
coordinatesID === "_UseDefault_",
(* Then *)
useCoords = TensorData[metricID]["DefaultCoords"],
(* Else *)
(* Check that the tensor object coordinatesID exists and represents a coordinate system. *)
CheckIfTensorExists[coordinatesID];
CheckIfCoordinates[coordinatesID];
useCoords = coordinatesID;
];
coordSymbols = TensorData[useCoords]["Components"][{{1}, useCoords}];
(* Define a temporary metric with any instance of the coordinate symbols replaced with coordinate functions in terms of the curve parameter. For example, with coordinates {x, y, z}, the expression z * f[x, y] will be replaced with z[c] * f[x[c], y[c]] where c is the curve parameter. *)
tempMetricID = NewTempID[];
tempMetricComponents = TensorData[metricID]["Components"][{{-1, -1}, useCoords}] /. (# -> #[parameter] & /@ coordSymbols);
SetTensorID[tempMetricID, Association[
"Components" -> Association[{{-1, -1}, useCoords} -> tempMetricComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {-1, -1},
"Metric" -> tempMetricID,
"Role" -> "Temporary",
"Symbol" -> TensorData[metricID]["Symbol"]
]];
(* Define a temporary vector representing the tangent vector to the curve. The vector consists of the derivatives of the coordinates with respect to the curve parameter. The vector will be associated to the metric tempMetricID. *)
tangentID = NewTempID[];
tangentComponents = coordSymbols /. ((# -> #'[parameter]) & /@ coordSymbols);
SetTensorID[tangentID, Association[
"Components" -> Association[{{1}, useCoords} -> tangentComponents],
"DefaultCoords" -> useCoords,
"DefaultIndices" -> {1},
"Metric" -> tempMetricID,
"Role" -> "Temporary",
"Symbol" -> OverDot[TensorData[useCoords]["Symbol"]]
]];
(* Calculate the Lagrangian and give it the symbol L. *)
TCalc[newID, tangentID["\[Mu]"] . tangentID["\[Mu]"], "L"];
(* Change the Lagrangian's associated metric back to the original metric. *)
ChangeTensorKey[newID, "Metric", metricID];
(* Set the Lagrangian's role. *)
ChangeTensorKey[newID, "Role", "Lagrangian"];
Return[newID];
];
CreateUsageMessage[TCalcNormSquared, "TCalcNormSquared[`tensorID` calculates the norm-squared of the tensor `tensorID` with respect to its metric, that is, the tensor contracted with itself in all indices. For example, for a vector v^a the norm-squared will be v^a v_a and for a rank-2 tensor T^ab the result will be T^ab T_ab. The result is stored in a new rank-0 tensor object with ID \"`tensorID`NormSquared\"."];
TCalcNormSquared[tensorID_String] := Module[
{
indices,
NormSquaredID,
rank
},
(* Check that tensorID exists. *)
CheckIfTensorExists[tensorID];
(* Get the tensor's rank so we know how many indices to contract, and prepare the index string. *)
rank = Length[TensorData[tensorID]["DefaultIndices"]];
indices = StringTake[DefaultIndexLetters, rank];
(* Calculate the norm squared, and give the tensor object the correct ID. The symbol will be left as the default. *)
NormSquaredID = TCalc[
(tensorID <> "NormSquared")[""],
tensorID[indices] . tensorID[indices]
];
(* Set the role of the tensor to NormSquared for future reference. *)
ChangeTensorKey[NormSquaredID, "Role", "NormSquared"];
Return[NormSquaredID];
];
CreateUsageMessage[TCalcRicciScalar, "TCalcRicciScalar[`metricID` calculates the Ricci scalar from the metric `metricID` and stores the result in a new tensor object with ID \"`metricID`RicciScalar\". If a tensor with ID \"`metricID`RicciTensor\" exists, it will be assumed to be the Ricci tensor of the metric, and will be used in the calculation. Otherwise, \"`metricID`RicciTensor\" will be created using TCalcRicciTensor[]."];
TCalcRicciScalar::ErrorNotMetric = "The Ricci scalar can only be calculated from a tensor object representing a metric.";
TCalcRicciScalar[metricID_String] := Module[
{
RicciScalarID
},
(* Check that metricID exists. *)
CheckIfTensorExists[metricID];
(* Check that metricID is indeed a metric. *)
If[
TensorData[metricID]["Role"] =!= "Metric",
(* Then *)
Message[TCalcRicciScalar::ErrorNotMetric];
Abort[];
];
(* If the Ricci tensor was not already calculated, calculate it now. *)
If[
!KeyExistsQ[TensorData, metricID <> "RicciTensor"],
(* Then *)
TCalcRicciTensor[metricID];
];
(* Calculate the Ricci scalar, and give the tensor object the correct ID and symbol. *)
RicciScalarID = TCalc[
(metricID <> "RicciScalar")[""],
(metricID <> "RicciTensor")["\[Mu]\[Mu]"],
"R"
];
(* Set the role of the tensor to Ricci Scalar for future reference. *)
ChangeTensorKey[RicciScalarID, "Role", "Ricci Scalar"];
Return[RicciScalarID];
];
CreateUsageMessage[TCalcRicciTensor, "TCalcRicciTensor[`metricID`] calculates the Ricci tensor from the metric `metricID` and stores the result in a new tensor object with ID \"`metricID`RicciTensor\". If a tensor with ID \"`metricID`Riemann\" exists, it will be assumed to be the Riemann tensor of the metric, and will be used in the calculation. Otherwise, \"`metricID`Riemann\" will be created using TCalcRiemannTensor[]."];
TCalcRicciTensor::ErrorNotMetric = "The Ricci tensor can only be calculated from a tensor object representing a metric.";
TCalcRicciTensor[metricID_String] := Module[
{
RicciTensorID
},
(* Check that metricID exists. *)
CheckIfTensorExists[metricID];
(* Check that metricID is indeed a metric. *)
If[
TensorData[metricID]["Role"] =!= "Metric",
(* Then *)
Message[TCalcRicciTensor::ErrorNotMetric];
Abort[];
];
(* If the Riemann tensor was not already calculated, calculate it now. *)
If[
!KeyExistsQ[TensorData, metricID <> "Riemann"],
(* Then *)
TCalcRiemannTensor[metricID];
];
(* Calculate the Ricci tensor, and give the tensor object the correct ID and symbol. *)
RicciTensorID = TCalc[
(metricID <> "RicciTensor")["\[Mu]\[Nu]"],
(metricID <> "Riemann")["\[Lambda]\[Mu]\[Lambda]\[Nu]"],
"R"
];
(* Set the role of the tensor to Ricci Tensor for future reference. *)
ChangeTensorKey[RicciTensorID, "Role", "Ricci Tensor"];
Return[RicciTensorID];
];
CreateUsageMessage[TCalcRiemannTensor, "TCalcRiemannTensor[`metricID`] calculates the Riemann tensor from the metric `metricID` and stores the result in a new tensor object with ID \"`metricID`Riemann\". If a tensor with ID \"`metricID`Christoffel\" exists, it will be assumed to be the Christoffel symbols of the metric, and will be used in the calculation. Otherwise, \"`metricID`Christoffel\" will be created using TCalcChristoffel[]."];
TCalcRiemannTensor::ErrorNotMetric = "The Riemann tensor can only be calculated from a tensor object representing a metric.";
TCalcRiemannTensor[metricID_String] := Module[
{
RiemannID
},
(* Check that metricID exists. *)
CheckIfTensorExists[metricID];
(* Check that metricID is indeed a metric. *)
If[
TensorData[metricID]["Role"] =!= "Metric",
(* Then *)
Message[TCalcRiemannTensor::ErrorNotMetric];
Abort[];
];
(* If the Christoffel symbols were not already calculated, calculate them now. *)
If[
!KeyExistsQ[TensorData, metricID <> "Christoffel"],
(* Then *)
TCalcChristoffel[metricID];
];
(* Calculate the Riemann tensor, and give the tensor object the correct ID, symbol, and default index configuration. *)
RiemannID = TChangeDefaultIndices[
TCalc[
(metricID <> "Riemann")["\[Rho]\[Sigma]\[Mu]\[Nu]"],
TPartialD["\[Mu]"] . (metricID <> "Christoffel")["\[Rho]\[Nu]\[Sigma]"] -
TPartialD["\[Nu]"] . (metricID <> "Christoffel")["\[Rho]\[Mu]\[Sigma]"] +
(metricID <> "Christoffel")["\[Rho]\[Mu]\[Lambda]"] . (metricID <> "Christoffel")["\[Lambda]\[Nu]\[Sigma]"] -
(metricID <> "Christoffel")["\[Rho]\[Nu]\[Lambda]"] . (metricID <> "Christoffel")["\[Lambda]\[Mu]\[Sigma]"],
"R"
],
{1,-1,-1,-1}
];
(* Set the role of the tensor to Riemann for future reference. *)
ChangeTensorKey[RiemannID, "Role", "Riemann"];
Return[RiemannID];
];
CreateUsageMessage[TChangeDefaultCoords, "TChangeDefaultCoords[`tensorID`, `coordinatesID`] changes the default coordinate system of the tensor object `tensorID` to `coordinatesID`."];
TChangeDefaultCoords::ErrorCoordTensor = "Cannot change the default coordinate system for a tensor object representing a coordinate system."
TChangeDefaultCoords[tensorID_String, coordinatesID_String] := (
(* Check that the tensor objects sourceID and coordinatesID exist. *)
CheckIfTensorExists[tensorID];
CheckIfTensorExists[coordinatesID];
(* Check that the tensor object tensorID does not itself represents a coordinate system. *)
If[
TensorData[tensorID]["Role"] === "Coordinates",
(* Then *)
Message[TChangeDefaultCoords::ErrorCoordTensor];
Abort[];
];
(* Check that the tensor object coordinatesID represents a coordinate system. *)
CheckIfCoordinates[coordinatesID];
(* Add a representation to the tensor in the new coordinate system with the default indices, if it doesn't already exist. *)
AddRepresentation[tensorID, TensorData[tensorID]["DefaultIndices"], coordinatesID];
(* Change the DefaultCoords key. *)
ChangeTensorKey[tensorID, "DefaultCoords", coordinatesID];
Return[tensorID];
);
CreateUsageMessage[TChangeDefaultIndices, "TChangeDefaultIndices[`ID`, `indices`] changes the default index configuration of the tensor object `ID` to `indices`.
`indices` must be a list of the form {\[PlusMinus]1, \[PlusMinus]1, ...}, where +1 corresponds to an upper index and -1 corresponds to a lower index."];
TChangeDefaultIndices::ErrorCoords = "Cannot change the default index configuration for a tensor object representing a coordinate system."
TChangeDefaultIndices::ErrorMetric = "Cannot change the default index configuration for a tensor object representing a metric."
TChangeDefaultIndices::ErrorChristoffel = "Cannot change the default index configuration for a tensor object representing a Levi-Civita connection (Christoffel symbols)."
TChangeDefaultIndices[ID_String, indices_List] := (
(* Check that the tensor object ID exists. *)
CheckIfTensorExists[ID];
(* Check that the tensor object does not represent a coordinate system. *)
If[
TensorData[ID]["Role"] === "Coordinates",
(* Then *)
Message[TChangeDefaultIndices::ErrorCoords];
Abort[];
];
(* Check that the tensor object does not represent a metric. *)
If[
TensorData[ID]["Role"] === "Metric",
(* Then *)
Message[TChangeDefaultIndices::ErrorMetric];
Abort[];
];
(* Check that the tensor object does not represent a connection. *)
If[
TensorData[ID]["Role"] === "Christoffel",
(* Then *)
Message[TChangeDefaultIndices::ErrorChristoffel];
Abort[];
];
(* Check that the list of indices is of the correct form. *)
CheckIndicesForm[indices];
(* Add a representation to the tensor with the new indices in the default coordinate system, if it doesn't already exist. *)
AddRepresentation[ID, indices, TensorData[ID]["DefaultCoords"]];
(* Change the DefaultIndices key. *)
ChangeTensorKey[ID, "DefaultIndices", indices];
Return[ID];
);
CreateUsageMessage[TChangeID, "TChangeID[`oldID` \[Rule] `newID`] changes the ID of the tensor object `oldID` to `newID`.
If the tensor is a metric or a coordinate system, all currently defined tensors will be scanned, and any references to `oldID` will be replaced with `newID`."];
TChangeID[oldID_String -> newID_String] := TChangeID[oldID, newID];
TChangeID[oldID_String, newID_String] := (
(* Check that the tensor object oldID exists. *)
CheckIfTensorExists[oldID];
(* Check that the tensor object newID doesn't already exist. *)
CheckIfOverwriting[newID];
(* Copy the old tensor data to the new ID and then remove the old ID. *)
SetTensorID[newID, TensorData[oldID]];
RemoveTensorID[oldID];
(* If the tensor is a metric, replace all references to it. *)
If[
this["Role"] === "Metric",