(* ::Package:: *)

<< SimulationTools`

numPoints = 200;

ParallelNeeds["SimulationTools`"];

SpinWeightedSphericalHarmonic[s_, l_, m_, th_, 
   ph_] := (Factorial[l + m] Factorial[
       l - m] (2 l + 1) Factorial[l + s]^(-1) Factorial[
        l - s]^(-1) (4 Pi)^(-1))^(1/2) (Sin[th/2])^(2 l) Sum[
    Binomial[l - s, r] Binomial[l + s, 
      r + s - m] (-1)^(l - r - s) Exp[
      I m ph] Cot[th/2]^(2 r + s - m), {r, Max[m - s, 0], 
     Min[l - s, l + m]}];

With[{s = -2}, 
  Table[SpinWeightedSphericalHarmonic[s, l, m, (0 | 0.), \[Phi]_] := 
    Evaluate[
     Limit[SpinWeightedSphericalHarmonic[s, l, 
       m, \[Theta], \[Phi]], \[Theta] -> 0]];
   SpinWeightedSphericalHarmonic[s, l, 
     m, (\[Pi] | N[\[Pi]]), \[Phi]_] := 
    Evaluate[
     Limit[SpinWeightedSphericalHarmonic[s, l, 
       m, \[Theta], \[Phi]], \[Theta] -> \[Pi]]], {l, Abs[s], 
    8}, {m, -l, l}]];

$SimulationPath = {"/scratch/barry/Simulations/LIGO"};

ParallelEvaluate[$SimulationPath = \
{"/scratch/barry/Simulations/LIGO"}];

sim = "LIGO_2_28";

M = ReadADMMass[sim];

coords1min = ReadAHColumns[sim, 1, {15, 17, 19}];
coords2min = ReadAHColumns[sim, 2, {15, 17, 19}];
coords3min = ReadAHColumns[sim, 3, {15, 17, 19}];
coords1max = ReadAHColumns[sim, 1, {16, 18, 20}];
coords2max = ReadAHColumns[sim, 2, {16, 18, 20}];
coords3max = ReadAHColumns[sim, 3, {16, 18, 20}];

mergertime = CoordinateRanges[coords3max][[1, 1]];

Table[\[Psi]4[l, m] = 
    Interpolation[
     PadRight[
      Shifted[500 ReadPsi4[sim, l, m, 500], -RadialToTortoise[500, 
         M]], 4000]];, {l, 2, 8}, {m, -l, l}];


(*Clear[h]*)


(*hlmax4=Sum[h[l,m][t-rstar]Expand[SpinWeightedSphericalHarmonic[-2,l,\
m,\[Theta],\[Phi]]],{l,Abs[2],4},{m,-l,l}];*)


Table[h[l, m] = 
   Interpolation[
    PadRight[
     Shifted[500 Slab[
        StrainFromPsi4[ReadPsi4[sim, l, m, 500], 0.005], ;; 
         1560], -RadialToTortoise[500, M]], 4000]], {l, 2, 8}, {m, -l,
    l}];


(*Manipulate[With[{x0=0.06,x1=0.09,x2=0.0525,x3=0.0825},Plot[Re[h[2,2]\
[t-RadialToTortoise[r,M]]SpinWeightedSphericalHarmonic[-2,2,2,\[Theta]\
,\[Phi]]+h[2,-2][t-RadialToTortoise[r,M]]\
SpinWeightedSphericalHarmonic[-2,2,-2,\[Theta],\[Phi]]],{r,0,200},\
ColorFunctionScaling\[Rule]False,ColorFunction\[Rule]Function[{r,f},\
If[f>0,ColorData[{"DeepSeaColors","Reverse"}][(f-x0)/(x1-x0)],\
ColorData[{"SolarColors","Reverse"}][(-f-x0)/(x1-x0)]]],Background\
\[Rule]Black,PlotRange\[Rule]All,PlotLabel\[Rule]t]],{\[Theta],0,\[Pi]\
},{\[Phi],0,2\[Pi]},{t,0,1500}]*)


