# Examplex 4 -- TS diagram (Colour)
#
# Contributed 17 Feb 1994 by Dan Kelley <Dan.Kelley@Dal.Ca>.


#
# On a light brown background, draw blue [T(t),S(t)] evolution
# diagrams for various years, with data from Labrador Sea, and
# a predicted slope in red.  Isopycnals are drawn in black. 
# (This is a color overhead for seminar, hence the fontsize etc)

# Define colors to use
\TS_trace = "blue"
\isopycnals = "black"
\background = "rgb 1.0000 0.9059 0.7294"
\axes = "black"
\guiding_line = "red"

`Set Up'
{
    \data = "./bravo1.25db"
    .dSdT. = -0.011             # Slope for reference line
    .start. = 8                 # Start month (wrt Jan=0; need negative)
    .stop. = 12                 # Stop  month
    .delta_t. = 1               # Interpolation time interval (days)
    \filter = "2_0.02"          # Butterworth filter file
    .first. = 1964              # First year to plot
    .last. = 1967               # Last year to plot
    set missing value -9
    set x margin 3
    set y margin 6
    set x name "Salinity, PSU"
    set y format %.0lf$\circ$
    set y name "Temperature, $\circ$C"
    set x axis 34.5 34.8 0.1
    set y axis 3 9 2 1
    set line width axis rapidograph 4x0
    set color \background
    draw polygon filled         \
        ..xleft.. ..ybottom..   \
        ..xright.. ..ybottom..  \
        ..xright.. ..ytop..     \
        ..xleft.. ..ytop..
    set color \axes
    set font size 16            # For overhead projection
    draw axes
    system rm -f tmp*[0-9]
}

`Draw Isopycnals'
{
    set clip on
    set line width rapidograph 1
    set color \isopycnals
    draw isopycnal 26.75 unlabelled
    draw isopycnal 27.00 unlabelled
    draw isopycnal 27.25 unlabelled
    draw isopycnal 27.50 unlabelled
    # Label manually, for nicer effect (laborious)
    set font size 14
    draw label "1027.00" at 34.731690 8.522833 rotated 16
    #draw label "1027.25" at 34.731690 6.808333 rotated 17
    draw label "1027.50" at 34.731690 4.797500 rotated 19
    set clip off
}

`Draw guiding line'
{
    set clip postscript on
    set color \guiding_line
    set line width rapidograph 1
    .S. = 34.7
    .T. = -1.8
    draw line from .S. .T. to {rpn .S. ..ytop.. .T. - .dSdT. * +} ..ytop..
    set dash off
    set clip postscript off
}

`Interpolate to standard times .days. \outfile_raw \outfile_smoothed'
#
#.days.   = sampling interval
#
{
    # Trim out the (alternate) header lines, bin in time intervals
    # (0.0208333year = 1wk; 0.0416667yr=2week), then put 
    # into file 
    new \time2 color=#CDAD00> .days.
    .days. = \.word4.
    show "Interpolating `\data' to " .days. " days"
    show "Using butterworth filter ~/filters/butterworth/\filter"
    sprintf \time2 "%.10lf" {rpn 1964 .days. 365.25 / +}
    system                                                      \
        cat \data                                               \
        | gawk ' {                                              \
        if (NR % 2)                                             \
        printf("%lf ", $1 + $2 / 365.0);                        \
        if (1 - NR % 2)                                         \
        printf("%lf %lf\n",$3,$2)                               \
        }'                                                      \
        | interpolate3col 1964 \time2 1974                      \
        | cat > \.word5.
    if ..exit_status..
        show "System call failed"
        quit
    end if
    system                                                      \
        cat \.word5.                                            \
        | butterfilter -c 3 $HOME/filters/butterworth/\filter   \
        | cat > \.word6.
    if ..exit_status..
        show "System call failed"
        quit
    end if
    delete .days. \time2
}

`Select years .start. .stop. \input_file \output_file'
{
    system gawk '\.word2. <= $1 && $1 <= \.word3.' < \.word4. > \.word5.
}

`Draw year labels'
{
    # Draw labels for curves manually, to avoid overwriting
    set color \TS_trace
    set font size 14
    draw label whiteunder "1964" at 34.643232 \
        {rpn 8.332266 yusertocm "M" ascent 2 / - ycmtouser}
    draw label whiteunder "1965" at 34.615939 \
        {rpn 7.463851 yusertocm "M" ascent 2 / - ycmtouser}
    draw label whiteunder "1966" at 34.653777 \
        {rpn 6.930395 yusertocm "M" ascent 2 / - ycmtouser}
    draw label whiteunder "1967" at                             \
        {rpn 34.545845 xusertocm "1967" width 2 / - xcmtouser}  \
        8.865722
}

# Main program

set line width axis rapidograph 00
Set Up
#Draw Isopycnals
Draw guiding line
Interpolate to standard times .delta_t. tmp1.\.pid. tmp2.\.pid.

# Plot timeseries
set color \TS_trace
set line width rapidograph 2
.year. = .first.
while {rpn .last. .year. <=}
    .t0. = {rpn .year. .start. 12 / +} 
    .t1. = {rpn .year.  .stop. 12 / +}
    show "Plotting " .t0. "-" .t1.
    Select years {rpn .t0.} {rpn .t1.} tmp2.\.pid. tmp3.\.pid.
    open tmp3.\.pid.
    read columns * x y
    close
    draw curve overlying
    .year. += 1
end while
set color \TS_trace
Draw year labels

# Draw label indicating months
set line width rapidograph 00
\label = "Aug-Dec"
set color \TS_trace
draw label boxed "\label" at                                            \
    {rpn ..xleft.. ..xright.. + 2 / xusertocm "\label" width 2 / -}     \
    {rpn ..ytop.. yusertocm "M" ascent 2.5 * +}                         \
    cm

.offset. = 3                    # distance of labels above plot
set color \axes
draw label "$OWS Bravo$,  Labrador Sea"                 \
    centered at                                         \
    {rpn ..xleft.. ..xright.. + 2 / xusertocm}          \
    {rpn ..ymargin.. ..ysize.. + .offset. + "M" ascent 6 * +} cm

system rm -f tmp*.\.pid.

quit