# choose "pdf" or "png", defaults to "pdf"
type = "pdf"

if(type eq "png") {
print "Make pngs"
set terminal pngcairo transparent font "Helvetica" fontscale 3 linewidth 5 size 1900,1200
set output "disc.pdf"
} else {
print "Make pdfs"
set terminal pdfcairo font "Helvetica" fontscale 2 linewidth 5 size 16,12
set output "disc.pdf"
}

print "Making 2D plots"

# 2D plots

set format y "%g"
set format x "%g"
set logscale x
set xlabel "Time / years"
set yrange[0:4]

if(type eq "png") set output "disc_Nzones.png"

plot "cbdisc.log" u 't':'nzones' w l title "N_{zones}" 

set logscale y
set yrange[*:*]
set key top left
set xrange[1000:1e5]

if(type eq "png") set output "disc_radii.png"
plot \
"cbdisc.log" u 't':'Rin' w l title "R_{in}/R_{⊙}",\
'cbdisc.log' u 't':'Rout' w l title "R_{out}/R_{⊙}"

if(type eq "png") set output "disc_mass.png"
plot \
'cbdisc.log' u 't':'M' w l title "Disc Mass / M_{⊙}"

if(type eq "png") set output "disc_angmom.png"
plot \
'cbdisc.log' u 't':'J' w l title "Disc Angular Momentum (cgs)"

if(type eq "png") set output "disc_moment_of_inertia.png"
plot \
'cbdisc.log' u 't':'I' w l title "Disc Moment of Inertia (cgs)"

if(type eq "png") set output "disc_luminosity.png"
plot \
'cbdisc.log' u 't':'L' w l title "Disc Luminosity (L_{⊙})"

if(type eq "png") set output "disc_luminosity_ratio.png"
plot \
'cbdisc.log' u 't':'L/L*' w l title "Disc Luminosity / Stellar Luminosity"

if(type eq "png") set output "disc_temperatures.png"
plot \
"cbdisc.log" u 't':'Tin' w l title "T_{in}/K",\
'cbdisc.log' u 't':'Tout' w l title "T_{out}/K"

if(type eq "png") set output "disc_rho.png"
plot \
"cbdisc.log" u 't':'rhoin' w l title "rho_{in}/K",\
'cbdisc.log' u 't':'rhoout' w l title "rho_{out}/K"

if(type eq "png") set output "disc_H_Rin.png"
plot \
'cbdisc.log' u 't':'Hin/Rin' w l title "H_{in}/R_{in}"

if(type eq "png") set output "disc_viscous_timescale_outer_edge.png"
plot \
'cbdisc.log' u 't':'tviscout' w l title "Viscous timescale at R_{out}"

############################################################
# 3D surface plots
############################################################
print "Making 3D plots"
set logscale xycb
set xlabel 'Time / years'
set ylabel 'Radius / R_{⊙}'
set format cb '10^{%T}'
set pm3d map clip1in
set palette cubehelix
set lmargin at screen 0.15
set rmargin at screen 0.87

t_filter(t) = MAX(1000,t)

if(type eq "png") set output "disc_map_T.png"
splot \
'< gawk "{if(\$1<1e5){print}}" cbdisc.dat' u (t_filter(column('t'))):'R':'T' title "Temperature / K"

if(type eq "png") set output "disc_map_L.png"
splot \
'< gawk "{if(\$1<1e5){print}}" cbdisc.dat' u (t_filter(column('t'))):'R':'L' title "L(R) / L_{⊙}"

if(type eq "png") set output "disc_map_Sigma.png"
splot \
'< gawk "{if(\$1<1e5){print}}" cbdisc.dat' u (t_filter(column('t'))):'R':'Sigma' title "{/Symbol S} / g/cm^{2}"

if(type eq "png") set output "disc_map_H_on_R.png"
splot \
'< gawk "{if(\$1<1e5){print}}" cbdisc.dat' u (t_filter(column('t'))):'R':'H/R' title "H / R"

if(type eq "png") set output "disc_map_cs.png"
splot \
'< gawk "{if(\$1<1e5){print}}" cbdisc.dat' u (t_filter(column('t'))):'R':(column('cs')/1e5) title "c_{s} / km/s"

if(type eq "png") set output "disc_map_P.png"
splot \
'< gawk "{if(\$1<1e5){print}}" cbdisc.dat' u (t_filter(column('t'))):'R':'P' title "Pressure"

if(type eq "png") set output "disc_map_F.png"
splot \
'< gawk "{if(\$1<1e5){print}}" cbdisc.dat' u (t_filter(column('t'))):'R':'F' title "Ang Mom Flux, F"

if(type eq "png") set output "disc_map_tvisc.png"
splot \
'< gawk "{if(\$1<1e5){print}}" cbdisc.dat' u (t_filter(column('t'))):'R':'tvisc' title "Viscous Timescale / s"