u/mccann: planet.m

(* ::Package:: *)

(* ::Input:: *)
(*ClearAll["Global`*"]*)


(* ::Subsubsection:: *)
(*Constants*)


(* ::Input:: *)
(*GN=6.67259*10^-8;*)
(*kb = 1.380658*10^-16;*)
(*kbEV=8.6173303*10^-5;*)
(*me=9.1093897*10^-28;*)
(*h=6.6260755*10^-27;*)
(*\[Sigma]=6.3*10^-18 (16/13.6)^-3;*)
(*mH=1.67*10^-24;*)
(*\[Gamma]=5/3;*)
(*n=1/(\[Gamma]-1); (*Scaling heights for a given power law n from: Subscript[X, n](r)=Subscript[X, 0]f(r)^n*)*)


(* ::Subsubsection:: *)
(*Problem setup*)


(* ::Input:: *)
(*M=  .5*10^30;*)
(*R0=1.5 10^10;*)
(*T0 =1100 ;*)
(*x0=0;*)
(*\[Mu]0=mH/(1+x0);*)
(*adist=(1 10^12)/R0*)
(*Ms= 1.989 10^33;*)
(*Rs=(6.957 10^10)/R0;*)
(*\[CapitalOmega] = Sqrt[(GN(Ms+M))/(adist R0)^3];*)
(*\[Phi][r_]=(-GN M)/(Sqrt[r^2] R0)- (GN Ms)/(Sqrt[(adist+r)^2]R0)-(\[CapitalOmega](adist+r)R0)^2/2;*)
(*T\[Phi][x_,y_,z_]=(-GN M)/(Sqrt[x^2+y^2+z^2] R0)- (GN Ms)/(Sqrt[(adist+x)^2+y^2+z^2]R0)-((\[CapitalOmega] R0)^2 ((adist+x)^2+y^2))/2;*)
(*dxMIN=1/64;(*In terms of R0*)*)
(*Courant=0.4;*)


(* ::Subsubsection:: *)
(*Radius defintions*)


(* ::Input:: *)
(*Rc = ((\[Gamma]-1)GN M \[Mu]0)/(\[Gamma] kb T0) 1/R0;*)
(*Rz=Rc/(Rc-1)*)
(*RHill=adist (M/(3 Ms))^(1/3);*)
(*TRz=x/.FindRoot[T\[Phi][x ,0,0]==T\[Phi][-1,0,0]+\[Gamma]/(\[Gamma]-1) (kb T0)/\[Mu]0,{x,-1}]*)


(* ::Subsubsection:: *)
(*Scale height functions*)


(* ::Input:: *)
(*Ne[a_,b_]=n*Log[a/b (b-Rz)/(a-Rz)];*)
(*F[r_]=r/n (1-r/Rz);*)
(*H[r_,s_]=((1-s^(-1/n))(Rz-r))/(s^(-1/n) (Rz-r)+r) r; (*folding factor >1*)*)
(*G[r_]=(Rz/r-1)^-n Rz((Gamma[1-n]*Gamma[1+n])/Gamma[2]-Beta[r/Rz,1-n,1+n]);*)


(* ::Subsubsection:: *)
(*Optical depth one defined at R0*)


(* ::Input:: *)
(*n0=1/(\[Sigma] G[1]R0)*)
(*\[Rho]0=\[Mu]0 n0 ;*)
(*P0 = n0 kb T0;*)
(*cs0 = Sqrt[( \[Gamma] P0)/\[Rho]0]; (*adiabatic*)*)


(* ::Subsubsection:: *)
(*Variable profiles*)


