(*Computation of the bending function and the maximum bending moment for point load*)
e=9*10^9;
ν=0.25;
g=0.113;
ρ=1020;
Te=5000;
V=-300000000;
d=((e*Te^3)/(12*(1-ν^2)));
α= Sqrt[Sqrt[4*d/(ρ*g)]];
xm=α*(π/4); 
w[x_]:=((V*α^3)/(4d))*Exp[-x/α]*(Cos[x/α]);

(*The following expression are identical, this is done for comparison*)
M[x_]:=(V/2)*α*Exp[((-x)/α)]Sin[(x/α)];
derM1[x_]:=w''[x]

point1={33898.87,w[33898.87]};

point1s={33898.87,s[33898.87]};
point2s={34228.91,s[34228.91]};
point1M={33898.87,-M[33898.87]};
point2M={34228.91,M[34228.91]};
point1derM1={33898.87,-derM1[33898.87]};
point2derM1={34228.91,derM1[34228.91]};

Plot[w[x],{x,0,200000},Epilog->{PointSize[Large],Red,Point[point1]},AxesLabel->{"x[m]","w(x)[m]"},AxesLabel->{"x[m]","w(x)[m]"}]

Plot[-M[x],{x,0,200000},Epilog->{PointSize[Large],Red,Point[point1M],Dashed,Black,Line[{{xm,0},{xm,-M[xm]}}],Text[Style["x_m",Black,8],{xm,-1},{0,1}]},AxesLabel->{"x[m]","M(x)[Nm]"}]

Plot[-derM1[x],{x,0,200000},PlotRange->All,Epilog->{PointSize[Large],Red,Point[point1derM1],Blue,Point[point2derM1],Dashed,Black,Line[{{xm,0},{xm,-derM1[xm]}}],Text[Style["x_m",Black,8],{xm,-1},{0,1}]},AxesLabel->{"x","M(x)"},PlotLabel->"Bending moment DM1"]

(* Computations for the bending function obtained by variation of constants for point load as an approximation of the dirac delta distribution by a Gaussian distribution as a goes to 0 and for the distributed load q(x) *)
ClearAll["Global`*"]

(* constants for q(x) *)
BB=-0.448;
CC=3.43*10^(-2);
DD=-1.39*10^(-3);
EE=2.64*10^(-5);
FF=-1.88*10^(-7);

rE=252000; (*radius of Enceladus*)
a=20; (* parameter for the Gaussian *)

(*physical parameters*)
e=9*10^9;
ν=0.25;
g=0.113;
ρ=1020;
Te=5000;
V=-10000;
d=((e*Te^3)/(12*(1-ν^2)));
α= Sqrt[Sqrt[4*d/(ρ*g)]];

(* Specified parameters for distributed load *)
Tev=5000;
dv=((e*Tev^3)/(12*(1-ν^2)));
αv= Sqrt[Sqrt[4*dv/(ρ*g)]];

(* Specified parameters for point load *)
Tepoint=5000;
dpoint=((e*Tepoint^3)/(12*(1-ν^2)));
αpoint= Sqrt[Sqrt[4*dpoint/(ρ*g)]];

(* function for distributed load in radians *)
v[x_]:=Piecewise[{{0,x>90},{Power[10,BB*x+CC*x*x+DD*x*x*x+EE*x*x*x*x+FF*x*x*x*x*x],x<=90}}];

(* function for distributed load in meters*)
q[x_]:=v[x/rE*180/π]*27615.08;

(*approximation of Dirac delta for point load*)
qpoint[x_]:=(300000000/(a*Sqrt[π]))*Exp[-((x/a)*(x/a))]; 


Plot[{q[x],qpoint[x]},{x,0,500000}]

