(* Ground Zero Model; Web: www.911history.de; As of: 2019-04-04; Contents: Creative Commons *)

(* Power Source: Uranium fission; Chain reaction inhibitor: Thorium; Cooling agent: liquid granite (natural convection); Final ignition: pressure-induced phase separation uranium/thorium mix *)

(* shows for WTC2: 'glacier' valley and hot melt-funnel after eruption and a progressive (100 seconds) 125 kt Energy Input, followed by a final 25 kt nuclear fizzle reaction (6 seconds)  *)
(* shows for WTC1: screenshot of the dynamics during the nuclear plasma eruption, acceleration and direction defined by the melt-funnel structure / pit below freight elevator 50 *)
(* shows for WTC7: creation of the eruptive melt-funnel and blocking of the eruption due to ingressing water (cooling of the liquefied granite) *)

(* Ignition (nuclear fizzle) is naturally achieved by: 'void coefficient of reactivity' changes from negative to positive under extreme pressure, as soon as the boiling uranium/thorium mixture has burned it's way deep enough down into the rock, with a heavy liquid granite dome above *)

(* As soon as the boiling stops under pressure, liquid thorium and liquid uranium can separate: https://en.wikipedia.org/wiki/Phase_(matter) *)

(* Uranium: Melting point: 1132 °C | Boiling point: 4131 °C *)
(* Thorium: Melting point: 1750 °C | Boiling point: 4788 °C *)
(* Granite: Melting point: 1250 °C | Boiling point: 2200 - 2900 °C [Silicon dioxide, as one of many components]*)

(* Theory: Silverstein Glacier Valley is the result of a displaced [far reaching] convection flow of liquid granite, a solidified testimony of the energy flow / energy landscape *)
(* Destruction of the Twin Towers is the result of x- and gamma-radiation induced weakening of the steel core, followed by a short-lived nuclear plasma eruption, upshot: 350 m *)

(* Map: overlay of Google maps and elevation data of the glacier valley; source1 (photos): https://www.nytimes.com/2008/09/22/nyregion/22rocks.html *)
(* Map: overlay of Google maps and elevation data of the glacier valley; source2 (elevations): https://dspace.sunyconnect.suny.edu/bitstream/handle/1951/47900/moss-09.pdf *)
(* NOTE: the valley may have indeed existed BEFORE the nuclear event, but was most likely contaminated by it *)
(*  *)




DynamicModule[{red=0,green=0,blue=0,t1=0,t2=0},Column[{Grid[{{Graphics[{Red,Rectangle[]},ImageSize->30],Slider[Dynamic[red]]},{Graphics[{Green,Rectangle[]},ImageSize->30],Slider[Dynamic[green]]},{Graphics[{Blue,Rectangle[]},ImageSize->30],Slider[Dynamic[t1]],"Slow Energy Input WTC2 =",Dynamic[Round[slowEnergyInput2[5400 t1]]],"[kt]","|| Radius Hot Granite Sphere Next Day=",Dynamic[Round[radiusHeatZone2[5400 t1]]],"[m]"},{Graphics[{Black,Rectangle[]},ImageSize->30],Slider[Dynamic[t2]],"Fast Energy Input WTC2=",Dynamic[Round[slowEnergyInput2[5400 t2]]],"[kt]","|| Radius Gaseous Pressure Sphere =",Dynamic[Round[radiusGasousSphere2[5400 t2]]],"[m]"}}],
Style["Resulting radius of a hot sphere:",20,Bold, FontFamily->"Helvetica"],
Style["brown = Energy",16,Brown, FontFamily->"Helvetica"],
Style["blue = size of sphere with 100 °C",16,Blue, FontFamily->"Helvetica"],
Style["orange = size of 100% solid sphere with 1200 °C",16,Orange, FontFamily->"Helvetica"],
Style["purple = size of 100% liquid sphere with 2000 °C",16,Purple, FontFamily->"Helvetica"],
Grid[{{Show[{Plot[{slowEnergyInput2[t1],radiusHeatZone2Temp100[t1],radiusHeatZone2Temp1200[t1],radiusHeatZone2Temp2000[t1]},{t1,0,5400},PlotStyle->{Brown,Blue,Orange,Purple}], Graphics[{PointSize[0.04],RGBColor[(217/255),(177/255),(1/255)],Point[Dynamic[{5400t2,slowEnergyInput2[t2*5400]}]]}]},ImageSize->400],
Dynamic[
Show[{Graphics3D[{Black,Thick,Dashed,bathtub,Blue,Opacity[.1],cone1,cone2,cone7,Opacity[.25],LightCyan,Air1,Red,Z1,Z2,Z7,Blue,river,waterFlowIn,waterFlowOut,waterFixpoint,fixpoint,Yellow,rockbed,Green,Opacity[0.1],soil,Gray,Antenna,WTC7roofBoxHigh, WTC7roofBoxLow,WTC1,WTC2,Orange,Opacity[0.2],valley4,LightPurple,Opacity[.5],WTC7base,WTC7roof,WTC7south,WTC7north,WTC7east,WTC7west,
RGBColor[Dynamic[red],Dynamic[green],Dynamic[blue]],
Red,Opacity[0.5],Yellow,Opacity[0.9],Sphere[{0,0,0},radiusGasousSphere2[5400 t2]]}, PlotRange->All,Axes->True],


(* Water Flow Glacier Valley *)
ParametricPlot3D[waterFlow,{u,0,0.5Pi}],

(* Glacier Valley; PotHole 1; OLD=Ellipsoid[{50,-80,75},{12.5,12.5,20}]*)
 Plot3D[valley1,{x,37.5,62.5},{y,-92.5,-67.5},Mesh->3,PlotPoints->10],

(* Glacier Valley; PotHole 2; OLD = Sphere[{92.5,-62.5,75},20] *)
Plot3D[valley2,{x,72.5,112.5},{y,-82.5,-42.5},Mesh->3,PlotPoints->10],

(* Glacier Valley; PotHole 3; OLD = Sphere[{132.5,-80,75},20] *)
Plot3D[valley3,{x,112.5,152.5},{y,-100,-60},Mesh->3,PlotPoints->10],

aE1=6;bE1=4;cE1=2;
aJ1=.5;bJ1=.5;

(* Plasma Eruption Cone 1 *)
 (* ContourPlot3D[{(x+12.5)^2/aJ1^2+(y-137)^2/bJ1^2== t2(480-z) },{x,-300,300},{y,-300,300},{z,75,500},Mesh->3,PlotPoints->20], *)

(* Plasma Eruption Ellipses
ContourPlot3D[(x+12.5)^2/aE1^2+(y-137)^2/bE1^2+(z-430)^2/cE1^2==80,{x,-500,500},{y,-500,500},{z,-500,500},Mesh->5,PlotPoints->10], *)

(* Plasma Eruption Cone 2 *)
 ContourPlot3D[{(x)^2/aJ1^2+(y)^2/bJ1^2== t2(480-z) },{x,-300,300},{y,-300,300},{z,75,500},Mesh->3,PlotPoints->20],


(* Melt Funnel WTC1 *)
ParametricPlot3D[{-12.5 - a*z* Sin[θ],137- a*z * Cos[θ],(d+
-1.5* (3.8+(38-3.8)t1)*(E^(-HiCurve(b * z)^3))-.5*(3.8+(38-3.8)t1)* (E^(-LoCurve(b * z)^2))-4* MinimaHi-4*MinimaLo)}, (* I've inserted (3.8+(38-3.8)t1) instead of "c" here. *)
{θ,FunnelCutOut,2π-FunnelCutOut},{z,0,1.5},PlotPoints->10,Mesh->Nmesh,ViewPoint->{0,3,1.5},Axes->True,Boxed->True,PlotRange->All,PlotLabel->"Click and drag to change viewpoint"],



(* Melt Funnel WTC2 *)
ParametricPlot3D[{a*z* Sin[θ],a*z * Cos[θ],(d+
-1.5* (3.8+(38-3.8)t1)*(E^(-HiCurve(b * z)^3))-.5*(3.8+(38-3.8)t1)* (E^(-LoCurve(b * z)^2))-4* MinimaHi-4*MinimaLo)}, (* I've inserted (3.8+(38-3.8)t1) instead of "c" here. *)
{θ,FunnelCutOut,2π-FunnelCutOut},{z,0,1.5},PlotPoints->10,Mesh->Nmesh,ViewPoint->{0,3,1.5},Axes->True,Boxed->True,PlotRange->All,PlotLabel->"Click and drag to change viewpoint"],



(* Melt Funnel WTC7 *)
ParametricPlot3D[{145 - a*z* Sin[θ],272.5- a*z * Cos[θ],(d+
-1.5* (3.8+(38-3.8)t1)*(E^(-HiCurve(b * z)^3))-.5*(3.8+(38-3.8)t1)* (E^(-LoCurve(b * z)^2))-4* MinimaHi-4*MinimaLo)}, (* I've inserted (3.8+(38-3.8)t1) instead of "c" here. *)
{θ,FunnelCutOut,2π-FunnelCutOut},{z,0,1.5},PlotPoints->10,Mesh->Nmesh,ViewPoint->{0,3,1.5},Axes->True,Boxed->True,PlotRange->All,PlotLabel->"Click and drag to change viewpoint"]


}, ImageSize->750]]}}]}],


Initialization:> (rockbed=Polygon[{{-500,-500,-200},{-500,500,-200},{500,500,-200},{500,-500,-200}}];


Air1= Cylinder[{{-12.5,137,455},{-12.5,137,510}},32];
Air2=Cylinder[{{0,0,420},{0,0,510}},32];
Air7= Cylinder[{{145,272.5,125},{145,272.5,275}},20];
Antenna=Cylinder[{{-12.5,137,517},{-12.5,137,627}},3];

bathtub=Line[{{-52.5,312.5,100},{-62.5,255,100},{-65,227.5,100},{-85,97.5,100},{-92.5,62.5,100},{-105,-12.5,100},{32,-87.5,100},{50,-27.5,100},{75,115,100},{90,180,100},{112.5,222,100},{-52.5,312,100}}];

cone1=Cone[{{-12.5,137,475},{-12.5,137,0}},10];
cone2=Cone[{{0,0,475},{0,0,0}},10];
cone7=Cone[{{145,272.5,350},{145,272.5,0}},10];
fixpoint=Sphere[{500,500,825},0.1];
river=Triangle[{{-500,-100,90},{-415,500,90},{-500,500,90}}];
soil=Cuboid[{-500,-500,75},{500,500,90}];

valley4=GeometricTransformation[Ellipsoid[{115,-20,75},{40,20,10}],RotationTransform[-25 °,{0,0,1},{115,-20,75}]];

waterFlow={120 + 200 Cos[u],-80 + 100Sin[u],75};
waterFlowIn=Sphere[{122.5,-100,75},5];
waterFlowOut=Sphere[{-2.5,-95, 75},5];
waterFixpoint=Sphere[{15,-205, 75},1];
WTC1=GeometricTransformation[Cuboid[{-44.5,105,100},{19.5,169,517}],RotationTransform[-29 Degree,{0,0,1},{-12.5,137,275}]];
WTC2=GeometricTransformation[ Cuboid[{-32,32,100},{32,-32,515}],RotationTransform[-29 Degree,{0,0,1},{0,0,275}]];
WTC7base=Polygon[{{65,337.5,100},{175,277.5,100},{137.5,232.5,100},{52.5,280,100},{65,337.5,100}}];
WTC7roof=Polygon[{{65,337.5,290},{175,277.5,290},{137.5,232.5,290},{52.5,280,290},{65,337.5,290}}];
WTC7roofBoxHigh=GeometricTransformation[Cuboid[{115,252.5,290},{155,287.5,300}],RotationTransform[-29 Degree,{0,0,1},{135,270,300}]];
WTC7roofBoxLow=GeometricTransformation[Cuboid[{77.5,285,290},{122.5,310,295}],RotationTransform[-29 Degree,{0,0,1},{100,297.5,295}]];
WTC7north=Polygon[{{65,337.5,100},{175,277.5,100},{175,277.5,290},{65,337.5,290},{65,337.5,100}}];
WTC7south=Polygon[{{137.5,232.5,100},{52.5,280,100},{52.5,280,290},{137.5,232.5,290},{137.5,232.5,100}}];
WTC7east=Polygon[{{65,337.5,100},{52.5,280,100},{52.5,280,290},{65,337.5,290},{65,337.5,100}}];
WTC7west=Polygon[{{175,277.5,100},{137.5,232.5,100},{137.5,232.5,290},{175,277.5,290},{175,277.5,100}}];
Z1= Sphere[{-12.5,137,0},5];
Z2= Sphere[{0,0,0},5];
Z7= Sphere[{145,272.5,0},5];
slowEnergyInput2[t1_]:=((k1 Q t1)/(4.18*^12));radiusHeatZone2[t1_]:=(((0.75 k1 k2 Q t1)/(π 790 500 2620))^(1/3));radiusGasousSphere2[t2_]:=(((0.75 k1 k2 Q t2)/(π 790 5000 2620))^(1/3));
k1=1;
k2=1;Q=5.8*^10;



(* Original title of the mathematical funnel model: frustrated energy landscapes (folding funnels); Source/Credit: Terry Oas (oas@duke.edu) *)
(* Source/Credit: https://oaslab.com/drawing_funnels.html  *)
(* *)

SeedRandom[5];(*to get the same funnel every time, leave seed fixed; to get a new funnel, change this number, to get a random funnel, remove number from []*)
       Nmesh=10;(*Number of mesh lines drawn on funnel, the more lines, the slower rendering and interactive rotation*)
HiCurve=16;(*Curvature of the upper part of funnel, try 2*)
LoCurve=8;(*Curvature of the lower part of funnel, try 8*)
BottomRadius=0.2; (*the radius (0-1.5) of the part of the funnel with no local minima to make a smooth native well*)
MiddleRadius=1.5; (*part of funnel where the number and size of minima changes from the Lo region to the Hi region*)
TopRadius=0.5;(*maximum radius where local minima occur.  To avoid having the rim of funnel cut through minima, make < 1.5*)
FunnelCutOut=0.1π; (*angular size (in radians) of arc cut out of funnel to allow side view of interior*)
StartAngle=0.25π; (*lowest angle where minima occur.  To avoid cutting through local wells with cut-out, make > FunnelCutOut*)
StopAngle=1.75π;(*highest angle where minima occur.  To avoid cutting through local wells with cut-out, make < 2π- FunnelCutOut*)

(*Define appearance of local minima on Lo region of funnel*)
NhiMinima=10; (*number of local minima, try 100*)
HiMinimaScale={1.6,2.6}; (*range of minima well depths*)
HiMinimaRadialWidth=0.08; (*width of minima in radial direction*)
HiMinimaTangentWidth=0.001; (*width of minima pepindicular to radius*)

(*Define appearance of local minima on Hi region of funnel*)
NloMinima=40; (*number of local minima, try 20*)
LoMinimaScale={0.2,0.3}; (*range of minima well depths*)
LoMinimaRadialWidth=0.01; (*width of minima in radial direction*)
LoMinimaTangentWidth=0.001; (*width of minima pepindicular to radius*)

MinimaHi=0;
       MinimaLo=0;
a = 65; (* increasing top view area *)
       b= 3.5; (* increasing suddenness of chasm *)

d= 75; (* increasing vertical position (height) of the entire funnel; 'Zero Point' is 75 m below bedrock *)

SetSystemOptions["CheckMachineUnderflow"->False];
SynchronousUpdating->Automatic;
Do[MinimaHi+=Random[Real, HiMinimaScale]((E^(-((((((θ-RandomReal[{StartAngle,StopAngle}])/(2 π))))^2)/HiMinimaTangentWidth)))) ((E^(-((((z-RandomReal[{MiddleRadius,TopRadius}]))^2)/HiMinimaRadialWidth)))),{NhiMinima}];
        Do[MinimaLo+=Random[Real, LoMinimaScale]((E^(-((((((θ-RandomReal[{StartAngle,StopAngle}])/(2 π))))^2)/LoMinimaTangentWidth)))) ((E^(-((((z-RandomReal[{BottomRadius,MiddleRadius}]))^2)/LoMinimaRadialWidth)))),{NloMinima}];

e1= 0.625; (* definition of ellipsoid 1 for PotHole 1; x and y radius are equal*)
valley1 =75 -Sqrt[400-(x-50)^2/e1^2-(y+80)^2/e1^2];
valley2  = 75 -Sqrt[400-(x-92.5)^2-(y+62.5)^2];
valley3 = 75 -Sqrt[400-(x-132.5)^2-(y+80)^2];

k1=1; k2=1; Q= 11.6*^10;

slowEnergyInput2[t1_]:= k1 *Q* t1/4.18*^12;

radiusHeatZone2Temp100[t1_]:=CubeRoot[(0.75*k2* k1* Q*t1)/(π* 790*100*2620)];
radiusHeatZone2Temp1200[t1_]:=CubeRoot[(0.75*k2* k1* Q*t1)/(π* 790*1200*2620)];
radiusHeatZone2Temp2000[t1_]:=CubeRoot[(0.75*k2* k1* Q*t1)/(π* (790*2000 + 418400)*2620)];

)]

(* Parameters *)
(* Density of granite: ro = 2620 kg/m³ | Heat capacity of granite: c = 790 J/kgK | Enthalpie of fusion of granite: H = 420000 J/kg  *)
(* Thermal conductivity of granite: lambda = 2.8 Watt/MeterKelvin; Source: https://de.wikipedia.org/wiki/Wärmeleitfähigkeit *)

(* Question 1: how many kg of granite at 20 °C can be heated up to 1220 °C  when very slowly absorbing 150 kt ? *)
(* Solution 1: Heat Q = 150 kt = 150 x 4.18 10^12 J | Delta T = 1200 °K *)
(* c = Q/(mass x Delta T)  | mass = Q/(c x Delta T) =  6.27 10^14 /(790 x 1200) = 6.6 x 10^8 kg *)

(* Question 2: this corresponds to which volume ? *)
(* Solution 2: ro = m/V  | V = 6.6 x 10^8 kg / (2620 kg/m³) = 252440 m³ *)

(* Question 3: this corresponds to which radius of a sphere ? *)
(* Solution 3: V = 4/3 * π r^3  | radiusHeatZone2[t_]:=CubeRoot[(0.75*V)/π; = 252440 m³] = 39 m*)

(* Question 4: how big is the sphere when the granite is in addition liquefied ? *)
(* Solution 4: __________Use Enthalpie of fusion of granite: H = 420000 J/kg  _____________________*)
(* https://ocw.mit.edu/courses/earth-atmospheric-and-planetary-sciences/12-109-petrology-fall-2005/lecture-notes/Nov3notes.pdf *)


(* Energy-Input WTC2 (PreFizzle): 75 kt in 90 min (= 0.0138889 kt per second *)
(* 0.0138889 kt per second = 5.8 10^10 Joule per second (58.000 MW) *)
(* 1 kt = 4.18 10^12 J *)

(* k1 = damping factor (loss by radiation); k1 ≤ 1) *)
(* k2 = damping factor (loss by fusion enthalpie); k2 ≤ 1)*)

(* Energy-Input WTC2 (Total): 150 kt in 90 min (= 0.02777 kt per second *)