(* ::Input:: *)
(*Tf[x_,y_,z_]=1+(\[Gamma]-1)/\[Gamma] \[Mu]0/(kb T0) (T\[Phi][-1,0,0]-T\[Phi][x ,y,z]);T\[Rho][x_,y_,z_] = \[Rho]0 Tf[x,y,z]^(1/(\[Gamma]-1));*)
(*TP[x_,y_,z_]=P0 Tf[x,y,z]^(\[Gamma]/(\[Gamma]-1));*)
(*TT[x_,y_,z_]=T0 Tf[x,y,z];*)
(*Tcs[x_,y_,z_]=cs0 Sqrt[Tf[x,y,z]];*)
(*f[r_]=1+(\[Gamma]-1)/\[Gamma] \[Mu]0/(kb T0) (\[Phi][-1]-\[Phi][r]);*)
(*\[Rho][r_] = \[Rho]0 f[r]^(1/(\[Gamma]-1));*)
(*P[r_]=P0 f[r]^(\[Gamma]/(\[Gamma]-1));*)
(*T[r_]=T0 f[r];*)
(*cs[r_]=cs0 Sqrt[f[r]];*)
(*rnorm=-1;*)
(*LogPlot[{T\[Rho][x,0,0]/T\[Rho][rnorm,0,0],TP[x,0,0]/TP[rnorm,0,0],TT[x,0,0]/TT[rnorm,0,0],Tcs[x,0,0]/Tcs[rnorm,0,0]},{x,.6,Rz},PlotLegends->"Expressions",GridLines->{{{Rz,Black}},{}}]*)


(* ::Subsubsection:: *)
(*Thermal ionization*)


(* ::Input:: *)
(*sahaA=T0^(3/2)/n0 ((2 \[Pi] me kb)/h^2)^(3/2) Rz/Rc;*)
(*sahaB=13.6/(kbEV T0) Rz/Rc;*)
(*X[r_]=(sahaA r/(Rz-r) Exp[-sahaB r/(Rz-r)])/2 (Sqrt[1+4/(sahaA r/(Rz-r)) Exp[sahaB r/(Rz-r)]]-1);*)
(*Plot[{X[r],((1-X[r])*T\[Rho][r,0,0])/T\[Rho][rnorm,0,0]},{r,.6,Rz}]*)
(*r/.FindRoot[X[r]==.5,{r,.8}]*)
(*ListPlot[Table[{i ,NIntegrate[\[Sigma]  ((1-X[r])T\[Rho][r,0,0])/\[Mu]0 R0,{r,i,Rz}]},{i,.8%,Rz,.025}]]*)


(* ::Subsubsection:: *)
(*Scale height functions*)


(* ::Input:: *)
(*Hiso[r_]=(kb (r R0)^2 T[r])/(GN M \[Mu]0) 1/R0;*)
(*Plot[{F[x],H[x,E],G[x],Hiso[x]},{x,0,Rz},PlotLegends->"Expressions"]*)


(* ::Subsubsection:: *)
(*Inner boundary & mask*)


(* ::Input:: *)
(*r/.NSolve[\[Rho][r]==\[Rho][-1] E^2,r][[6]]*)
(*rib=Floor[Ceiling[2 %/dxMIN]/2]dxMIN*)
(*rmask=rib+5dxMIN*)


(* ::Input:: *)
(*N[rib]*)
(*N[rmask]*)


(* ::Subsubsection:: *)
(*outter boundary*)


(* ::Input:: *)
(*rob=TRz+dxMIN/2*)
(*T\[Rho][rob,0,0]/\[Rho][-1]*)


(* ::Subsubsection:: *)
(*Ratio of density between cell centers*)


(* ::Input:: *)
(*\[CapitalDelta]R =Table[Style[{i dxMIN,\[Rho][(i-1) dxMIN]/\[Rho][i dxMIN]},Hue[-.125+.375 \[Rho][(i-1) dxMIN]/\[Rho][i dxMIN]]],{i,Floor[-rmask/dxMIN]-.5,-rob/dxMIN}]*)
(*Show[Plot[{F[x]/F[Rz/2]+1},{x,-rmask,Rz},PlotLegends->"Expressions",GridLines->{{{-rmask,Thick},{-rib,Blue},{1,Orange},{-rob,Green}},{{1.2397783001618599,Cyan},{1.137042892629939,Cyan},{1.5003432045110259,Red},{3.015323228152034,Red}}},GridLinesStyle->Directive[ Dashed]],ListPlot[\[CapitalDelta]R],PlotRange->{{-rmask-dxMIN,-rob+dxMIN},{1,Max[3,\[CapitalDelta]R[[Floor[-rob/dxMIN]-Floor[-rmask/dxMIN]+1]][[1]][[2]]]}}]*)