(* integrands for the parameters obtained by the variation (distributed load)*)
intdA[t_]:=((αv^3)/(8d))*Exp[-t/αv]*(Cos[t/αv]+Sin[t/αv])*q[t];
intdB[t_]:=-((αv^3)/(8d))*Exp[-t/αv]*(Cos[t/αv]-Sin[t/αv])*q[t];
intdC[t_]:=((αv^3)/(8d))*Exp[t/αv]*(Cos[t/αv]-Sin[t/αv])*q[t];
intdD[t_]:=((αv^3)/(8d))*Exp[t/αv]*(Cos[t/αv]+Sin[t/αv])*q[t];
(*integrands for the integral constants*)
intdC0[t_]:=((αv^3)/(8d))*Exp[-t/αv]*(3*Cos[t/αv]-Sin[t/αv])*q[t];
intdD0[t_]:=-((αv^3)/(8d))*Exp[-t/αv]*(Cos[t/αv]-Sin[t/αv])*q[t];

(* integrands for the parameters obtained by the variation (point load) *)
intdApoint[t_]:=((αpoint^3)/(8d))*Exp[-t/αpoint]*(Cos[t/αpoint]+Sin[t/αpoint])*qpoint[t];
intdBpoint[t_]:=-((αpoint^3)/(8d))*Exp[-t/αpoint]*(Cos[t/αpoint]-Sin[t/αpoint])*qpoint[t];
intdCpoint[t_]:=((αpoint^3)/(8d))*Exp[t/αpoint]*(Sin[t/αpoint]+Cos[t/αpoint])*qpoint[t];
intdDpoint[t_]:=((αpoint^3)/(8d))*Exp[t/αpoint]*(Sin[t/αpoint]-Cos[t/αpoint])*qpoint[t];
intdC0point[t_]:=((αpoint^3)/(8d))*Exp[-t/αpoint]*(3*Cos[t/αpoint]-Sin[t/αpoint])*qpoint[t];
intdD0point[t_]:=-((αpoint^3)/(8d))*Exp[-t/αpoint]*(Cos[t/αpoint]-Sin[t/αpoint])*qpoint[t];

(* corresponding integrals and integration constants for distributed load *)
cC0=NIntegrate[intdC0[t],{t,0,300000}];
cD0=NIntegrate[intdD0[t],{t,0,300000}];
cA[x_]:=NIntegrate[intdA[t],{t,x,300000}];
cB[x_]:=NIntegrate[intdB[t],{t,x,300000}];
cC[x_]:=NIntegrate[intdC[t],{t,0,x}]+cC0;
cD[x_]:=NIntegrate[intdD[t],{t,0,x}]+cD0;

(* corresponding integrals and integration constants for point load *)
cC0point=NIntegrate[intdC0point[t],{t,0,300000}];
cD0point=NIntegrate[intdD0point[t],{t,0,300000}];
cApoint[x_]:=NIntegrate[intdApoint[t],{t,x,300000}];
cBpoint[x_]:=NIntegrate[intdBpoint[t],{t,x,300000}];
cCpoint[x_]:=NIntegrate[intdCpoint[t],{t,0,x}]+cC0point;
cDpoint[x_]:=NIntegrate[intdDpoint[t],{t,0,x}]+cD0point;

(* bending functions for distributed and point load respectively *)
W[x_]:=(cA[x]*Exp[(x/(αv))]*Cos[(x/(αv))]+cB[x]*Exp[(x/(αv))]*Sin[(x/(αv))]+cC[x]*Exp[((-x)/(αv))]*Cos[(x/(αv))]+cD[x]*Exp[((-x)/(αv))]*Sin[(x/(αv))]);
Wpoint[x_]:=(cApoint[x]*Exp[(x/(αpoint))]*Cos[(x/(αpoint))]+cBpoint[x]*Exp[(x/(αpoint))]*Sin[(x/(αpoint))]+cCpoint[x]*Exp[((-x)/(αpoint))]*Cos[(x/(αpoint))]+cDpoint[x]*Exp[((-x)/(αpoint))]*Sin[(x/(αpoint))]);
(*fundamental system*)
w1[x_]:=Exp[x/αv]*Cos[x/αv]
w2[x_]:=Exp[x/αv]*Sin[x/αv]
w3[x_]:=Exp[-x/αv]*Cos[x/αv]
w4[x_]:=Exp[-x/αv]*Sin[x/αv]

