(*

BE CAREFUL - ϕ is not SQRT[ξ ξ]!!!

*)

OUTPUT=0;

LichRUN=0;

SetDirectory["~/Work/SpecTest4/"];

data={
{"test_even_04_44"},
{"test"},
{"test_05R_20l"},
{"test_05R_24l"},
{"test_05R_30l"},
{"test_05R_34l"},
{"test_05R_40l"},
{"test_05R_44l"},
{"test_05R_50l"}
};

item=2;

Print[data[[item,1]]];

SetDirectory[data[[item,1]]];

(*  1 - if even,   2 - if not *)

EvenR=1;

EvenA=1;

Coordinates

file1=OpenRead["patch.dat",BinaryFormat->True];

{NR,NA}=BinaryReadList[file1,"Integer32",2]

tmplist=Partition[BinaryReadList[file1,{"Real64","Real64","Integer32","Integer32","Integer32","Integer32","Integer32"}],NA];

Close[file1];

rarr=tmplist[[All,1,1]]

aarr=tmplist[[1,All,2]]

rpatch=tmplist[[All,All,1]];

apatch=tmplist[[All,All,2]];

Read data

garr=Import["data.out","Real64"];

garr=Partition[Partition[Partition[garr,NA],NR],5];

Length[garr]

ϕ2arr=Import["phi.out","Real64"];

ϕ2arr=Partition[Partition[ϕ2arr,NA],NR];

Ricciarr=Import["Ricci.out","Real64"];

Ricciarr=Partition[Partition[Ricciarr,NA],NR];

R2sqrarr=Import["R2sqr.out","Real64"];

R2sqrarr=Partition[Partition[R2sqrarr,NA],NR];

R4sqrarr=Import["R4sqr.out","Real64"];

R4sqrarr=Partition[Partition[R4sqrarr,NA],NR];

If[LichRUN==1,

LichRow=Import["LichRow.out","Integer64"];

LichCol=Import["LichCol.out","Integer64"];

LichVal=Import["LichVal.out","Real64"];

Lich=SparseArray[Thread[Transpose[{LichRow+1,LichCol+1}]->LichVal]];

maxeval=Eigenvalues[Lich,1][[1]];

mineval=Eigenvalues[Lich,-1][[1]];

condition=Abs[maxeval/mineval];

];

Select

tval=-1;

gdat=garr[[tval]];

sqrϕ=Sqrt[ϕ2arr[[tval]]];

ϕ=ϕ2arr[[tval]];

Tdat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[gdat[[1]]]}];

Adat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[gdat[[2]]]}];

Bdat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[gdat[[3]]]}];

Fdat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[gdat[[4]]]}];

Sdat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[gdat[[5]]]}];

ϕdat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[ϕ]}];

Riccidat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[Ricciarr[[tval]]+20]}];

R2sqrdat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[R2sqrarr[[tval]]-80]}];

R4sqrdat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[R4sqrarr[[tval]]]}];

ϕdatCUT=Transpose[{Flatten[rpatch[[EvenR;;NR-1,EvenA;;NA]]],Flatten[apatch[[EvenR;;NR-1,EvenA;;NA]]],Flatten[ϕ[[EvenR;;NR-1,EvenA;;NA]]]}];

RiccidatCUT=Transpose[{Flatten[rpatch[[EvenR;;NR-1,EvenA;;NA-1]]],Flatten[apatch[[EvenR;;NR-1,EvenA;;NA-1]]],Flatten[Ricciarr[[tval]][[EvenR;;NR-1,EvenA;;NA-1]]+20]}];

R2sqrdatCUT=Transpose[{Flatten[rpatch[[EvenR;;NR-1,EvenA;;NA-1]]],Flatten[apatch[[EvenR;;NR-1,EvenA;;NA-1]]],Flatten[R2sqrarr[[tval]][[EvenR;;NR-1,EvenA;;NA-1]]-80]}];

R4sqrdatCUT=Transpose[{Flatten[rpatch[[EvenR;;NR-1,EvenA;;NA-1]]],Flatten[apatch[[EvenR;;NR-1,EvenA;;NA-1]]],Flatten[R4sqrarr[[tval]][[EvenR;;NR-1,EvenA;;NA-1]]]}];

A2dat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[rpatch^2*(1-rpatch^2)*gdat[[2]]]}];

B2dat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[apatch^2*gdat[[3]]]}];

F2dat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[rpatch*apatch*gdat[[4]]]}];

S2dat=Transpose[{Flatten[rpatch],Flatten[apatch],Flatten[apatch*gdat[[5]]]}];

Display - metric

plt1=ListPlot3D[R4sqrdat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[R4sqrdat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[Tdat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[Tdat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[Sdat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[Sdat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[Bdat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[Bdat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[Adat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[Adat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[A2dat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[A2dat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[B2dat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[B2dat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[F2dat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[F2dat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[S2dat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[S2dat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

Display - error  no cut

plt1=ListPlot3D[ϕdat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[ϕdat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[Riccidat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[Riccidat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[R2sqrdat,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[R2sqrdat,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

Display - error  CUT

plt1=ListPlot[Transpose[{rarr,ϕ[[All,NA]]}]];

plt2=Plot[1.5*10^-9(1-r^2)^(2),{r,0,1}];

Show[plt1,plt2,PlotRange->All]

plt1=ListPlot3D[ϕdatCUT,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[ϕdatCUT,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[RiccidatCUT,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[RiccidatCUT,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]

plt1=ListPlot3D[R2sqrdatCUT,PlotRange->All,AxesLabel->{"r","a"},Mesh->False];

plt2=ListPointPlot3D[R2sqrdatCUT,PlotRange->All,AxesLabel->{"r","a"}];

Show[plt1,plt2]