strainPlot[t_] := 
 Module[{coordmin1, coordmin2, coordmin3, coordmax1, coordmax2, 
   coordmax3, hors, hor1, hor2, x, y, z}, 
  With[{x0 = 1.2 0.06, x1 = 1.2 0.0825, x2 = 1.2 0.0525, 
    x3 = 1.2 0.0825}, 
   DensityPlot3D[
    Evaluate[
     With[{r = Sqrt[x^2 + y^2 + z^2], \[Theta] = 
        ArcCos[z/Sqrt[x^2 + y^2 + z^2]], \[Phi] = ArcTan[x, y], 
       rstar = RadialToTortoise[Sqrt[x^2 + y^2 + z^2], M]}, 
      Erf[Abs[r]/(Sqrt[2] 20)] Max[
        Erf[(TortoiseToRadial[t - 20, M] - r)/(Sqrt[2] 20)], 0.] Max[
        Erf[(210 - Sqrt[x^2 + y^2])/(10 Sqrt[2])], 0] Max[
        Erf[(210 - Abs[z])/(10 Sqrt[2])], 0] Re[
        h[2, 2][t - rstar] 1/2 E^(2 I \[Phi]) Sqrt[
           5/\[Pi]] Cos[1/2 \[Theta]]^4 + 
         h[2, -2][t - rstar] 1/2 E^(-2 I \[Phi]) Sqrt[
           5/\[Pi]] Sin[1/2 \[Theta]]^4]]], {x, -200, 200}, {y, -200, 
     200}, {z, -200, 200}, ColorFunctionScaling -> False, 
    ColorFunction -> 
     Function[{f}, 
      If[f > 0, 
       ColorData[{"DeepSeaColors", "Reverse"}][(f - x0)/(x1 - x0)], 
       ColorData[{"SolarColors", "Reverse"}][(-f - x0)/(x1 - x0)]]], 
    OpacityFunctionScaling -> False, 
    OpacityFunction -> 
     Function[{f}, Max[Min[(Abs[f] - x2)/(x3 - x2), 1.0], 0]], 
    Background -> Black, Boxed -> False, Axes -> False, 
    ViewPoint -> Front, ViewAngle -> 31 Degree, 
    PlotPoints -> numPoints, 
    PlotRange -> 200 {{-1, 1}, {-1, 1}, {-1, 1}}]]]

its = ReadIterations[sim, "Psi4r", "xyz", Map -> 0, 
   RefinementLevel -> 3];


dt = ReadTimeStep[sim]/2^7;


contourPlotRe[it_] := 
 Module[{coordmin1, coordmin2, coordmin3, coordmax1, coordmax2, 
   coordmax3, psi4r}, 
  psi4r = ReadGridFunction[sim, "Psi4r", "xyz", Map -> 0, 
    Iteration -> it, RefinementLevel -> 3]; 
  ListContourPlot3D[psi4r, 
   Contours -> {-0.01, -0.006, -0.003, -0.001, 0.001, 0.003, 0.006, 
     0.01}, ContourStyle -> Opacity[0.3], Boxed -> False, 
   Background -> Black, Mesh -> None, Axes -> False, 
   ColorFunction -> 
    Function[{x, y, z, f}, 
     If[f > 0, ColorData[{"DeepSeaColors", "Reverse"}][100 f], 
      ColorData[{"SolarColors", "Reverse"}][-100 f]]], 
   ColorFunctionScaling -> False, 
   PlotRange -> 200 {{-1, 1}, {-1, 1}, {-1, 1}}]]


