(* Define Parameters *)
(SubStar[τ])=1;
R0=1.5;
q=0;
β=1;

(* Optical Depth *)
t[p_,z_]:=Piecewise[{{∞,Abs[p]≤1&&z≤(Sqrt[1-(p^2)])}},NIntegrate[((SubStar[τ])/(((p^2)+(Z^2)) (((1-(1/(Sqrt[(p^2)+(Z^2)]))))^β))),{Z,z,∞}]]

(* Plot contours of optical depth = .25, .5, 1, 2, and 4 *)
contours=ContourPlot[t[p,-z],{z,-10,10},{p,-10,10},Contours->{.33,1,3},PlotRange->{0,7},ContourShading->None,ContourStyle->{{Green},{White,Dashed},{White,Dotted}}]

(* Radial Velocity, as seen by observer on the left *)
vrad[x_,y_]:=Piecewise[{{0,(Sqrt[(x^2)+(y^2)])<1}},(((1-(1/(Sqrt[(x^2)+(y^2)]))))^β) ((-x)/(Sqrt[(x^2)+(y^2)]))]
(* Plot radial velocity in terms of Doppler shift *)
shift=DensityPlot[vrad[x,y],{x,-10,10},{y,-10,10},ColorFunction->(Blend[{CMYKColor[0,1,1,0],CMYKColor[0,1,0,0],CMYKColor[1,1,0,0]},#]&),PlotPoints->70,Exclusions->None]
(* For Blue/Magenta/Red: 
(Blend[{CMYKColor[0,1,1,0],CMYKColor[0,1,0,0],CMYKColor[1,1,0,0]},#]&) *)

(* Define a color scheme that scales opacity *)
BlackWhite[z_]:=RGBColor[0,0,0,1-z]
(* Emissivity *)
emissivity[x_,y_]:=(((1-(1/(Sqrt[(x^2)+(y^2)]))))^(-2 β)) ((((x^2)+(y^2)))^(-2-q))
(* Plot the emissivity, log 10 *)
extinction=DensityPlot[Log[10,emissivity[x,y]],{x,-10,10},{y,-10,10},ColorFunction->BlackWhite,PlotPoints->70]

(* Create a white screen to make the plot circular *)
screen=RegionPlot[(Sqrt[(x^2)+(y^2)])≥10,{x,-15,15},{y,-15,15},PlotStyle->White,PerformanceGoal->"Quality"];
(* Create a border *)
border=Graphics[{Gray,Circle[{0,0},10.25]}];

(* Create the star image and the black circle inside R0 *)star=Graphics[Table[{Blend[{RGBColor[.4,.8,1],White},(x^(.7))],Disk[{0,0},1-x]},{x,0,1,1/10}]];
r0=Graphics[{Black,Disk[{0,0},R0]}];

(* Create and position labels for the optical depth contours *)
labels=Graphics[{White,Text[Style["4",FontSize->14],{18.5,2}],Text[Style["2",FontSize->14],{18,4.1}],Text[Style["1",FontSize->14],{16.5,7.5}],Text[Style["0.5",FontSize->14],{12,13}],Text[Style["0.25",FontSize->14],{1,18}]}];

(* Create parameter table *)
plotlabel=Graphics[{Black,Text[Style[Column[{Row[{"(SubStar[τ]) = ",(SubStar[τ])}] ,Row[{"(Subscript[R, 0]) = ",R0}],Row[{"β = ",β}],Row[{"q = ",q}]}],FontSize->13],{-9,8}]}];

(* All together *)
Show[shift,extinction,contours,r0,star,screen,labels,border,plotlabel,Frame->None]