User: barry.wardell Date: 2016/06/05 05:07 PM
Added: /about/gallery/gw150914/ Curvature.m, Strain.m
Log: GW150914 gallery: Add scripts for generating images with SimulationTools
File Changes:
Directory: /about/gallery/gw150914/ ===================================
File [added]: Curvature.m Delta lines: +71 -0 =================================================================== --- about/gallery/gw150914/Curvature.m (rev 0) +++ about/gallery/gw150914/Curvature.m 2016-06-05 22:07:20 UTC (rev 1545) @@ -0,0 +1,71 @@ +<<SimulationTools` + + +simDirs = ReplaceList[$ScriptCommandLine, + {___, "--simulations-directory", val_, ___} :> val]; + +$SimulationPath = Join[simDirs, $SimulationPath]; + +outDir = Directory[]; + +sim="GW150914_32"; + + +coords1min=ReadAHColumns[sim,1,{15,17,19}]; +coords2min=ReadAHColumns[sim,2,{15,17,19}]; +coords1max=ReadAHColumns[sim,1,{16,18,20}]; +coords2max=ReadAHColumns[sim,2,{16,18,20}]; +{coords1mini,coords2mini,coords1maxi,coords2maxi}=Interpolation/@{coords1min,coords2min,coords1max,coords2max}; +horizons[t_]:=Show[Graphics3D[{Black,Opacity[1.0],Specularity[0.7],Ellipsoid[(coords1mini[t]+coords1maxi[t])/2,(coords1maxi[t]-coords1mini[t])/2]}],Graphics3D[{Black,Opacity[1.0],Specularity[0.7],Ellipsoid[(coords2mini[t]+coords2maxi[t])/2,(coords2maxi[t]-coords2mini[t])/2]}]] + + +(* This function takes a DataRegion on the upper half plane (z[GreaterEqual]0) and reproduces the data on the full domain assuming the quantity is reflection symmetric about the plane z=0 *) +UpperHalfToFull[d_DataRegion, parity_:1]:= +Module[{dx,dy,dz,xrange,yrange,zrange,fulldata}, +{dx,dy,dz}=CoordinateSpacings[d]; +{xrange,yrange,zrange}=CoordinateRanges[d]; +fulldata=Join[ +parity Map[Reverse,ToListOfData[Slab[d,All,All,0;;]],{2}], +ToListOfData[Slab[d,All,All,(0+dz);;]], +3]; +ToDataRegion[fulldata,{xrange[[1]],yrange[[1]],-zrange[[2]]},{dx,dy,dz},"VariableName"->VariableName[d]] +]; + + +Re[Psi]4=ReadGridFunction[sim,"Psi4r","xyz","RefinementLevel"->3,"Map"->0,Iteration->0]; + + +Re[Psi]4full=UpperHalfToFull[Re[Psi]4]; + + +its=ReadIterations[sim,"curvIr","xyz",RefinementLevel->3,Map->0]; + + +Re[ScriptCapitalI]=ReadGridFunction[sim,"curvIr","xyz","RefinementLevel"->3,"Map"->0,"Iteration"->0]; +Im[ScriptCapitalI]=ReadGridFunction[sim,"curvIi","xyz","RefinementLevel"->3,"Map"->0,"Iteration"->0]; +Re[ScriptCapitalJ]=ReadGridFunction[sim,"curvJr","xyz","RefinementLevel"->3,"Map"->0,"Iteration"->0]; +Im[ScriptCapitalJ]=ReadGridFunction[sim,"curvJi","xyz","RefinementLevel"->3,"Map"->0,"Iteration"->0]; + + +Re[ScriptCapitalI]full=UpperHalfToFull[Re[ScriptCapitalI]]; +Im[ScriptCapitalI]full=UpperHalfToFull[Im[ScriptCapitalI],-1]; +Re[ScriptCapitalJ]full=UpperHalfToFull[Re[ScriptCapitalJ]]; +Im[ScriptCapitalJ]full=UpperHalfToFull[Im[ScriptCapitalJ],-1]; + + +Psi4Plot=Show[ListContourPlot3D[Re[Psi]4full,Contours->{-0.01,-0.006,-0.003,-0.001,-0.0005,0.0005,0.001,0.003,0.006,0.01},ContourStyle->Opacity[0.3],Background->Black,Mesh->None,ColorFunction->Function[{x,y,z,f},If[f>0,ColorData[{"DeepSeaColors","Reverse"}][200f],ColorData[{"SolarColors","Reverse"}][-200f]]],ColorFunctionScaling->False,BoxRatios->Automatic],horizons[0],Boxed->False,Axes->False]; + + +Export[FileNameJoin[{outDir,"RePsi4.png"}],Psi4Plot,ImageSize->{Automatic,2160}] + + +ReIPlot=Show[ListContourPlot3D[Re[ScriptCapitalI]full,Contours->10^{-1,-2,-3,-4,-5},ContourStyle->Opacity[0.3],Background->Black,Mesh->None,BoxRatios->Automatic,ColorFunction->ColorData[{"DeepSeaColors","Reverse"}]],horizons[0],Boxed->False,Axes->False]; + + +Export[FileNameJoin[{outDir,"ReI.png"}],ReIPlot,ImageSize->{Automatic,2160}] + + +ImIPlot=Show[ListContourPlot3D[Im[ScriptCapitalI]full,Contours->Join[-10^{-1,-2,-3,-4,-5,-5.5},10^{-1,-2,-3,-4,-5,-5.5}],ContourStyle->Opacity[0.3],Background->Black,Mesh->None,BoxRatios->Automatic,ColorFunction->Function[{x,y,z,f},If[f>0,ColorData[{"DeepSeaColors","Reverse"}][30000f],ColorData[{"SolarColors","Reverse"}][-30000f]]],ColorFunctionScaling->False],horizons[0],Boxed->False,Axes->False]; + + +Export[FileNameJoin[{outDir,"ImI.png"}],ImIPlot,ImageSize->{Automatic,2160}]
File [added]: Strain.m Delta lines: +69 -0 =================================================================== --- about/gallery/gw150914/Strain.m (rev 0) +++ about/gallery/gw150914/Strain.m 2016-06-05 22:07:20 UTC (rev 1545) @@ -0,0 +1,69 @@ +<<SimulationTools` + + +simDirs = ReplaceList[$ScriptCommandLine, + {___, "--simulations-directory", val_, ___} :> val]; + +$SimulationPath = Join[simDirs, $SimulationPath]; + +outDir = Directory[]; + + +sim="GW150914_32"; + + +MADM=ReadADMMass[sim]; + + +coords1min=ReadAHColumns[sim,1,{15,17,19}]; +coords2min=ReadAHColumns[sim,2,{15,17,19}]; +coords1max=ReadAHColumns[sim,1,{16,18,20}]; +coords2max=ReadAHColumns[sim,2,{16,18,20}]; +{coords1mini,coords2mini,coords1maxi,coords2maxi}=Interpolation/@{coords1min,coords2min,coords1max,coords2max}; +horizons[t_]:=Show[Graphics3D[{Black,Opacity[1.0],Specularity[0.7],Ellipsoid[(coords1mini[t]+coords1maxi[t])/2,(coords1maxi[t]-coords1mini[t])/2]}],Graphics3D[{Black,Opacity[1.0],Specularity[0.7],Ellipsoid[(coords2mini[t]+coords2maxi[t])/2,(coords2maxi[t]-coords2mini[t])/2]}]] + + +r0=500; + + +Do[[Psi]4[l,m]=ReadPsi4[sim,l,m,r0],{l,2,8},{m,-l,l}]; +Do[r[Psi]4[l,m]=Shifted[r0 [Psi]4[l,m],-RadialToTortoise[r0,MADM]],{l,2,8},{m,-l,l}]; + + +tMerger=LocateMaximumPoint[Abs[r[Psi]4[2,2]]]; + + +modes=Sort[Flatten[Table[{l,m,Interpolation[Abs[r[Psi]4[l,m]],tMerger]},{l,2,8},{m,-l,l}],1],#1[[3]]>#2[[3]]&]; + + +strongestmodes=modes[[1;;8,1;;2]]; + + +[Omega]0=0.025; + + +Do[rh[l,m]=Psi4ToStrain[r[Psi]4[l,m],[Omega]0],{l,2,8},{m,-l,l}] + + +Table[rhinterp[l,m]=Interpolation[PadRight[Slab[rh[l,m],MinCoordinate[rh[l,m]];;1050],4000]],{l,2,8},{m,-l,l}]; + + +SpinWeightedSphericalHarmonic[s_,l_,m_,[Theta]_,[Phi]_]:=Sqrt[(2 l+1)/(4[Pi]) ((l+m)!(l-m)!)/((l+s)!(l-s)!)] E^(I m [Phi]) Sin[[Theta]/2]^(2 l) !( +*UnderoverscriptBox[([Sum]), (r = Max[m - s, 0]), (Min[l - s, l + m])](Binomial[l - s, r]\ Binomial[l + s, r + s - m]\ +*SuperscriptBox[(((-1))), (l - r - s)]\ \ +*SuperscriptBox[(Cot[[Theta]/2]), (2\ r + s - m)])); + + +rhr[Theta][Phi][r_,[Theta]_,[Phi]_]:=With[{t=700},Re[Sum[If[MemberQ[strongestmodes,{l,m}],rhinterp[l,m][t-TortoiseToRadial[r,MADM]]SpinWeightedSphericalHarmonic[-2,l,m,[Theta],[Phi]],0],{l,2,8},{m,-l,l}]]]; + + +strainPlot=With[{x0=1.15 0.06,x1=1.15 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]},rhr[Theta][Phi][r,[Theta],[Phi]]]], +{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->200,PlotRange->200{{-1,1},{-1,1},{-1,1}}]]; + + +strainHorizonPlot=Show[strainPlot,horizons[745],ViewAngle->28Degree,ImageSize->{400,300}]; + + +Export[FileNameJoin[{outDir,"StrainHorizon.png"}],strainHorizonPlot,ImageSize->{400,300}2160/300]