(*bending moment and the bending moment derivative in terms of the components of the fundamental system for distributed load*)
BM[x_]:=d*(cA[x]*w1''[x]+cB[x]*w2''[x]+cC[x]*w3''[x]+cD[x]*w4''[x]); 
derBM[x_]:=cA[x]*w1'''[x]+cB[x]*w2'''[x]+cC[x]*w3'''[x]+cD[x]*w4'''[x]; 

(* comparison of two ways how to obtain the maximum bending moment*)
root = FindRoot[W'''[x],{x,10000,100000}]
root2=FindRoot[derBM[x],{x,10000,100000}]

Plot[-W[x],{x,0,200000},Axes->True,AxesLabel->{"x[m]","w(x)[m]"},(*Epilog->{PointSize[Large],Red,Point[{40550.7,-W[40550.7]}]}*)];
(* bending moment function and the corresponding point of the maximum*)
Plot[BM[x],{x,0,100000},PlotRange->All,AxesLabel->{"x[m]","M(x)[Nm]"},Epilog->{PointSize[Large],Red,Point[{40550.7,BM[40550.7]}],Dashed,Black,Line[{{40550.7,0},{40550.7,BM[40550.7]}}],Text[Style["x_r",Black,8],{40550.7,-1},{0,1}]}]


Plot[derBM[x],{x,0,100000},PlotRange->All,AxesLabel->{"x","DM(x)"},PlotLabel->"Bending moment derivative",Epilog->{PointSize[Large],Red,Point[{48537.61,-derBM[48537.61]}]}];
Plot[W'''[x],{x,0,100000},PlotRange->All,AxesLabel->{"x","W'''(x)"},PlotLabel->"Third derivative of W(x)"]

(* bending function and its moment for the mass of the load from Hemingway*)
bendingf[x_]:=((300000000*αv^3)/(4d))*Exp[-x/αv]*(Cos[x/αv]);
M[x_]:=d*bendingf''[x];

Plot[{-W[x],-bendingf[x]},{x,0,200000},Axes->True,AxesLabel->{"x[m]","w(x)[m]"},PlotLegends->{"Distributed load","Point load"},Epilog->{PointSize[Large],Blue,Point[{40550.7,-W[40550.7]}],PointSize[Large],Orange,Point[{33898.87,-bendingf[33898.87]}]}]
Plot[{BM[x],M[x]},{x,0,200000},Axes->True,AxesLabel->{"x[m]","M(x)[Nm]"},PlotLegends->{"Distributed load","Point load"},Epilog->{PointSize[Large],Blue,Point[{40550.7,BM[40550.7]}],PointSize[Large],Orange,Point[{33898.87,M[33898.87]}]}]

(* Plot of the dependence of the position of the new crack on the elastic shell thickness*)
ClearAll["Global`*"]
e=9*10^9;
ν=0.35;
g=0.113;
ρ=930;
K=((π/4)*(e/(3*ρ*g*(1-ν^2))))^((1/4));
c=(((4/π))^4(ρ*g*3*(((1-ν^2))/e)))^((1/3));
xm[Te_]:=Te^((3/4))*K;
Te[x_]:=x^((4/3))*c;

(* Values of the positions of the ew crack for elastic shell thickness from 1-10km  *)
values={{15233.9,1000},{22927.5,2000},{29361.7,3000},{35157,4000},{40550,5000},{45662.1,6000},{50560,7000},{55288.3,8000},{59876.1,9000},{64344,10000}};

listplot = ListPlot[values,AxesLabel->{"Te","xm"},PlotStyle->{Red,PointSize[Medium]}];
functionplot = Plot[Te[x],{x,0,70000},PlotRange->All,Epilog->{Red,Line[{{35000,0},{35000,Te[70000]}}],Dashed,Black,Line[{{0,Te[35000]},{35000,Te[35000]}}]}];

combinedPlot=Show[listplot,functionplot,PlotRange->All,AxesLabel->{"xm[m]","Te[m]"},AxesOrigin->{0,0},Epilog->{PointSize[Medium],Blue,Point[{35000,Te[35000]}],Red,Line[{{35000,0},{35000,Te[50000]}}],Dashed,Black,Line[{{0,Te[35000]},{35000,Te[35000]}}],Dashed,Black,Line[{{0,4000},{35157,4000}}]}];
combinedPlot