#* Draw a TS diagram, possibly choosing isopycnals automatically
#* @param $rho_inc (optional), density step between isopycnals
cmd draw TS diagram (;$) { # [$rho_inc]
my $rho; my $rho_min; my $rho_max; my $rho_inc; my $rho_inc_range;
my $S; my $dS;
set x name "Salinity, PSU";
set y name "Temperature, \(\circ\)C";
draw axes;
$rho_min = sw_density($_xleft, $_ytop, 0); # min density on diagram
$rho_max = sw_density($_xright, $_ybottom, 0); # max density on diagram
if ($#_ == 0) {
# Use specified density increment.
$rho_inc = abs($_[0]);
} else {
# Didn't give an increment, so pick one, using 1-2-5 scaling
# to give of the order of 5 isopycnals in range of data.
$rho_inc = ($rho_max - $rho_min) / 5;
$rho_inc_range = exp10(floor(log10($rho_inc))); # a power of 10
if ($rho_inc / $rho_inc_range < 2) {
$rho_inc = $rho_inc_range;
} else {
if ($rho_inc / $rho_inc_range < 5) {
$rho_inc = $rho_inc_range * 2;
} else {
$rho_inc = $rho_inc_range * 5;
}
}
}
$rho_min = $rho_inc * floor($rho_min / $rho_inc);
$rho_max = $rho_inc * ceil($rho_max / $rho_inc);
#DEBUG: print "$rho_min $rho_max $rho_inc\n";
for ($rho = $rho_min; $rho <= $rho_max; $rho += $rho_inc) {
draw isopycnal $rho;
}
if (sw_freezing_temperature($_xleft, 0) > $_ybottom) {
draw line from $_xleft sw_freezing_temperature($_xleft, 0)
to $_xright sw_freezing_temperature($_xright, 0);
}
set color blue;
set symbol size 0.05;
draw symbol bullet;
set color black;
# Draw freezing-point line
$dS = ($_xright - $_xleft) / 50;
for ($S = $_xleft; $S <= $_xright; $S += $dS) {
draw line from ($S-$dS) sw_freezing_temperature($S-$dS, 0)
to $S sw_freezing_temperature($S, 0);
}
}
$file = "./TS-stn308.gz"; # from EUBEX dataset
open(IN, "zcat $file |") or die "cannot access $file\n";
read columns IN x=3 y=2;
draw TS diagram;
#draw TS diagram 1;