Riemann problem

Riemann solver (this part is done)

(* define function giving Mach number as a function of z ≡ ((Δv)/c) ≥ 0 *)
Clear[mach];
mach[z_,γ_]:=((γ+1)/4)z+(Sqrt[(((((γ+1)/4)z))^2)+1])

(* define function giving ratio ((Subscript[c, c])/c) ≡ ((sound speed at contact interface)/(sound speed outside boundary)) as a function of z ≡ ((Δv)/c), for both -ve and +ve z *)
Clear[ccc];
ccc[z_,γ_]:=If[z≤0,1+((γ-1)/2)z,(Sqrt[1+((γ-1)/2)z(-z+2mach[z,γ])])]

(* define function giving ratio ((Subscript[p, c])/p) ≡ ((pressure at contact interface)/(pressure outside boundary)) given z ≡ ((Δv)/c), for both -ve and +ve z *)
Clear[pcp];
pcp[z_,γ_]:=If[z≤0,((ccc[z,γ])^((2γ)/(γ-1))),1+γ z mach[z,γ]]

(* ((Subscript[p, c])/p) ≥ 0 for z ≥ zmin *)
Clear[zmin];
zmin[γ_]:=-(2/(γ-1))

(* define function that plots ((Subscript[p, c])/p) as a function of z ≡ ((Δv)/c) *)
Clear[plotpcp];
plotpcp[γ_]:=Plot[((pcp[z,γ])^((γ-1)/(2γ))),{z,zmin[γ],8}]

(* plot ((Subscript[p, c])/p) as a function of z ≡ ((Δv)/c) *)
plotpcp[(5/3)]

(* Riemann solver: solve for velocity (Subscript[v, c]) and pressure (Subscript[p, c]) at contact interface by requiring that (Subscript[p, c]) be same for left and right sides *)
Clear[vcpc];
vcpc[{vl_,cl_,pl_},{vr_,cr_,pr_},γ_]:=(
Clear[vc];
vc=vc/.FindRoot[pl pcp[((vl-vc)/cl),γ]==pr pcp[((vc-vr)/cr),γ],{vc,vr+cr zmin[γ],vl-cl zmin[γ]}];
{vc,pl pcp[((vl-vc)/cl),γ]}
)

Solution of the Riemann problem (this part needs editing)

(* define function giving velocity v in non-constant part of rarefaction wave *)
vrarefaction[ξ_,{vl_,cl_},γ_]:=

(* define function giving sound speed c in non-constant part of rarefaction wave *)
crarefaction[ξ_,{vl_,cl_},γ_]:=

(* define function giving pressure p in non-constant part of rarefaction wave *)
prarefaction[ξ_,{vl_,cl_,pl_},γ_]:=

(* value of ξ ≡ (x/t) at boundary, for both rarefaction (vl ≤ vc) and shock (vl ≥ vc) *)
Clear[ξl];
ξl[vc_,{vl_,cl_},γ_]:=If[vl≤vc,,]

(* define function giving {velocity v, sound speed c, pressure p, density ρ} left of the contact interface as a function of ξ ≡ (x/t) *)
Clear[vcpρl];
vcpρl[ξ_,vc_,pc_,{vl_,cl_},γ_]:=(
Clear[v,c,p,ρ,zl,ccl,pl];
zl=((vl-vc)/cl);
ccl=;
pl= ;
{v,c,p}=If[zl≤0,
(* rarefaction wave *)
Which[
ξ≤vl-cl,{vl,cl,pl},
ξ≤vc-ccl,{vrarefaction[ξ,{vl,cl},γ],crarefaction[ξ,{vl,cl},γ],prarefaction[ξ,{vl,cl,pl},γ]},
ξ≤vc,{vc,ccl,pc}],
(* shock wave *)
Which[
ξ≤ξl[vc,{vl,cl},γ],{vl,cl,pl},
ξ≤vc,{vc,ccl,pc}]
];
(* density ρ *)
ρ=((γ p)/(c^2));
(* the list {v,c,p,ρ} returned by this function vcpρl[ξ,vc,pc,{vl,cl},γ] *)
{v,c,p,ρ}
)

