# Draw prob-density-function TS diagram for Bravo data
# This first call to perl is specific to the
# particular (weird) dataset. All that matters
# is that a file of (x,y) data is created, and
# stored in the file called `tmp-xy.\.pid.'
system perl <<"EOF"
#
# Slurp in x[], y[] data
$dir = "/users/dek/kelley/Labrador/bravo/data";
$Sfile = "$dir/S_25db_1day";
$Tfile = "$dir/T_25db_1day";
open(S, "$Sfile") || die "Can't open input $Sfile";
open(T, "$Tfile") || die "Can't open input $Tfile";
open(ST, ">tmp-xy.\.pid.")
|| die "Can't open tmp-xy.\.pid.";
$day = 5;
$year = 1964;
while(<S>) {
@S = split;
$_ = <T>;
@T = split;
if (240 < $day && $day < 360) {
for ($i = 0; $i < $#S; $i++) { #all depths
print ST "$S[$i] $T[$i]\n";
}
}
$day += 1;
if ($day > 365) {
$year++;
$day = 0;
}
if ($year > 1967) { last; }
}
EOF
#
# Create 2D PDF for (x,y) data in file \infile
# storing the results in \outfile. X ranges
# between the indicated limits, with the indicated
# binsize, as does y. The synonyms defined
# on the next 4 lines are the only input this
# perlscript needs; the perlscript itself is
# quite general. For details of what it does,
# particularly the scaling of the PDF, see
# the perl comments.
\xmin = "33.5"; \xmax = "35.5"; \xinc = "0.05";
\ymin = "-2.0"; \ymax = "11.0"; \yinc = "0.25";
\infile = "tmp-xy.\.pid."
\outfile = "tmp-grid.\.pid."
system perl <<"EOF"
#
# Prepare 2 dimensional probability-density-function
# dataset for contouring by Gri. This reads (x,y)
# data from a file called $infile (defined below)
# and creates an output file called $outfile
# (also defined below) containing the
# x-grid, the y-grid, and then the grid data,
# suitable for reading/contouring by the Gri
# commands
# open tmp-grid.\.pid.
# read .number_of_data.
# read grid x
# read grid y
# read grid data
# draw contour
#
# The values in the output grid are defined so
# that they sum to the reciprocal of the
# product of the x binsize and y binsize (see
# definition of $factor below).
#
$xmin = \xmin; $xmax = \xmax; $xinc = \xinc;
$ymin = \ymin; $ymax = \ymax; $yinc = \yinc;
$infile = "\infile";
$outfile = "\outfile";
#
# Slurp in the x[], y[] data, storing
# the total number of data in $n.
open(IN, "$infile") || die "Can't open infile";
open(OUT, ">$outfile") || die "Can't open outfile";
$n = 0;
while(<IN>) {
($x[$n], $y[$n]) = split;
$n++;
}
#
# Zero out matrix (stored in a linear array scanned
# as one reads a book).
$cols = int(1 + ($xmax - $xmin) / $xinc);
$rows = int(1 + ($ymax - $ymin) / $yinc);
for ($y = $ymax; $y > $ymin; $y -= $yinc) {
for ($x = $xmin; $x < $xmax; $x += $xinc) {
$l = int(($x - $xmin) / $xinc
+ $cols * int(($y - $ymin) / $yinc));
$sum[$l] = 0;
}
}
#
# Cumulate (x,y) data into the matrix
$inside = 0;
for ($i = 0; $i < $n; $i++) {
if ($ymin <= $y[$i] && $y[$i] <= $ymax
&& $xmin <= $x[$i] && $x[$i] <= $xmax) {
$l = int(($x[$i] - $xmin) / $xinc
+ $cols * int(($y[$i] - $ymin) / $yinc));
$sum[$l]++;
$inside++;
} else {
print STDERR "($y[$i], $x[$i]) clipped\n";
}
}
#
# Print number of points (to allow renormalization
# if the user wishes)
print OUT "$inside\n";
#
# Print x grid, y grid, then grid itself
for ($x = $xmin; $x < $xmax; $x += $xinc) {
printf OUT "%lg\n", $x;
}
print OUT "\n";
for ($y = $ymax; $y > $ymin; $y -= $yinc) {
printf OUT "%lg\n", $y;
}
print OUT "\n";
#
# The $factor variable scales the PDF.
$factor = 1 / $xinc / $yinc / $inside;
for ($y = $ymax; $y > $ymin; $y -= $yinc) {
for ($x = $xmin; $x < $xmax; $x += $xinc) {
$l = int(($x - $xmin) / $xinc
+ $cols * int(($y - $ymin) / $yinc));
printf OUT "%lg ", $factor * $sum[$l];
}
print OUT "\n";
}
EOF
# Axes
set x margin 3
set x margin 6
set x name "Salinity [PSU]"
set y name "Potential Temperature [$\circ$C]"
set x axis 34.5 34.8 0.1
set y axis 5 9 1
draw axes
# PDF
open tmp-grid.\.pid.
read .number_of_data.
read grid x
read grid y
read grid data
.smooth. = 0
# Contours. Use clipping, since the axes are
# zooming in on part of the PDF.
set font size 8
set contour label position centered
set clip postscript on
set line width rapidograph 4x0
draw contour 0.2 2.2 0.4 unlabelled
set line width rapidograph 0
draw contour 0.4 2.4 0.4
set clip postscript off
end if