horizonPlot[time_] := 
  Module[{coordmin1, coordmin2, coordmin3, coordmax1, coordmax2, 
    coordmax3, hor1, hor2, hors}, 
   If[time >= mergertime, 
    coordmin3 = Interpolation[coords3min, time];
    coordmax3 = Interpolation[coords3max, time];
    hors = 
     Graphics3D[{Black, Opacity[1.0], Specularity[0.3], 
       Ellipsoid[(coordmin3 + coordmax3)/2, (coordmax3 - coordmin3)/
         2]}];, coordmin1 = Interpolation[coords1min, time];
    coordmin2 = Interpolation[coords2min, time];
    coordmax1 = Interpolation[coords1max, time];
    coordmax2 = Interpolation[coords2max, time];
    hor1 = 
     Graphics3D[{Black, Opacity[1.0], Specularity[0.3], 
       Ellipsoid[(coordmin1 + coordmax1)/2, (coordmax1 - coordmin1)/
         2]}];
    hor2 = 
     Graphics3D[{Black, Opacity[1.0], Specularity[0.3], 
       Ellipsoid[(coordmin2 + coordmax2)/2, (coordmax2 - coordmin2)/
         2]}];
    hors = Sequence[hor1, hor2];];
   Show[hors, Boxed -> False, Background -> Black, Mesh -> None, 
    Axes -> False]];


Quiet[\[CapitalOmega] = 
   Interpolation[NDerivative[1][Slab[ReadBinaryPhase[sim], 0 ;; 800]],
     110.1362637364`]];


contourPlotRe201 = 
  Show[contourPlotRe[its[[201]]], horizonPlot[dt its[[201]]]];


(*strainAndContourPlotRe201=Block[{numPoints=400},Show[strainPlot[dt \
its[[201]]],contourPlotRe[its[[201]]],horizonPlot[dt its[[201]]]]];*)


fullPlot[it_, dx_: 0] := Module[{t, \[Theta], plot}, t = it dt;
   Which[it < its[[1]], 
    plot = Show[contourPlotRe201, ViewAngle -> 1.2 Degree, 
       ViewPoint -> {dx Cos[\[CapitalOmega] t] - 
          2 Sin[\[CapitalOmega] t], -2 Cos[\[CapitalOmega] t] - 
          dx Sin[\[CapitalOmega] t], 1.5}];, 
    its[[1]] <= it <= its[[201]], 
    plot = Show[contourPlotRe201, 
       ViewAngle -> 
        Min[Max[(32 (-it + its[[1 + 24]]) + 
              1.2 (it - its[[201 - 24]]))/(its[[1 + 24]] - 
              its[[201 - 24]]), 1.2], 32] Degree(*(32 (-it+
       its[[1]])+1.2 (it-its[[201]]))/(its[[1]]-its[[201]])*), 
       ViewPoint -> {dx, -2, 
         Min[Max[(0 (-it + its[[1 + 24]]) + 
              1.5 (it - its[[201 - 24]]))/(its[[1 + 24]] - 
              its[[201 - 24]]), 0], 1.5](*(0 (-it+its[[1]])+(it-
         its[[201]]) 1.5)/(its[[1]]-its[[201]])*)}];, 
    it > its[[201]], \[Theta] = 
     Min[Max[(0 (-it + its[[201 + 4 24]]) + \[Pi]/
            2 (it - its[[2269 - 4 24]]))/(its[[201 + 4 24]] - 
          its[[2269 - 4 24]]), 0.9813290485093916`], \[Pi]/2];
    plot = 
     Show[strainPlot[t], contourPlotRe[it], horizonPlot[t], 
      ViewAngle -> 32 Degree, 
      ViewPoint -> {dx, -2 Sin[\[Theta]], 2 Cos[\[Theta]]}];];
   plot];


