#* Draw a TS diagram, possibly choosing isopycnals automatically
#* @param $rho_inc (optional), density step between isopycnals
cmd draw TS diagram (;$) { # [$rho_inc]
        my ($rho, $rho_min, $rho_max, $rho_inc, $rho_inc_range, $S, $dS);
        set x name "Salinity, PSU";
        set y name "Temperature, \(\circ\)C";
        $rho_min = eos_dens($_xleft, $_ytop, 0); # min density on diagram
        $rho_max = eos_dens($_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);
        for ($rho = $rho_min; $rho <= $rho_max; $rho += $rho_inc) {
                draw isopycnal $rho;
        }
        if (eos_freezing_temp($_xleft, 0) > $_ybottom) {
                draw line connecting $_xleft eos_freezing_temp($_xleft, 0)
                        $_xright eos_freezing_temp($_xright, 0);
        }
        draw symbols "bullet";
        # Draw freezing-point line
        $dS = ($_xright - $_xleft) / 50;
        for ($S = $_xleft; $S <= $_xright; $S += $dS) {
                draw line connecting ($S-$dS) eos_freezing_temp($S-$dS, 0)
                                                            $S eos_freezing_temp($S, 0);
        }
}

$file = "/users/dek/bmay/thesis/eubex/data/stn308";
open(IN, "$file") or die "cannot access $file\n";
read columns IN x=3 y=2;
draw TS diagram;
#draw TS diagram 1;