\cancel mode verify
! bar_chart2.jnl
! 4/94
! bug fixes 8/94: use "ppl yaxis" to set lower plot bound
!		  use "ppl xaxis" to make X axis wide enough
! bug fix  12.98: overlay/vs -> overlay>overlay/vs *jd* per *sh*

! Description: Make a bar chart using hollow rectangles

! Usage:

!	GO bar_chart2 plot_command expression [axis] ["overlay"] [lower_bound]
! 1 plot_command   - "PLOT" with its qualifiers
! 2 expression     - a line of data (oriented along "axis") to be bar-charted
! 3 axis	   - "X", "Y", "Z", or "T"
! 4 overlay	   - "OVERLAY" or " " or "CLEAR"
! 5 lower_bound	   - the bottom of the bar chart (default = below min data)

! Example:
!	yes? SET DATA coads_climatology
!	yes? SET REGION/L=1/Y=40S:40N/X=110W
!	yes? GO  bar_chart2 PLOT/LINE=1 sst Y
!	yes? GO  bar_chart2 PLOT/LINE=2 sst[x=90w] Y overlay

! Note: To use a T axis requires CANCEL MODE CALENDAR

! Note: as of 4/25/94 this routine has a bug and the first bar is missing.
!       Kludge sol'n is to plot one extra bar at the beginning

! check for required args
QUERY/IGNORE $1%<Usage: GO bar_chart2 "PLOT"[/quals] expression [axis] ["OVERLAY"] [lower_bound]%
QUERY/IGNORE $2%<Usage: GO bar_chart2 "PLOT"[/quals] expression [axis] ["OVERLAY"] [[lower_bound]%
QUERY/IGNORE $3%X|X|Y|Z|T|<Usage: GO bar_chart2 "PLOT"[/quals] expression [axis] ["OVERLAY"] [[lower_bound]%
QUERY/IGNORE $4%clear|overlay|clear|<Usage: GO bar_chart2 "PLOT"[/quals] expression [axis] ["OVERLAY"] [[lower_bound]%

! determine the grid for bar_expr (so SET GRID works)
LET/QUIET bar_expr = missing($2,-9999.)	! attach a known missing data flag to the data
LOAD bar_expr
SET GRID/SAVE
SET GRID bar_expr

! determine the start and end of the X axis of the plot
LET/QUIET x_expr = $3"X|X|Y|Z|T"
LET/QUIET x_ax_end   = x_expr[$3"X|X|Y|Z|T"=@max]
LET/QUIET x_ax_start = x_expr[$3"X|X|Y|Z|T"=@min]
LET/QUIET x_ax_delta = (x_ax_end - x_ax_start)/(x_expr[$3"X|X|Y|Z|T"=@ngd]+x_expr[$3"X|X|Y|Z|T"=@nbd]-1)

! set the lower edge of the boxes
LET/QUIET bar_min = y3[X=@MIN]
LET/QUIET bar_y0 = $5"bar_min-0.2*ABS(bar_min)"   ! a bit below the min

! write the data to a file so it can be placed on a simple counting axis
! write the first and last X values to help complete the edge bars
sp rm -f bar_plot.dat
LIST/NOHEAD/FILE=bar_plot.dat/FORMAT=(2(1PG15.7))/APPEND x_ax_start - x_ax_delta, 0
LIST/NOHEAD/FILE=bar_plot.dat/FORMAT=(2(1PG15.7))/APPEND $3"X|X|Y|Z|T", bar_expr 
LIST/NOHEAD/FILE=bar_plot.dat/FORMAT=(2(1PG15.7))/APPEND x_ax_end+x_ax_delta, 0
LIST/NOHEAD/FILE=bar_plot.dat/FORMAT=(2(1PG15.7))/APPEND x_ax_end+2*x_ax_delta, 0

! prepare to read it back
SET DATA/SAVE
CANCEL DATA/NOERR bar_plot.dat
FILE/VARIABLES=bar_x,bar_y bar_plot.dat
SET DATA/RESTORE

! replicate the points x 3 by regridding to an axis from 1/3 to n+1/3 by 1/3
! then delete all except the integer point in each group of three
DEFINE AXIS/X=.66666667:10000:.3333333333    axtimes3  ! 2/3:big:1/3
DEFINE GRID/X=axtimes3 gtimes3
LET/QUIET bar_index = MOD(i[g=gtimes3]-1,3) + 1
LET/QUIET midmask = IF bar_index EQ 2 THEN 1
LET/QUIET x_mid   = midmask * bar_x[G=gtimes3,d=bar_plot.dat]
LET/QUIET y_mid   = midmask * bar_y[G=gtimes3,d=bar_plot.dat]
LET/QUIET x3   = x_mid[x=@FNR]
LET/QUIET y3_9999   = y_mid[x=@FNR]
LET/QUIET y3 = IF y3_9999 NE (-9999) THEN y3_9999 !screen out REAL missing data


! Define the region "bar" to be the number of points x3
LIST/NOHEAD/APPEND/FORMAT=("DEFINE REGION/I=1:",1PG11.6," bar")/FILE=bar_commands.jnl (bar_expr[$3"X|X|Y|Z|T"=@ngd]+1)*3 + 1
GO bar_commands.jnl
sp rm -f bar_commands.jnl     ! clean up

! define the pieces of each "bar" - 3 points
! for a data point i (Xi,Yi), the 3 points are 
! P1 = (Xi + Xi-1)/2, bar_y0
! P2 = (Xi + Xi-1)/2, Yi
! P3 = (Xi+1 + Xi)/2, Yi

LET/QUIET p1x = IF bar_index EQ 1 THEN (x3 + x3[x=@shf:-1])/2
LET/QUIET p2x = IF bar_index EQ 2 THEN (x3 + x3[x=@shf:-2])/2
LET/QUIET p3x = IF bar_index EQ 3 THEN (x3[x=@shf:+1] + x3)/2
LET/QUIET p1y = IF bar_index EQ 1 THEN bar_y0
LET/QUIET p2y = IF bar_index EQ 2 THEN y3
LET/QUIET p3y = IF bar_index EQ 3 THEN y3

! and now the fully assembled plot
LET/QUIET bar_px = MISSING(p1x,MISSING(p2x,p3x))
LET/QUIET bar_py = MISSING(p1y,MISSING(p2y,p3y))
SET VARIABLE/TITLE="Value" bar_px
SET VARIABLE/TITLE="Histogram Magnitude" bar_py

! plot it
DEFINE REGION/DEFAULT save
cancel region    ! probably unnecessary ...
$1/$4"vs|clear>vs|overlay>overlay/vs|"/@bar/D=bar_plot.dat bar_px,bar_py
SET GRID/SAVE
SET REGION save