fullPlot[originalPlot_, it_, dx_] := 
  Module[{t, \[Theta], plot}, t = it dt;
   Which[it < its[[1]], 
    plot = Show[originalPlot, ViewAngle -> 1.2 Degree, 
       ViewPoint -> {dx Cos[\[CapitalOmega] t] - 
          2 Sin[\[CapitalOmega] t], -2 Cos[\[CapitalOmega] t] - 
          dx Sin[\[CapitalOmega] t], 1.5}];, 
    its[[1]] <= it <= its[[201]], 
    plot = Show[originalPlot, 
       ViewAngle -> 
        Min[Max[(32 (-it + its[[1 + 24]]) + 
              1.2 (it - its[[201 - 24]]))/(its[[1 + 24]] - 
              its[[201 - 24]]), 1.2], 32] Degree(*(32 (-it+
       its[[1]])+1.2 (it-its[[201]]))/(its[[1]]-its[[201]])*), 
       ViewPoint -> {dx, -2, 
         Min[Max[(0 (-it + its[[1 + 24]]) + 
              1.5 (it - its[[201 - 24]]))/(its[[1 + 24]] - 
              its[[201 - 24]]), 0], 1.5](*(0 (-it+its[[1]])+(it-
         its[[201]]) 1.5)/(its[[1]]-its[[201]])*)}];, 
    it > its[[201]], \[Theta] = 
     Min[Max[(0 (-it + its[[201 + 4 24]]) + \[Pi]/
            2 (it - its[[2269 - 4 24]]))/(its[[201 + 4 24]] - 
          its[[2269 - 4 24]]), 0.9813290485093916`], \[Pi]/2];
    plot = 
     Show[originalPlot, ViewAngle -> 32 Degree, 
      ViewPoint -> {dx, -2 Sin[\[Theta]], 2 Cos[\[Theta]]}];];
   plot];


tOfIt[it_ /; it > its[[201]]] := dt it;
tOfIt[it_ /; its[[1]] <= it <= its[[201]]] := dt its[[201]];
tOfIt[it_ /; it < its[[1]]] := dt (it + its[[201]]);


it0 = -its[[1601]];
itf = its[[2269]];
dit = 512;


fullPlotWithLabel[it_] := 
  Module[{plot, image, left, right, imageleft, imageright}, 
   plot = fullPlot[it];
   left = fullPlot[plot, it, -0.2];
   right = fullPlot[plot, it, 0.2];
   image = 
    CriticalSection[{rasterlock}, 
     Show[Rasterize[plot, ImageSize -> {1080, 1080}], 
      Graphics[{White, 
        Style[Text[
          
          " " <> ToString[
            NumberForm[(tOfIt[it] - tOfIt[it0]) 1000 M (29 + 
                36) SolarMassInSeconds, {Infinity, 1}]] <> " ms ", 
          Scaled[{0.5, 0.98}], 
          Background -> Directive[Black, Opacity[0.3]]], 20]}]]];
   imageleft = 
    CriticalSection[{rasterlock}, 
     Show[Rasterize[left, ImageSize -> {1080, 1080}], 
      Graphics[{White, 
        Style[Text[
          " " <> ToString[
            NumberForm[(tOfIt[it] - tOfIt[it0]) 1000 M (29 + 
                36) SolarMassInSeconds, {Infinity, 1}]] <> " ms ", 
          Scaled[{0.5, 0.98}], 
          Background -> Directive[Black, Opacity[0.3]]], 20]}]]];
   imageright = 
    CriticalSection[{rasterlock}, 
     Show[Rasterize[right, ImageSize -> {1080, 1080}], 
      Graphics[{White, 
        Style[Text[
          " " <> ToString[
            NumberForm[(tOfIt[it] - tOfIt[it0]) 1000 M (29 + 
                36) SolarMassInSeconds, {Infinity, 1}]] <> " ms ", 
          Scaled[{0.5, 0.98}], 
          Background -> Directive[Black, Opacity[0.3]]], 20]}]]];
   {imageleft, image, imageright}];


CreateDirectory["FullMovie-centre"];


CreateDirectory["FullMovie-left"];


CreateDirectory["FullMovie-right"];


DistributeDefinitions["Global`*"];


ParallelDo[Module[{plots}, Print[{(it - it0)/dit, (itf - it0)/dit}];
   plots = fullPlotWithLabel[it];
   MapThread[
    If[! FileExistsQ[#1 <> "/FullMovie" <> 
         IntegerString[(it - it0)/dit, 10, 4] <> ".png"], 
      Export[#1 <> "/FullMovie" <> 
        IntegerString[(it - it0)/dit, 10, 4] <> 
        ".png", #2]] &, {{"FullMovie-left", "FullMovie-centre", 
      "FullMovie-right"}, plots}]];, {it, it0 + 1300 dit, itf, dit}, 
 Method -> "FinestGrained"]