(* define function giving {velocity v, sound speed c, pressure p, density ρ} as a function of ξ ≡ (x/t) *)
Clear[vcpρ];
vcpρ[ξ_,vc_,pc_,{vl_,cl_},{vr_,cr_},γ_]:=(
If[ξ≤vc,
(* left of contact interface *)
vcpρl[ξ,vc,pc,{vl,cl},γ],
(* right of contact interface; equations are same as those on left, subject to ξ -> -ξ, v -> -v *)
{-1,1,1,1}vcpρl[-ξ,-vc,pc,{-vr,cr},γ]]
)

Plot (this part is done, but you may prefer to write your own)

(* define function to plot several things side by side *)
Clear[plots];
plots[plotlist_]:=GraphicsRow[plotlist,0,ImageSize->340Length[plotlist]]

(* define a function to plot any of {v, c, p, ρ} as a function of ξ *)
Clear[plotvcpρ];
plotvcpρ[list_,{vl_,cl_,pl_},{vr_,cr_,pr_},γ_,plotlabel_]:=(
Clear[vc,pc];
{vc,pc}=vcpc[{vl,cl,pl},{vr,cr,pr},γ];
Plot[vcpρ[ξ,vc,{vl,cl,pl},{vr,cr,pr},γ][[list]],{ξ,ξl[vc,{vl,cl},γ]-1,-ξl[-vc,{-vr,cr},γ]+1},PlotRange->All,PlotLabel->plotlabel]
)

(* Sod's (1978, JCP 27, 1) parameters *)
plotvcpρ[{1,3},{0,(Sqrt[1.4]),1},{0,(Sqrt[1.12]),.1},1.4,"Velocity and sound speed"]

Module[{γ,vc,vl,cl,pl,vr,cr,pr},
γ=(5/3);
vc=1;
vl=0;
cl=1;
pl=1;
vr=0;
cr=1;
(* this choice of (Subscript[p, r]) yields the specified (Subscript[v, c]) *)
pr=pl((pcp[((vl-vc)/cl),γ])/(pcp[((vc-vr)/cr),γ]));
Print["((Subscript[p, c])/(Subscript[p, l])) = ",Simplify[pcp[((vl-vc)/cl),γ]]," = ",N[pcp[((vl-vc)/cl),γ]]];
Print["((Subscript[p, c])/(Subscript[p, r])) = ",Simplify[pcp[((vc-vr)/cr),γ]]," = ",N[pcp[((vc-vr)/cr),γ]]];
Print[plots[{plotvcpρ[{1},{vl,cl,pl},{vr,cr,pr},γ,"Velocity v"],plotvcpρ[{2},{vl,cl,pl},{vr,cr,pr},γ,"Sound speed c"]}]];
Print[plots[{plotvcpρ[{3},{vl,cl,pl},{vr,cr,pr},γ,"Pressure p"],plotvcpρ[{4},{vl,cl,pl},{vr,cr,pr},γ,"Density ρ"]}]];]

Module[{γ,vc,vl,cl,pl,vr,cr,pr},
γ=(5/3);
vc=-1;
vl=0;
cl=1;
pr=1;
vr=0;
cr=1;
(* this choice of (Subscript[p, r]) yields the specified (Subscript[v, c]) *)
pl=pr((pcp[((vc-vr)/cr),γ])/(pcp[((vl-vc)/cl),γ]));
Print["((Subscript[p, c])/(Subscript[p, l])) = ",Simplify[pcp[((vl-vc)/cl),γ]]," = ",N[pcp[((vl-vc)/cl),γ]]];
Print["((Subscript[p, c])/(Subscript[p, r])) = ",Simplify[pcp[((vc-vr)/cr),γ]]," = ",N[pcp[((vc-vr)/cr),γ]]];
Print[plots[{plotvcpρ[{1},{vl,cl,pl},{vr,cr,pr},γ,"Velocity v"],plotvcpρ[{2},{vl,cl,pl},{vr,cr,pr},γ,"Sound speed c"]}]];
Print[plots[{plotvcpρ[{3},{vl,cl,pl},{vr,cr,pr},γ,"Pressure p"],plotvcpρ[{4},{vl,cl,pl},{vr,cr,pr},γ,"Density ρ"]}]];]