forked from Obinna/xLightCone
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathxLightCone.m
2445 lines (1682 loc) · 137 KB
/
xLightCone.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:: *)
(************************************************************************)
(* 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. *)
(************************************************************************)
xAct`xLightCone`$Version={"0.0.1",{2014,09,29}};
xAct`xLightCone`$xTensorVersionExpected={"1.1.1",{2014,9,28}};
xAct`xLightCone`$xPertVersionExpected={"1.0.5",{2014,9,28}};
(* xLightCone: *)
(* Copyright (C) 2014- Obinna Umeh, Cyril Pitrou *)
(* This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License,or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place-Suite 330, Boston, MA 02111-1307,
USA.
*)
(* :Title: xLightCone *)
(* :Author: Obinna Umeh & Cyril Pitrou*)
(* :Context: xAct`xLightCone` *)
(* :Copyright: Obinna Umeh & Cyril Pitrou (2014-) *)
(* :Keywords: *)
(* :Source: xLightCone.nb *)
(* :Warning: *)
(* :Mathematica Version: 9.0 and later *)
(* :Limitations: *)
If[Unevaluated[xAct`xCore`Private`$LastPackage]===xAct`xCore`Private`$LastPackage,xAct`xCore`Private`$LastPackage="xAct`xLightCone`"];
Off[General::nostdvar]
Off[General::nostdopt]
BeginPackage["xAct`xLightCone`",{"xAct`xPert`","xAct`xTensor`","xAct`xPerm`","xAct`xCore`","xAct`ExpressionManipulation`"}]
If[Not@OrderedQ@Map[Last,{$xTensorVersionExpected,xAct`xTensor`$Version}],Throw@Message[General::versions,"xTensor",xAct`xTensor`$Version,$xTensorVersionExpected]]
If[Not@OrderedQ@Map[Last,{$xPertVersionExpected,xAct`xPert`$Version}],Throw@Message[General::versions,"xPert",xAct`xPert`$Version,$xPertVersionExpected]]
Print[xAct`xCore`Private`bars];
Print["Package xAct`xLightCone` version ",$Version[[1]],", ",$Version[[2]]];
Print["CopyRight (C) 2015-, Obinna Umeh under the General Public License."];
Off[General::shdw]
xAct`xLightCone`Disclaimer[]:=Print["These are points 11 and 12 of the General Public License:\n\nBECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM `AS IS\.b4 WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.\n\nIN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES."]
On[General::shdw]
If[xAct`xCore`Private`$LastPackage==="xAct`xLightCone`",
Unset[xAct`xCore`Private`$LastPackage];
Print[xAct`xCore`Private`bars];
Print["These packages come with ABSOLUTELY NO WARRANTY; for details type Disclaimer[]. This is free software, and you are welcome to redistribute it under certain conditions. See the General Public License for details."];
Print[xAct`xCore`Private`bars]];
$CovDFormat="Prefix";
Message[General::nostdvar,"$CovDFormat","Prefix"];
(*** VERSIONS ***)
$Version::usage="$Version is a global variable giving the version of the package xLightCone in use.";
$xTensorVersionExpected::usage="$xTensorVersionExpected is a global variable giving the oldest possible version of the package xTensor which is required by the version of the package xLightCone in use.";
$xPertVersionExpected::usage="$xPertVersionExpected is a global variable giving the oldest possible version of the package xPert which is required by the version of the package xLightCone in use.";
a::usage = "";
H::usage = "";
\[Theta]::usage = "";
\[ScriptK]::usage = "";
\[Phi]::usage = ".";
\[Psi]::usage = ".";
Bs::usage = ".";
Es::usage = ".";
Bvp::usage = ".";
Bvt::usage = ".";
Evp::usage = ".";
Evt::usage = ".";
Etpp::usage = ".";
Etpt::usage = ".";
Ett::usage = ".";
T::usage = ".";
Ls::usage = ".";
Lvp::usage = ".";
Lvt::usage = ".";
Vs::usage = ".";
Vvp::usage = ".";
Vvt::usage = ".";
V0::usage = ".";
Vspat::usage = ".";
\[Rho]::usage = ".";
P::usage = ".";
Conformal::usage = "";
ConformalWeight::usage = "";
DefMatterFields::usage = "";
DefMetricFields::usage = "";
DefScreenProjectedTensor::usage = "";
DefScreenSpaceMetric::usage = "";
ExtractComponents::usage = "";
$FirstOrderVectorPerturbations::usage = "";
$FirstOrderTensorPerturbations::usage = "";
InducedDecompositionLightCone::usage = "";
MyOrthogonalToVectorQ::usage = ""; (* Eventually this should be private *)
SetSlicingUpToScreenSpace::usage = "";
SetSlicingUpToScreenSpaceObinna::usage = "";
SplitMetric ::usage = "";
SplitMatter ::usage = "";
SplitPerturbations::usage = "";
ToInducedDerivativeScreenSpace::usage="";
ToLightConeFromRules::usage = "";
ToLightCone::usage = "";
(*$DebugInfoQ::usage = "";*)
$ConformalTime::usage ="";
$RadialLieDerivative::usage ="";
VisualizeTensorScreenSpace::usage ="";
$GeodesicEquation = "";
\[Epsilon]::usage = "\[Epsilon] is the default perturbation parameter. It is automatically defined when calling either the function DefMetricFields or the function DefMatterFields.";
Begin["xAct`xLightCone`Private`"]
(**Why a scalar head inside a scalar function?**)(**We choose that the Scalar head should be removed**)Unprotect[xAct`xTensor`NoScalar];
xAct`xTensor`NoScalar[f_?ScalarFunctionQ[expr_]]:=f[NoScalar@expr]
Protect[xAct`xTensor`NoScalar];
(**To avoid complaints from MakeRule when the rule has already been defined and the lhs is zero.**)
(*xAct`xTensor`MakeRule[{lhs_?(Evaluate[#]===0&),rhs_,conditions___},options:OptionsPattern[]]:={};*)
(*We create our own function based on MakeRule,but separating the special problematic case.This avoids overloading a definition from xTensor.It does not change the way xLightCone works,but that way,the code is clearly separated from xTensor.This is a much more polite way to modify the behaviour of xTensor.This method was advised by Leo Stein**)
SetAttributes[BuildRule,HoldFirst]
BuildRule[{lhs_?(Evaluate[#]===0&),rhs_,conditions___},options:OptionsPattern[]]:={};
BuildRule[{lhs_,rhs_,conditions___}]:=MakeRule[{lhs,rhs,conditions}];
BuildRule[{lhs_,rhs_,conditions___},options:OptionsPattern[]]:=MakeRule[{lhs,rhs,conditions},options];
(*In our case this definition is much better as we know there is no acceleration on the background.
Careful this is valid only in FL*)
Unprotect[OrthogonalToVectorQ];
OrthogonalToVectorQ[vector_][LieD[vector_?xTensorQ[i_]][expr_]]=.
OrthogonalToVectorQ[vector_][LieD[vector2_?xTensorQ[i_]][expr_]]:=OrthogonalToVectorQ[vector][expr]&&((*IndicesOf[Free,xAct`xTensor`Up][expr]===IndexList[]||*)LieD[vector[i]][vector2[-i]]===0);
(* This is very very Addhoc. It is because when there is a spatial derivative on a terme + a lie derivative, then OrthogonalToVectorQ is saying that it is not orthogonal because it cannot treat the LieDerivative case. This is a dirty cheat.*)
OrthogonalToVectorQ[vector_][expr1_+expr2_]:=(xAct`xTensor`Private`OToVcheck[vector,expr1+expr2, List@@FindFreeIndices[expr1+expr2]])||(OrthogonalToVectorQ[vector][expr2]&&OrthogonalToVectorQ[vector][expr2]);
Protect[OrthogonalToVectorQ];
(* We define our own OrthogonalToVectorQ where we add a rule for a list of vectors. Then if there is a list, it is distributed first*)
MyOrthogonalToVectorQ[vectors_List][tensor_?xTensorQ]:=xUpSet[
MyOrthogonalToVectorQ[vectors][tensor],And@@((OrthogonalToVectorQ[#][tensor]&)/@vectors)];
MyOrthogonalToVectorQ[vecs_List][expr_]:=((*Print["Hello!"];*)And@@((OrthogonalToVectorQ[#][expr]&)/@vecs));
MyOrthogonalToVectorQ[vec_][expr_]:=OrthogonalToVectorQ[vec][expr];
(***DEFAULT OPTIONS AND PROTECTED NAMES***)
(* We should clean that to make sure there are no usueless definitions *)
BackgroundFieldMethod=False;
$DebugInfoQ=False;
$CommutecdRules={};
$ConformalTime=True;
$FirstOrderVectorPerturbations=False;
$FirstOrderTensorPerturbations=False;
$OpenConstantsOfStructure=True;
$PrePrint=ScreenDollarIndices;
$SortCovDAutomatic=True;
Off[RuleDelayed::rhs];
$RadialLieDerivative=False;
$GeodesicEquation = False;
(*** COLORATION OF THE PERTURBATION ORDERS ***)
DerCharacter:=If[$ConformalTime,"\[Prime]","."];
RadialCharacter:="*"
(* The output format of the Lie derivative with respect to the vector normal to the hypersurfaces is a prime when one considers the conformal time, and a dot when one considers the cosmic time. *)
PerturbationOrderColor[1]:=RGBColor[0.85,0,0];
PerturbationOrderColor[2]:=RGBColor[0,0.6,0];
PerturbationOrderColor[3]:=RGBColor[0,0,0.85];
PerturbationOrderColor[4]:=RGBColor[0.6,0.4,0.2];
PerturbationOrderColor[5]:=RGBColor[0.7,0,0.7];
PerturbationOrderColor[_]:=RGBColor[1,0.6,0];
xTensorFormStop[Tensor]
MakeBoxes[Tens_?xTensorQ[LI[p_?((IntegerQ[#]&&#>=0)&)],LI[q_?((IntegerQ[#]&&#>=0)&)],LI[r_?((IntegerQ[#]&&#>=0)&)],inds___],StandardForm]:=
xAct`xTensor`Private`interpretbox[Tens[LI[p],LI[q],LI[r],inds],
(*If[(AnisotropyBool[SpaceType@InducedMetricOf[Tens]]||BianchiBool[SpaceType@InducedMetricOf[Tens]])&&Cases[{inds},_?UpIndexQ]=!={}&&q>=1,
Message[xPand::makeboxes];
RowBox[{OverscriptBox["\[Null]",RowBox[{"(",ToString[p],")"}]],"\[NegativeThinSpace]",OverscriptBox[MakeBoxes[Tens[inds],StandardForm],ToString[q]]}],*)
Switch[p,
0,RowBox[{OverscriptBox[MakeBoxes[Tens[inds],StandardForm],StringJoin@Join[Table[DerCharacter,{i,1,q}],Table[RadialCharacter,{j,1,r}]]]}],
_,RowBox[{OverscriptBox["\[Null]",RowBox[{StyleBox["("<>ToString[p]<>")",FontColor->PerturbationOrderColor[p]]}]],"\[NegativeThinSpace]",OverscriptBox[MakeBoxes[Tens[inds],StandardForm],StringJoin@Join[Table[DerCharacter,{i,1,q}],Table[RadialCharacter,{j,1,r}]]]}]
(*]*)
]
];
xTensorFormStart[Tensor]
(*** HANDLING EXPRESSIONS ***)
collect[expr_]:=Collect[expr,$PerturbationParameter,Identity]
(** fix by Jolyon Bloomfield. No Idea if this works or not. This should correct the bug found by Adam Solomon in August 2014 **)
FixScalar:=Scalar[expr_]:>Scalar[xAct`xTensor`Private`MathInputExpand[expr]]
org[expr_]:=Collect[ContractMetric[expr],$PerturbationParameter,ToCanonical[#/.FixScalar]&]
(* Old org before the fix*)
(* org[expr_]:=Collect[ContractMetric[expr],$PerturbationParameter,ToCanonical[#/.fixScalar]&]*)
TCnoCM[expr_]:=ToCanonical[expr,UseMetricOnVBundle->None]
(***PRIVATE BOOLEAN FUNCTIONS***)(**Miscellaneous**)$DefInfoQ=True;
(*If set to'False',the information messages are not printed.*)
(**Testing expressions**)
AnyIndicesListQ[inds_List]:=Cases[inds,AnyIndices[_]]=!={}
(*If the list of indices contains the head'AnyIndices',then AnyIndicesListQ[inds] returns'True';otherwise it returns'False'.*)
DefTensorQ[symb_]:=Cases[$Tensors,symb]==={symb}
(*If'symb' has been defined as a tensor,then DefTensorQ[symb] returns'True';otherwise it returns'False'.*)
GaugeQ[strg_]:=Cases[{"AnyGauge","FluidComovingGauge","ScalarFieldComovingGauge","FlatGauge","IsoDensityGauge","NewtonGauge","SynchronousGauge"},strg]==={strg}
(*If'strg' matches one of the element in the list,then GaugeQ[strg] returns'True';otherwise it returns'False'.*)
InducedMetricQ[symb_]:=If[MetricQ[symb],InducedFrom[symb]=!=Null,False]
InducedFromInducedMetricQ[symb_]:=False;
(*If'symb' is not defined as a metric,then InducedMetricQ[symb] returns'False'.Otherwise,it checks whether'symb' is an induced metric.If this is the case,then InducedMetricQ[symb] gives'True';otherwise it gives'False'.*)
SpaceTimeQ[strg_]:=Cases[{"FLCurved","FLFlat","Minkowski"},strg]==={strg}
(*If'strg' matches one of the element of the list,then SpaceTimeQ[strg] returns'True';otherwise it returns'False'.*)
TensorNullQ[tens_]:=If[Cases[SlotsOfTensor[tens],Labels]=!={},Print["** Warning: the function TensorNullQ is only suited to test tensors defined without label-indices. ",tens],With[{inds=DummyIn/@SlotsOfTensor[tens]},tens@@inds===0]]
(*TensorNullQ[tens] returns'True' if'tens' is a null tensor and'False' otherwise.This function is not suited to test tensors defined with label-indices.*)
(**Type of manifolds**)
(* Will probably disappear since we will only do curved and flat FL *)
FlatSpaceBool[Spacetype_]:=(Spacetype==="FLFlat"||Spacetype==="Minkowski")
CurvedSpaceBool[Spacetype_]:=(Spacetype==="FLCurved")
(*** CONVENIENT FUNCTIONS TO BUILD RULES ***)
PatternLeft[(left_:> right_),ElementsToBePatterned_List]:=(left/.((#->Pattern[#,_])&/@ ElementsToBePatterned)):>right
PatternLeft[(left_-> right_),ElementsToBePatterned_List]:=(left/.((#->Pattern[#,_])&/@ ElementsToBePatterned))->right
(* Note that only down indices can be passed to this function *)
IsElementInList[el_,list_]:=Length@Cases[list,el]>=1;
(*TranverseTensorQ[tens_]:=xTensorQ[tens]&&IsElementInList[Transverse,PropertiesList[tens]];*)
ScalarTensorQ[tens_]:=xTensorQ[tens]&&IsElementInList["Scalar",PropertiesList[tens]];
VectorTensorQ[tens_]:=xTensorQ[tens]&&IsElementInList["Vector",PropertiesList[tens]];
TensorTensorQ[tens_]:=xTensorQ[tens]&&IsElementInList["Tensor",PropertiesList[tens]];
(* Rules used to build Automaticrules. We first use MakeRule, and then we use these rule to put patterns on the Lable indices of the left hand side.*)
PatternTensorLeftScalar[(lhs_:>rhs_),TensDummy_?ScalarTensorQ[inds___]]:=Block[{TScal},
(Evaluate[lhs/.(TensDummy->PatternTest[Pattern[TensDummy,_],ScalarTensorQ])]:>(rhs))/.(TensDummy->TScal)
]
PatternTensorLeftVector[(lhs_:>rhs_),TensDummy_?VectorTensorQ[inds___]]:=Block[{TVect},
(Evaluate[lhs/.(TensDummy->PatternTest[Pattern[TensDummy,_],VectorTensorQ])]:>(rhs))/.(TensDummy->TVect)
]
PatternTensorLeftTensor[(lhs_:>rhs_),TensDummy_?TensorTensorQ[inds___]]:=Block[{TTens},
(Evaluate[lhs/.(TensDummy->PatternTest[Pattern[TensDummy,_],TensorTensorQ])]:>(rhs))/.(TensDummy->TTens)
]
(* ::Code:: *)
DefScreenSpaceMetric[metric_[inda_, indb_], Manifold_, cd2_, {cdpost_String, cdpre_String}, InducedHypersurface_, SpaceTimeType_?SpaceTimeQ] :=
(* Extracting the specifications of the problem (metric manifold normal vector etc...)*)
If[(MetricQ[First@InducedHypersurface]) && (xTensorQ[Last@InducedHypersurface]),
Module[{prot,
vbundle = VBundleOfIndex[inda],
h = First@InducedHypersurface,
n = Last@InducedHypersurface,
dim = DimOfManifold[Manifold],
indexlist = GetIndicesOfVBundle[VBundleOfIndex[inda], 3]},
With[{g = First@InducedFrom[First@InducedHypersurface],u = First@Rest@InducedFrom[First@InducedHypersurface],
ind1 = DummyIn[Tangent[Manifold]], ind2 = DummyIn[Tangent[Manifold]], ind3 = DummyIn[Tangent[Manifold]], ind4 = DummyIn[Tangent[Manifold]]},
(* Definition of the screen-space metric. Projected orthogonally to u and n.*)
DefTensor[metric[inda, indb], Manifold, If[$TorsionSign === 1, Symmetric[{inda, indb}], {}],
OrthogonalTo -> {u[-inda], u[-indb], n[-inda], n[-indb]}, ProjectedWith -> {h[-inda, ind1], h[-indb, ind1]}];
(* We defined a covariant derivative associated to this screen metric and we will load later all its properties. *)
DefCovD[cd2[inda], vbundle, {cdpost, cdpre}, OrthogonalTo -> {u[-inda], n[-inda]}, ProjectedWith -> {h[-inda, ind2], metric[-inda, ind2]}];
(* Even though the function DfCovD is called saying that the orthogonality should be wrt to u and n, wrt h and metric,
the function in xTensor is implemented sucha that only the orthogonality wrt u and h is implemented *)
(* For other references we need a function to extract the radial vector out of the screen space metric*)
(* TODO check if this is useful or not*)
RadialVectOrthToTheScreenSpace[metric] := RadialVectOrthToTheScreenSpace[metric] = n;
(* Warning and Error messages? I am not sure what is the purpose of this*)
OrthogonalVectors[x_] :=
If[x === h, Rest@InducedFrom[h],
If[x === metric, Flatten[{Rest@InducedFrom[h], Rest@{h, n}}],"\!\(" <> ToString[x] <>"\&-\) do not have orthogonal vectors "]];
InducedFromHyperSurface[x_] :=InducedFromHyperSurface[x] =
If[x === h, InducedFrom[h],
If[x === metric, {h, n},"\!\(" <> ToString[x] <> "\&-\) is not an induced metric"]];
(* We need to stecify that the trace of the screen space metric is not dim but dim-2.*)
(* We also ensure automatic contraction of the metric with itself *)
AutomaticRules[metric, BuildRule[{metric[ind3, ind1] metric[-ind3, -ind1], metric[ind1, -ind1]}]];
AutomaticRules[metric, BuildRule[{metric[-ind3, -indb] metric[ind3, ind1], metric[ind1, -indb]}, MetricOn -> All]];
AutomaticRules[metric, BuildRule[{metric[-ind1, ind1] , dim-2}]];
(* A few other automatic contractions *)
AutomaticRules[metric, BuildRule[{metric[ind3, ind2] g[-ind3, ind1], metric[ind1, ind2]}, MetricOn -> All]];
(*This line assigns further properties to the metric and the radial vector, We are following xTensor. *)
(* It calls a function which is implemented below*)
PropertiesOfInducedScreenSpaceMetric[metric[inda, indb], Manifold, cd2, {n, h, CovDOfMetric[h]}];
With[{KNSS=GiveSymbol[ExtrinsicK,metric]},
AutomaticRules[KNSS, BuildRule[{KNSS[ind3, ind2] g[-ind3, ind1], KNSS[ind1, ind2]}, MetricOn -> All]];
AutomaticRules[KNSS, BuildRule[{KNSS[ind3, ind2] g[ind1,-ind3], KNSS[ind1, ind2]}, MetricOn -> All]];
];
(*We need a series of obvious upvalues for the screen space metric. *)
(* These relations are to explain that the screen space metric is induced from the space and to say that it has a covariant derivative cd2*)
InducedFrom[metric] ^= {h, n};
VBundleOfMetric[metric] ^= VBundleOfMetric[g];
MetricOfCovD[cd2] ^= metric;
AppendTo[$Metrics, metric];
MetricQ[metric] ^= True;
InducedMetricQ[metric]^= True;
InducedFromInducedMetricQ[metric]^=True;
CovDOfMetric[metric] ^= cd2;
(* Ad hoc rule for OrthogonalToVectorQ. Not working TODO*)
Unprotect[OrthogonalToVectorQ];
OrthogonalToVectorQ[u][Evaluate[Projector[metric]][expr_]]:=True;
Protect[OrthogonalToVectorQ];
(* We also need to say that cd2 applied on the time vector is null*)
(* I have tried to build this rule with AutomaticRule/MakeRule but it fails...*)
(* Anyway this is an adhoc patch but if things were done correctly we would never have to put this rule*)
cd2[ind1_][u[ind2_]]:=0;
(* There is also this special definition which ensures tha the g metric can be contracted through cd2...*)
(*Special definitions suggested by Cyril*)(* This is patch which had beend added for cd but which is not added for cd2 since it has two master metrics (h and g)*)
(* So I need to cancel the Leibniz rule, define some particular cases and redefine the Leibniz rule*)
(* However now there is no need for such patch since this is already implemented in PropertiesOfInducedScreenSpaceMetric in the newest versions*)
(*
prot=Unprotect[cd2];
(* We undefine*)
cd2[i1_][x_ y_]=.;
cd2[_?GIndexQ][expr_]=.;
(* So as to be able to introduce a rule when there is the ultra master metric*)
cd2[i1_][g[a_?AIndexQ, b_?AIndexQ]] := 0;
cd2[i1_][x_ g[a_?AIndexQ, b_?AIndexQ]] := metric[a, b] cd2[i1][x]/; MyOrthogonalToVectorQ[{n,u}][x];
(* And we redefine *)
cd2[i1_][x_ y_]:=Module[{res},
res=Which[
(* Both are orthogonal in all their indices. We can use the Leibnitz rule *)
MyOrthogonalToVectorQ[{n,u}][x]&&MyOrthogonalToVectorQ[{n,u}][y],
cd2[i1][x]y+cd2[i1][y]x,
(* The expression is not globally orthogonal: complain and return unevaluated *)
Not@MyOrthogonalToVectorQ[{n,u}][x y],
Message[Validate::nonproj,x y];$Failed,
(* Expression is orthogonal, but factors are not. Avoid infinite recursion with this hack *)
FreeQ[{x,y},n],
(* CP Notice here that I do a double decomposition. First a 1+3 and then a 1+2 of the 3 space.*)
cd2[i1][Expand@GradNormalToExtrinsicK@Expand[InducedDecomposition[InducedDecomposition[x,{h,u}],{metric,n}]InducedDecomposition[InducedDecomposition[y,{h,u}],{metric,n}]]],
(* This should never happen *)
True,
$Failed
];
res/;res=!=$Failed
];
cd2[_?GIndexQ][expr_]:=$Failed/;Head[expr]=!=Times&&Not@MyOrthogonalToVectorQ[{n,u}][expr]&&Message[Validate::nonproj,expr];
(* OLD crap *)
(*cd2[i1_][x_ y_]:=cd2[i1][Expand@GradNormalToExtrinsicK@Expand[InducedDecomposition[x,{metric,n}]InducedDecomposition[y,{metric,n}]]]/;Head@x=!=g&&Head@y=!=g&&OrthogonalToVectorQ[n][x y];*)
(*cd2[i1_][x_ y_]:=cd2[i1][x]y+cd2[i1][y]x/;And[OrthogonalToVectorQ[n][x],OrthogonalToVectorQ[n][y]];*)
(*cd2[_?GIndexQ][expr_]:=Throw@Message[Validate::error,"Induced derivative acting on non-projected expression."]/;Not@OrthogonalToVectorQ[n][expr];*)
Protect[Evaluate[prot]];
(* ******* End of Patch *)
*)
(* Finally we ensure that the covd d2 we have define is the Levi Civita derivative associated woth the screen space metric*)
AutomaticRules[metric, BuildRule[{cd2[ind3][metric[ind1, ind2] ], 0}]];
(* A series of rules which states taht the derivative cd2 is induced (orthogonality to n of a cd2 applied to a prokected tensor)*)
(* These rules are implemented in xTensor and called when we define an induced derivative.
These were called automatically for h, but not for the screen space metric *)
xAct`xTensor`Private`MakeOrthogonalDerivative[cd2,metric[inda, -ind2],n[inda]];
xAct`xTensor`Private`MakeProjectedDerivative[cd2,metric[inda, -ind2],n[inda]];
(*g/:xAct`xTensor`Private`ZeroDerOnMetricQ[xAct`xTensor`Private`deronvbundle[cd2,TangentM],g]=.*)
(* We gather many many rules for the commutation of induced derivatives. These are in general made to make sure that the transverse conditions of vectors and tensors are used. It is also made to gather Laplacians.*)
(* TODO Corriger ca.*)
(* TODO. Problem the extrinsic curvature is automatically replace. So How can we choose this to be 0 in the geodesic equation?*)
(* Commente pour l'instant*)
(* TODO Gauss Codazzi should be done backwards. That is from Riemanncd2 to Riemanncd. How to do that?*)
(*
DummyS[Sym_]:=SymbolJoin[DS,Sym];
DummyV[Sym_]:=SymbolJoin[DV,Sym];
DummyT[Sym_]:=SymbolJoin[DT,Sym];
(* In general the tensors should not be assumed to be transverse. However they are symmetric and they should be perturbed.*)
Block[{Print},
DefScreenProjectedTensor[Evaluate[DummyS[metric]][],metric,SpaceTimesOfDefinition->{"Perturbed"}];
DefScreenProjectedTensor[Evaluate[DummyV[metric]][-ind1],metric,SpaceTimesOfDefinition->{"Perturbed"},TensorProperties->{(*"Transverse"*)}];
DefScreenProjectedTensor[Evaluate[DummyT[metric]][-ind1,-ind2],metric,TensorProperties->{"SymmetricTensor",(*"Transverse",*)"Traceless"},SpaceTimesOfDefinition->{"Perturbed"}];
];
$CommutecdRules=Join[$CommutecdRules,Flatten@Join[
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftScalar[#,Evaluate[DummyS[metric]][LI[p],LI[q],LI[r]]]&/@BuildRule[Evaluate[{cd2[-ind1][cd2[ind2][cd2[ind3][cd2[ind1][Evaluate[DummyS[h]][LI[p],LI[q],LI[r]]]]]],org[CommuteCovDs[CommuteCovDs[cd2[-ind1][cd2[ind2][cd2[ind3][cd2[ind1][Evaluate[DummyS[h]][LI[p],LI[q],LI[r]]]]]],cd2,{ind2,-ind1}],cd2,{ind3,-ind1}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftScalar[#,Evaluate[DummyS[metric]][LI[p],LI[q],LI[r]]]&/@BuildRule[Evaluate[{cd2[ind2][cd2[-ind1][cd2[ind3][cd2[ind1][Evaluate[DummyS[h]][LI[p],LI[q],LI[r]]]]]],org[CommuteCovDs[cd2[ind2][cd2[-ind1][cd2[ind3][cd2[ind1][Evaluate[DummyS[h]][LI[p],LI[q],LI[r]]]]]],cd2,{ind3,-ind1}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftScalar[#,Evaluate[DummyS[metric]][LI[p],LI[q],LI[r]]]&/@BuildRule[Evaluate[{cd2[-ind1][cd2[ind1][cd2[ind3][cd2[ind2][Evaluate[DummyS[h]][LI[p],LI[q],LI[r]]]]]],org[CommuteCovDs[CommuteCovDs[CommuteCovDs[CommuteCovDs[cd2[-ind1][cd2[ind1][cd2[ind3][cd2[ind2][Evaluate[DummyS[h]][LI[p],LI[q],LI[r]]]]]],cd2,{ind3,ind1}],cd2,{ind3,-ind1}],cd2,{ind2,ind1}],cd2,{ind2,-ind1}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftScalar[#,Evaluate[DummyS[metric]][LI[p],LI[q],LI[r]]]&/@BuildRule[Evaluate[{cd2[-ind2][cd2[-ind1][cd2[ind2][Evaluate[DummyS[h]][LI[p],LI[q],LI[r]]]]],org[CommuteCovDs[cd2[-ind2][cd2[-ind1][cd2[ind2][Evaluate[DummyS[h]][LI[p],LI[q],LI[r]]]]],cd2,{-ind1,-ind2}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftScalar[#,Evaluate[DummyS[metric]][LI[p],LI[q],LI[r]]]&/@BuildRule[Evaluate[{cd2[-ind2][cd2[ind2][cd2[-ind1][Evaluate[DummyS[h]][LI[p],LI[q],LI[r]]]]],org[CommuteCovDs[CommuteCovDs[cd2[-ind2][cd2[ind2][cd2[-ind1][Evaluate[DummyS[h]][LI[p],LI[q],LI[r]]]]],cd2,{-ind1,ind2}],cd2,{-ind1,-ind2}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[metric]][LI[p],LI[q],LI[r],ind2]]&/@BuildRule[Evaluate[{cd2[-ind2][cd2[-ind1][cd2[ind1][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind2]]]],org[CommuteCovDs[CommuteCovDs[cd2[-ind2][cd2[-ind1][cd2[ind1][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind2]]]],cd2,{-ind1,-ind2}],cd2,{ind1,-ind2}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[metric]][LI[p],LI[q],LI[r],ind2]]&/@BuildRule[Evaluate[{cd2[-ind2][cd2[-ind1][cd2[ind1][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind2]]]],org[CommuteCovDs[CommuteCovDs[cd2[-ind2][cd2[-ind1][cd2[ind1][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind2]]]],cd2,{-ind1,-ind2}],cd2,{ind1,-ind2}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[metric]][LI[p],LI[q],LI[r],ind3]]&/@BuildRule[Evaluate[{cd2[-ind3][cd2[ind1][cd2[ind2][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind3]]]],org[CommuteCovDs[CommuteCovDs[cd2[-ind3][cd2[ind1][cd2[ind2][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind3]]]],cd2,{ind1,-ind3}],cd2,{ind2,-ind3}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[metric]][LI[p],LI[q],LI[r],ind3]]&/@BuildRule[Evaluate[{cd2[ind1][cd2[-ind3][cd2[ind2][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind3]]]],org[CommuteCovDs[cd2[ind1][cd2[-ind3][cd2[ind2][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind3]]]],cd2,{ind2,-ind3}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[metric]][LI[p],LI[q],LI[r],ind2]]&/@BuildRule[Evaluate[{cd2[-ind3][cd2[-ind1][cd2[ind3][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind2]]]],org[CommuteCovDs[cd2[-ind3][cd2[-ind1][cd2[ind3][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind2]]]],cd2,{ind3,-ind1}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftVector[#,Evaluate[DummyV[metric]][LI[p],LI[q],LI[r],ind1]]&/@BuildRule[Evaluate[{cd2[-ind1][cd2[-ind2][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind1]]],org[CommuteCovDs[cd2[-ind1][cd2[-ind2][Evaluate[DummyV[h]][LI[p],LI[q],LI[r],ind1]]],cd2,{-ind2,-ind1}]]}]]),
PatternLeft[#,{p,q,r}]&/@(PatternTensorLeftTensor[#,Evaluate[DummyT[metric]][LI[p],LI[q],LI[r],ind2,ind3]]&/@BuildRule[Evaluate[{cd2[-ind2][cd2[-ind1][Evaluate[DummyT[h]][LI[p],LI[q],LI[r],ind2,ind3]]],org[CommuteCovDs[cd2[-ind2][cd2[-ind1][Evaluate[DummyT[h]][LI[p],LI[q],LI[r],ind2,ind3]]],cd2,{-ind1,-ind2}]]}]])
]];
Block[{Print},
UndefTensor[Evaluate[DummyS[metric]]];
UndefTensor[Evaluate[DummyV[metric]]];
UndefTensor[Evaluate[DummyT[metric]]];
];
*)
]];
, Print["** DefMetric:: You have to ensure first that the following objects are defined: metric induced from the super metric, a hypersurface specifying four vector and a \
screen space specifying vector"]];
(*I am not quite sure for now about the sign of ExtrinsicKSign on the \
screeen space, so i have defined the following for now, we may end up \
merging with xTensor*)
$ExtrinsicKOnSSSign = $ExtrinsicKSign;
$AccelerationOfnSign = $AccelerationSign;
(*This function is called in DefScreenSpaceMetric and sets various properties for the screen space metric. It is adapted from xTensor*)
(* Maybe it is useless to copy everything here.*)
(* Maybe that would be enough to call DefInducedMetric. To be debated*)
PropertiesOfInducedScreenSpaceMetric[metric_[-ind1_, -ind2_],
dependencies_, covd_, {vector_, supermetric_, superCD_}] :=
With[{vbundle = VBundleOfIndex[ind1]},
With[{extrinsicKname = GiveSymbol[ExtrinsicK, metric],
accelerationname = GiveSymbol[Acceleration, vector],
projectorname = GiveSymbol[Projector, metric],
superprojectorname = GiveSymbol[Projector, supermetric],
epsilonname = GiveSymbol[epsilon, metric],
superepsilonname = GiveSymbol[epsilon, supermetric],
proj = ProjectWith[metric],
u=Last@InducedFrom@supermetric,
g=First@InducedFrom@supermetric,
norm = Scalar@Simplify@ContractMetric[supermetric[-ind1, -ind2] vector[ind1] vector[ind2],supermetric], indexlist = GetIndicesOfVBundle[vbundle, 3]},
With[{i1 = indexlist[[1]], i2 = indexlist[[2]],
i3 = indexlist[[3]]},(*Register pair metric/vector*)
xUpSet[VectorOfInducedMetric[metric], vector];
(* Definition of extrinsic curvature *)
DefTensor[extrinsicKname[i1, i2], dependencies,
Symmetric[{1, 2}],
PrintAs :> GiveOutputString[ExtrinsicK, metric],
OrthogonalTo -> {vector[-i1], u[-i1](* I add this as well*)},
ProjectedWith -> {metric[i3, -i2], supermetric[i3,-i2](*Here I add this*)}, ProtectNewSymbol -> False,
Master -> metric,
DefInfo -> {"extrinsic curvature tensor",
"Associated to vector " <> ToString[vector]},
TensorID -> {ExtrinsicK, metric}];
(* Definition of Acceleration *)
DefTensor[accelerationname[i1], dependencies,
PrintAs :> GiveOutputString[Acceleration, vector],
OrthogonalTo -> {vector[-i1]},
ProjectedWith -> {metric[i2, -i1]}, ProtectNewSymbol -> False,
Master -> metric,
DefInfo -> {"acceleration vector",
"Associated to vector " <> ToString[vector]},
TensorID -> {Acceleration, vector}];
(*Relations among them and the derivatives.Improved by Thomas,
to use HasOrthogonalIndexQ*)
xAct`xTensor`Private`GradNormalToExtrinsicKRules[metric] = {superCD[a_][vector[b_]] :> $ExtrinsicKOnSSSign extrinsicKname[a, b]
+ $AccelerationOfnSign vector[a] accelerationname[b],
vector[-a_] superCD[b_][expr_] :> -superCD[b][vector[-a]] expr /;
xAct`xTensor`Private`HasOrthogonalIndexQ[expr, vector[-a]],
vector[a_] superCD[b_][expr_] :> -superCD[b][vector[a]] expr /;
xAct`xTensor`Private`HasOrthogonalIndexQ[expr, vector[a]],
LieD[vector[_]][expr_] vector[-a_] :> -$AccelerationOfnSign norm accelerationname[-a] expr /; xAct`xTensor`Private`HasOrthogonalIndexQ[expr, vector[-a]](*,
LieD[vector[_]][expr_]vector[a_]:>0/;
HasOrthogonalIndexQ[expr,vector[a]]*)};
xAct`xTensor`Private`ExtrinsicKToGradNormalRules[metric] =
extrinsicKname[a_, b_] :>
Module[{c =
DummyIn@vbundle}, $ExtrinsicKOnSSSign (supermetric[a,
c] superCD[-c][vector[b]] - $AccelerationOfnSign vector[
a] accelerationname[b])];
(*Projectors and metrics*)
xAct`xTensor`Private`ProjectorToMetricRules[metric] =
metric[i1_, i2_] ->
supermetric[i1, i2] - vector[i1] vector[i2]/norm;
xAct`xTensor`Private`MetricToProjectorRules[metric] =
supermetric[i1_, i2_] ->
metric[i1, i2] + vector[i1] vector[i2]/norm;
(*Define projector inert-head*)
DefInertHead[projectorname, LinearQ -> True, Master -> metric,
PrintAs :> GiveOutputString[Projector, metric],
ProtectNewSymbol -> False,
DefInfo -> {"projector inert-head", ""}];
(* I add this contraction of the screen projector with h*)
xTagSet[{projectorname, ContractThroughQ[projectorname, supermetric]},True];
projectorname[supermetric[a_, b_]] := metric[a, b];
(*The metric, but not the supermetric, can be contracted through the projector*)
(* TODO CHeck that these contraction rules are consistent *)
xTagSet[{projectorname, ContractThroughQ[projectorname, metric]},True];
(*The supermetric is converted into metric when contracted with the projector or covd*)
xTagSetDelayed[{projectorname,supermetric[i1_, i2_] projectorname[expr_]},
metric[i1, i2] projectorname[expr]/;Or[IsIndexOf[expr, -i1, metric], IsIndexOf[expr, -i2, metric]]];
xTagSetDelayed[{covd, supermetric[i1_, i2_] covd[i3_][expr_]},
metric[i1, i2] covd[i3][expr]/;Or[IsIndexOf[covd[i3][expr], -i1, metric],IsIndexOf[covd[i3][expr], -i2, metric]]];
(*Projection rule with vector*)
xTagSetDelayed[{projectorname, vector[i_] projectorname[expr_]},0 /; IsIndexOf[expr, -i, metric]];
(* CP We add the orthogonality with respect to the master vector u*)
xTagSetDelayed[{projectorname, u[i_] projectorname[expr_]},0 /; IsIndexOf[expr, -i, metric]];
(*This should be checked that it works. TODO*)
(*Particular cases*)
projectorname[1] := 1;
projectorname[rest_. x_?ScalarQ] := Scalar[x] projectorname[rest];
projectorname[vector[i_] expr_.] := 0 /; Not@IsIndexOf[expr, -i, supermetric];
(* CP I add this rule because now this projector is orthogonal to both vectors. *)
projectorname[u[i_] expr_.] := 0 /; Not@IsIndexOf[expr, -i, supermetric];
projectorname[projectorname[expr_]] := projectorname[expr];
projectorname[tensor_?xTensorQ[inds__]] := tensor[inds] /; MyOrthogonalToVectorQ[{vector,u}][tensor];
(* CP: Here we specify that it should be orthogonal to both vector to deserve the removal of the Projector Head.*)
(* This lead to a bug previously...*)
(* CP: I add this rule. It should be OK*)
projectorname[superprojectorname[expr_]] := projectorname[expr];
projectorname[covd[k_][expr_]] := covd[k][expr];
xAct`xTensor`Private`ProjectDerivativeRules[
covd] = {covd[i_][expr_] :>
If[IsIndexOf[expr, -i],
With[{dummy = DummyAs[i]},
metric[i, -dummy] projectorname[superCD[dummy][expr]]],
projectorname[superCD[i][expr]]]};
Module[{prot = Unprotect[covd]},(*Leibnitz rule.Three cases considered*)
covd[i1_][x_Scalar y_Scalar] := x covd[i1][y] + y covd[i1][x];
(*Special definitions suggested by Cyril*)
covd[i1_][supermetric[a_?AIndexQ, b_?AIndexQ]] := 0;
(*Cyril suggests removing the MyOrthogonalToVectorQ check to handle the many-supermetrics case*)
covd[i1_][x_ supermetric[a_?AIndexQ, b_?AIndexQ]] := metric[a, b] covd[i1][x] /; MyOrthogonalToVectorQ[{u,vector}][x];
(* CP ANd we add the same type of rule for the supersupermetric.*)
covd[i1_][g[a_?AIndexQ, b_?AIndexQ]] := 0;
covd[i1_][x_ g[a_?AIndexQ, b_?AIndexQ]] := metric[a, b] covd[i1][x]/; MyOrthogonalToVectorQ[{vector,u}][x];
covd[i1_][x_ delta[a_?AIndexQ, b_?AIndexQ]] := metric[a, b] covd[i1][x] /; MyOrthogonalToVectorQ[{u,vector}][x];
covd[i1_][vector[a_?AIndexQ] vector[b_?AIndexQ] x_.] :=
0 /; And[!IsIndexOf[x, ChangeIndex@a], !IsIndexOf[x, ChangeIndex@b], !PairQ[a, b]];
(* CP I add this. This should be consistent. Why not? *)
covd[i1_][u[a_?AIndexQ] vector[b_?AIndexQ] x_.] :=
0 /; And[!IsIndexOf[x, ChangeIndex@a], !IsIndexOf[x, ChangeIndex@b], !PairQ[a, b]];
covd[i1_][u[a_?AIndexQ] u[b_?AIndexQ] x_.] :=
0 /; And[!IsIndexOf[x, ChangeIndex@a], !IsIndexOf[x, ChangeIndex@b], !PairQ[a, b]];
(*Product of two,perhaps contracted,expressions*)
covd[i1_][x_ y_] :=
Module[{res}, res = Which[(*Both are orthogonal in all their indices.We can use the Leibnitz rule*)
MyOrthogonalToVectorQ[{vector,u}][x] && MyOrthogonalToVectorQ[{vector,u}][y],
covd[i1][x] y + covd[i1][y] x,(*The expression is not globally orthogonal: complain and return unevaluated*)
Not@MyOrthogonalToVectorQ[{vector,u}][x y],
Message[Validate::nonproj, x y]; $Failed,(*Expression is orthogonal, but factors are not.Avoid infinite recursion with this hack*)
FreeQ[{x, y}, vector],
covd[i1][Expand@GradNormalToExtrinsicK@Expand[InducedDecomposition[InducedDecomposition[x,{supermetric,u}], {metric, vector}] InducedDecomposition[InducedDecomposition[y,{supermetric,u}], {metric, vector}]]],
(*This should never happen*)True, $Failed];
res /; res =!= $Failed];
(*Induced derivatives of non-spatial objects are not accepted,
not even divergencies*)
covd[_?GIndexQ][expr_] := $Failed /;Head[expr] =!= Times && Not@MyOrthogonalToVectorQ[{vector,u}][expr] && Message[Validate::nonproj, expr];
Protect[Evaluate[prot]];];
(*Special definitions*)
metric /:LieD[vector[_]][metric[-a_Symbol, -b_Symbol]] := $ExtrinsicKOnSSSign ToCanonical[(extrinsicKname[-a, -b] + extrinsicKname[-b, -a]),UseMetricOnVBundle->None];
metric /:LieD[vector[_]][metric[a_Symbol, -b_Symbol]] := -$AccelerationOfnSign accelerationname[-b] vector[a];
metric /: LieD[vector[_]][metric[-a_Symbol,b_Symbol]] := -$AccelerationOfnSign accelerationname[-a] vector[b];
metric /: LieD[vector[_]][metric[a_Symbol,b_Symbol]] := -$ExtrinsicKOnSSSign ToCanonical[(extrinsicKname[a, b] +
extrinsicKname[b,a]),UseMetricOnVBundle->None] - $AccelerationOfnSign (vector[a] accelerationname[b] + vector[b] accelerationname[a]);
vector /:LieD[vector[_]][vector[a_Symbol]] := 0;
vector /:LieD[vector[_]][vector[-a_Symbol]] := $AccelerationOfnSign norm accelerationname[-a];
Module[{prot = Unprotect[{superepsilonname, epsilonname}]},
superepsilonname /:
LieD[vector[_]][superepsilonname[inds__?DownIndexQ]] :=
Module[{dummy = DummyIn[vbundle]}, $ExtrinsicKOnSSSign extrinsicKname[dummy, -dummy] superepsilonname[inds]];
superepsilonname /:
LieD[vector[_]][superepsilonname[inds__?UpIndexQ]] :=
Module[{dummy = DummyIn[vbundle],first = First[{inds}]}, -$ExtrinsicKOnSSSign extrinsicKname[dummy, -dummy] superepsilonname[inds]];
epsilonname /: LieD[vector[_]][epsilonname[inds__?DownIndexQ]] :=
Module[{dummy = DummyIn[vbundle]}, $ExtrinsicKOnSSSign extrinsicKname[dummy, -dummy] epsilonname[inds]];
epsilonname /: LieD[vector[_]][epsilonname[inds__?UpIndexQ]] :=
Module[{dummy = DummyIn[vbundle],
first = First[{inds}]}, -$ExtrinsicKOnSSSign extrinsicKname[dummy, -dummy] epsilonname[inds]
+ $AccelerationOfnSign norm accelerationname[-dummy] superepsilonname[dummy, inds]];
Protect[Evaluate[prot]]];];
(*Gauss Codazzi Rules,
for abstract indices.Only for Riemann.Norms are wrong*)
With[{riemann = Riemann[covd], superRiemann = Riemann[superCD],
superRicci = Ricci[superCD],
superRicciScalar = RicciScalar[superCD], K = extrinsicKname,
AA = accelerationname},
xAct`xTensor`Private`GaussCodazziRules[metric] :=
{superRiemann[a_?AIndexQ, b_?AIndexQ, c_?AIndexQ,
d_?AIndexQ] :>
Module[{e = DummyIn@vbundle, PDK},
PDK[x_, y_] := projectorname[vector[e] superCD[-e]@K[x, y]];
riemann[a, b, c,
d] + $RiemannSign (-K[a, c] K[b, d]/norm + K[a, d] K[b, c]/norm -
AA[b] AA[d] vector[a] vector[c]/norm - K[b, e] K[d, -e] vector[a] vector[c]/norm^2 +
AA[a] AA[d] vector[b] vector[c]/norm + K[a, e] K[d, -e] vector[b] vector[c]/norm^2 +
AA[b] AA[c] vector[a] vector[d]/norm + K[b, e] K[c, -e] vector[a] vector[d]/norm^2 -
AA[a] AA[c] vector[b] vector[d]/norm - K[a, e] K[c, -e] vector[b] vector[d]/
norm^2 + $ExtrinsicKOnSSSign (vector[b] vector[c] PDK[a, d]/norm^2 +
vector[a] vector[d] PDK[b, c]/norm^2 - vector[a] vector[c] PDK[b, d]/norm^2 -
vector[b] vector[d] PDK[a, c]/norm^2 + vector[d] covd[a]@K[b, c]/norm -
vector[c] covd[a]@K[b, d]/norm - vector[d] covd[b]@K[a, c]/norm +
vector[c] covd[b]@K[a, d]/norm + vector[b] covd[c]@K[a, d]/norm -
vector[a] covd[c]@K[b, d]/norm - vector[b] covd[d]@K[a, c]/norm +
vector[a] covd[d]@K[b, c]/norm) + $AccelerationOfnSign (vector[b] vector[
d] covd[c]@AA[a]/norm - vector[a] vector[d] covd[c]@AA[b]/norm -
vector[b] vector[c] covd[d]@AA[a]/norm + vector[a] vector[c] covd[d]@AA[b]/norm))],
superRicci[a_?AIndexQ, b_?AIndexQ] :>
Module[{c = DummyIn@vbundle}, ReleaseHold[Hold[superRiemann[a, -c, b, c]] /.xAct`xTensor`Private`GaussCodazziRules[metric]]],
superRicciScalar[] :>Module[{a = DummyIn@vbundle, b = DummyIn@vbundle},
Expand[(metric[a, b] + vector[a] vector[b]/norm) ReleaseHold[Hold[superRicci[-a, -b]] /. xAct`xTensor`Private`GaussCodazziRules[metric]]]]}];
If[$ProtectNewSymbols,
Protect[extrinsicKname, accelerationname, projectorname]];]];
DirectionVectorQ[expr_]:=False;
SetSlicingUpToScreenSpaceObinna[g_?MetricQ, u_, normu_: - 1, h_,
cd_, {cdpost_String, cdpre_String}, n_, normn_: 1, NSS_,
cd2_, {cd2post_String, cd2pre_String}, SpaceTimeType_?SpaceTimeQ] :=
Module[{m, p, q, DummyS, DummyV, DummyT, ui, indsdimminustwo,
indsdim, dim, prot,prot2},
With[{Manifold = ManifoldOfCovD@CovDOfMetric[g],
CD = CovDOfMetric[g],ah=a[h],Hh=H[h]},
dim = DimOfManifold[Manifold];
With[{ind1 = DummyIn[Tangent[Manifold]],ind2 = DummyIn[Tangent[Manifold]],ind3 = DummyIn[Tangent[Manifold]],
ind4 = DummyIn[Tangent[Manifold]], ind5 = DummyIn[Tangent[Manifold]], ind6 = DummyIn[Tangent[Manifold]], ind7 = DummyIn[Tangent[Manifold]],
i1 = DummyIn[Tangent[Manifold]], i2 = DummyIn[Tangent[Manifold]], i3 = DummyIn[Tangent[Manifold]], i4 = DummyIn[Tangent[Manifold]], i5 = DummyIn[Tangent[Manifold]],
dummy = DummyIn[Tangent[Manifold]]},
DefTensor[u[-ind1], {Manifold},PrintAs -> "\!\(" <> ToString[u] <> "\&-\)"]
Off[DefMetric::old];
(* Definition of the space induced metric. Standard in xTensor.*)
DefMetric[1, h[-ind1, -ind2], cd, {cdpost, cdpre},InducedFrom -> {g, u},PrintAs -> "\!\(" <> ToString[h] <> "\&-\)"];
On[DefMetric::old];
(* Definition of the direction vector *)
DefTensor[n[-ind1], {Manifold}, OrthogonalTo -> {u[ind1]},ProjectedWith -> {h[ind1, -ind2]},PrintAs -> "\!\(" <> ToString[n] <> "\&-\)"];
(*A boolean value which is used for furtehr checks and avoids the user to create problems if the screen space splitting had not been defined.*)
DirectionVectorQ[n]^=True;
(* The Screen space metric definition. We call the function defined above which contains most of the definitions.
This is the central part of our splitting here.*)
DefScreenSpaceMetric[NSS[-ind1, -ind2], Manifold,cd2, {cd2post, cd2pre}, {h, n}, SpaceTimeType];
(*We remove automatic Leibniz rule when there is a Scalar Head.This is to ensure that the induced derivative does not spoil an'InducedDecomposition'*)
prot = Unprotect[cd];
prot2 = Unprotect[cd2];
cd[a_][Scalar[expr_]] =.;
cd2[a_][Scalar[expr_]] =.;
Protect[prot];
Protect[prot2];
(* Functions used in the post processing of the SplitPerturbations*)
$Rulecdh[h1_]:={
h1[-a_,b_] cd[a_][expr1_]:>cd[b][expr1],
h1[a_,b_] cd[-a_][expr1_]:>cd[b][expr1],
h1[b_,-a_] cd[a_][expr1_]:>cd[b][expr1],
h1[b_,a_] cd[-a_][expr1_]:>cd[b][expr1],
h1[-a_,b_] cd[c_]@cd[a_][expr1_]:>cd[c]@cd[b][expr1],h1[a_,b_] cd[c_]@cd[-a_][expr1_]:>cd[c]@cd[b][expr1],
h1[b_,-a_] cd[c_]@cd[a_][expr1_]:>cd[c]@cd[b][expr1],h1[b_,a_] cd[c_]@cd[-a_][expr1_]:>cd[c]@cd[b][expr1]};
$RulecdNSS[NSS1_]:={NSS1[-a_,b_] cd2[a_][expr1_]:>cd2[b][expr1],
NSS1[a_,b_] cd2[-a_][expr1_]:>cd2[b][expr1],
NSS1[b_,-a_] cd2[a_][expr1_]:>cd2[b][expr1],
NSS1[b_,a_] cd2[-a_][expr1_]:>cd2[b][expr1],
NSS1[-a_,b_] cd2[c_]@cd2[a_][expr1_]:>cd2[c]@cd2[b][expr1],
NSS1[a_,b_] cd2[c_]@cd2[-a_][expr1_]:>cd2[c]@cd2[b][expr1],
NSS1[b_,-a_] cd2[c_]@cd2[a_][expr1_]:>cd2[c]@cd2[b][expr1],
NSS1[b_,a_] cd2[c_]@cd2[-a_][expr1_]:>cd2[c]@cd2[b][expr1]};
(*The default positionof indices for the extrinsic curvature and the acceleration is down*)
(SlotsOfTensor[#] ^:= {-Tangent[Manifold], -Tangent[Manifold]}) & /@ {ExtrinsicK[h],ExtrinsicK[NSS]};
(SlotsOfTensor[#] ^:= {-Tangent[Manifold]}) & /@ {Acceleration[u],Acceleration[n]};
(*The acceleration of ushould vanish for homogeneous spacetimes.*)
Acceleration[u][ind1_] = 0;
(* And it should be of unit norm or at least the norm specified by the user *)
AutomaticRules[u, MakeRule[{u[ind1] u[-ind1], normu}]];
AutomaticRules[u, MakeRule[{u[-ind1] g[ind1, ind2], u[ind2]}]];
(* And similarly for n where the norm should also be unity or the one specified.*)
AutomaticRules[n, MakeRule[{n[ind1] n[-ind1], normn}]];
AutomaticRules[n, MakeRule[{n[-ind1] g[ind1, ind2], n[ind2]}]];
(* The acceleration of the background n should vanish as well*)
Acceleration[n][ind1_] = 0;
(*Rules for the antisymmetric tensor and its projected version *)
(* Currently not implemented on the sphere, but we should ! TODO*)
If[IntegerQ@dim && dim >= 2,
indsdim = GetIndicesOfVBundle[Tangent@Manifold, dim, {ind5}];
AutomaticRules[epsilon[g], MakeRule[
Evaluate[{epsilon[g]@@indsdim u[-indsdim[[1]]] h[-indsdim[[2]], ind5],
ReplaceIndex[Evaluate[epsilon[g] @@ indsdim],indsdim[[2]] -> ind5] u[-indsdim[[1]]]}]
]];];
(*We specify the spacetime. Here it is stollen from basic 1+3*)
(DefScreenProjectedTensor[#[[1]],NSS,SpaceTimesOfDefinition->{"Background"},TensorProperties->{"Traceless","Transverse","SymmetricTensor"},PrintAs->#[[2]]])&/@$ListFieldsBackgroundOnly[h,NSS];
ah::usage=a::usage;
Hh::usage=H::usage;
(* Patch to have a nice output for the Hubble factor *)
Hh/:PrintAs[Hh]=.;
PrintAs[Hh]^:=If[$ConformalTime,"\[ScriptCapitalH]","H"];
(* TODO find a way to remove all these useless LI indices*)
Unprotect[NoScalar];
(* We know it is safe to remove the scalar heads on scale factors and Hubble factors, because we shall never replace a Hubble factor nor a scale factor by something else which is not a pure scalar field...*)
NoScalar[Power[Scalar[Hh[LI[0],LI[0],LI[0]]],ni_Integer]]:=Power[Hh[(*LI[0],LI[0],LI[0]*)],ni];
NoScalar[Power[Scalar[ah[LI[0],LI[0],LI[0]]],ni_Integer]]:=Power[ah[(*LI[0],LI[0],LI[0]*)],ni];
Protect[NoScalar];
ah/:ah[LI[0],LI[0],LI[r_?((IntegerQ@#)&&(#>0)&)]]:=0;
Hh/:Hh[LI[0],LI[q_],LI[LI[r_?((IntegerQ@#)&&(#>0)&)]]]:=0;
If[SpaceTimeType==="Minkowski",
ah[LI[0],LI[0],LI[0]]=1;
Hh[LI[0],LI[0],LI[0]]=0;,
ah[LI[0],LI[1],LI[0]]:=Hh[LI[0],LI[0],LI[0]]ah[LI[0],LI[0],LI[0]];
ah[LI[0],LI[q_?((IntegerQ[#]&&#>=2)&)],LI[0]]:=NoScalar@org@Nest[LieD[u[ind1]][#]&,ah[LI[0],LI[0],LI[0]],q];
DefConformalMetric[g,ah];
];
(* Obvious. Should be automatic in xPert*)
cd[_][$PerturbationParameter]=0;
cd2[_][$PerturbationParameter]=0;
(* The vector used for the background slicing is indeed background so should not be perturbed *)
(* However we want to be able to perturb the normal vector to constan time hypersurfaces.*)
(* In order to do so, we will use the vector \[ScriptCapitalN][h][\[Mu]] which is equal to the background normal vector n^\[Mu].
Howeverm it is ano rule attached to it, so it can be perturbed. Its perturbation is then d\[ScriptCapitalN][LI[p],\[Mu]] *)
(* See the section devoted to the perturbation of the normal vector for more explanations*)
Perturbation[u[ind1_],ni_]^:=0/;ni>=1;
(* A rule to force the natural appearance of the extrinsic curvature of the direction vector*)
CD[ind1_][n[ind2_]]:=cd[ind1][n[ind2]];
(* Lie Derivatives of metric with up indices should be automatic*)
(*AutomaticRules[g,BuildRule[Evaluate[{LieD[u[dummy]][g[ind1,ind2]],LieD[u[dummy]][g[ind1,ind2]]//MetricToProjector//ToCanonical}],MetricOn->None]];*)
AutomaticRules[g,BuildRule[Evaluate[{LieD[u[dummy]][g[-ind1,-ind2]],LieD[u[dummy]][g[-ind1,-ind2]]//MetricToProjector//ToCanonical}],MetricOn->None]];
(*AutomaticRules[g,BuildRule[Evaluate[{LieD[n[dummy]][g[ind1,ind2]],LieDToCovD[LieD[n[dummy]][g[ind1,ind2]],CD]//GradNormalToExtrinsicK//ToCanonical}],MetricOn->None]];*)
AutomaticRules[g,BuildRule[Evaluate[{LieD[n[dummy]][g[-ind1,-ind2]],LieDToCovD[LieD[n[dummy]][g[-ind1,-ind2]],CD]//GradNormalToExtrinsicK//ToCanonical}],MetricOn->None]];
Which[
FlatSpaceBool[SpaceTimeType],
Riemann[cd][i1_,i2_,i3_,i4_]=0;
Ricci[cd][i1_,i2_]=0;
RicciScalar[cd][]=0;,
CurvedSpaceBool[SpaceTimeType],
DefTensor[\[ScriptK][h][],{Manifold},PrintAs->StringJoin["\[ScriptK]"(*,ToString[h]*)]];
CD[ind1_][\[ScriptK][h][]]:=0(*-1*(-2/3)*u[ind1]\[ScriptK][h][]CD[-ind2][u[ind2]]*); (* This is because Subscript[\[ScriptCapitalL], u](\[ScriptK]^(1/2) Subscript[\[Epsilon]^\[Mu], \[Nu]\[Sigma]]) must be zero since Subscript[\[ScriptCapitalL], u](Subscript[C^\[Mu], \[Nu]\[Sigma]])=0*)
(* Given that we first work on the conformal metric where the trace of the extrinsic curvature is 0, then this point does not matter at all for us.*)
cd[ind1_][\[ScriptK][h][]]=0;
If[SpaceTimeType==="FLCurved",
IndexSet[Riemann[cd][i1_,i2_,i3_,i4_],\[ScriptK][h][](h[i1,i3]h[i2,i4]-h[i1,i4]h[i2,i3])];
IndexSet[Ricci[cd][i1_,i2_],\[ScriptK][h][]h[i1,i2](h[i3,-i3]-1)];