#!/bin/sh
#
# drft06rr_dredgemap will take a dredge ID from the command line and make
# a contour map of bathymetry at the dredge site with the dredge
# track overlaid.  If 2-4 IDs are provided on the command line
# we will instead make a 4-up page with up to 4 dredge maps at
# 50% (for cruise report).  Needs drft06rr_dredgetime.

# Set the width of the box in arc minutes:
BOXWIDTH=2
TL=0.2
#---------------------------------------------------------
cd dredges
gmtset OUTPUT_DEGREE_FORMAT D
W=`gmtmath -Q $BOXWIDTH 2 DIV 60 DIV 6371.008 PI MUL 2 MUL 360 DIV MUL =`

if [ $# -eq 0 ]; then
	echo "usage: drft06rr_dredgemap ID [ID2 ID3 ID4] [-f]"
	exit
fi

if [ $# -gt 1 ] && [ ! $2 = "-f" ]; then
	multi=1
	PS=dredge_chart.ps
	psxy -R0/15/0/15 -Jx1 -P -X1.25 -Y2 -K /dev/null > $PS
	echo "0.5 dup scale" >> $PS
	Wc=-Wc0.5p
	Wa=-Wa1p
else
	multi=0
	PS=drft06rr.D-${1}.map.ps
	EPS=drft06rr.D-${1}.map.eps
	GIF=../drft06rr/dredges/drft06rr.D-${1}.map.gif
	if [ $# -eq 1 ] && [ -f $GIF ]; then
		echo "drft06rr_dredgemap: Map already exists, use -f to force a new map" >&2
		exit
	fi
	PDF=../drft06rr/dredges/drft06rr.D-${1}.map.pdf
	psxy -R0/6/0/7 -Jx1 -P -K -X1.25 -Y2 /dev/null > $PS
	Wc=-Wc0.25p
	Wa=-Wa0.75p
fi
i=0
x=0
y=0

for ID in $*; do
	if [ $ID = "-f" ]; then
		continue
	fi
	if [ ! -f $ID.xyt ]; then
		echo "dredgemap: Cannot find file $ID.xyt"
		exit
	fi
	
	i=`expr $i + 1`
	if [ $multi -eq 1 ]; then
		if [ $i -eq 1 ]; then
			x=0
			y=7
		elif [ $i -eq 2 ]; then
			x=6.75
			y=0
		elif [ $i -eq 3 ]; then
			x=-6.75
			y=-7
		else
			x=6.75
			y=0
		fi
	fi
	gmtconvert $ID.xyt -fi0x,1y,2T > tmp
	x0=`head -1 tmp | cut -f1`
	y0=`head -1 tmp | cut -f2`
	x1=`tail -1 tmp | cut -f1`
	y1=`tail -1 tmp | cut -f2`
	az=`gmtmath -Q $y1 $y0 SUB $x1 $x0 SUB DIV ATAN R2D 90 ADD DUP 90 GT 180 MUL SUB =`
	DX=`gmtmath -Q $az COSD $TL MUL =`
	DY=`gmtmath -Q $az SIND $TL MUL =`

# Find approximate location of dredge using mean dredge location

	gmtmath -S -Ca $ID.xyt -f0x,1y,2T MEAN = tmp
	lon=`cut -f1 tmp`
	lat=`cut -f2 tmp`
	rm -f tmp

# Find which bathymetry box this dredge fell in (A-G)

	grep -v '^#' ../drft06rr_boxes.lis | awk '{printf "%lg\t%lg\t%s\n", 0.5 * ($2 + $3), 0.5 * ($4 + $5), $1}' > tmp
	box=`mapproject -G$lon/$lat/k tmp -f0x,1y | sort -k 4,4 -g | awk '{if (NR == 1) print $3}'`

# Find -R for a small square of 2*$W km centered on the mean location

	range=`(echo -$W -$W; echo $W $W) | mapproject -JE$lon/$lat/1 -R0/360/-90/90 -I -Fk -C`
	R=`echo $range | awk '{printf "-R%f/%f/%f/%f\n",$1,$3,$2,$4}'`
	west=`echo $range | awk '{print $1}'`
	north=`echo $range | awk '{print $4}'`
	H=`echo $west $north | mapproject $R -JM | cut -f2`
	H=`gmtmath -Q $H 0.05 SUB =`
	HS=`gmtmath -Q $H 0.4 SUB =`
	t_on=`head -1 $ID.xyt | cut -f3`
	t_off=`tail -1 $ID.xyt | cut -f3`
# Lay down background bathymetry contours, first pull out subset and grdmath away the SS

	grdmath ../gridding/grids/drft06rr.box$box.mb+ss.0.001.grd 1 FMOD ABS 0.02 GT 0 NAN ../gridding/grids/drft06rr.box$box.mb+ss.0.001.grd MUL = mb.grd
	grdcontour $R mb.grd -G6 -JM6 -O -Z-1 -Nm -A100 -C25 -K -S16 -T0.15/0.03 $Wa $Wc-X$x -Y$y >> $PS
	psxy -R -JM -O -K $ID.xyt -Sc0.025 -G0 >> $PS
	head -1 $ID.xyt | psxy -R -JM -O -K -Sc0.08 -W0.5p -G0 >> $PS
	tail -1 $ID.xyt | psxy -R -JM -O -K -Sc0.08 -W0.5p -G255 >> $PS
	../drft06rr_dredgetime $ID | awk '{if ($NF < 0.02) print $0}' | psxy -R -JM -O -K -SA0.08 -W0.5p -G255 >> $PS
	../drft06rr_dredgetime $ID | mapproject -R -JM | awk '{printf ">\n%s\t%s\n%lg\t%lg\n", $1, $2, $1+'$DX', $2+'$DY'}' | psxy -R0/6/0/6 -Jx1 -O -K -M -W0.5p >> $PS
	../drft06rr_dredgetime $ID | awk '{printf "%s %s 10 %lg 0 LM %s:%s\n", $1, $2, '$az', $3, $4}' | pstext $R -JM -O -K -W255o0.5p -C0.02/0.02 -D$DX/$DY >> $PS
	gmtset ANOT_FONT_SIZE 9 LABEL_FONT_SIZE 10
	pslegend -R -JM -Dx0.05/$H/1.8/1.2/LT -O -K -F -G255 << EOF >> $PS
G -0.05i
H 12 5 Drift Leg 6 D-$ID
G 0.5i
S 0.05 c 0.1i 0/0/0 0.5p 0.2 On bottom: $t_on
S 0.05 c 0.1i 255/255/255 0.5p 0.2 Off bottom: $t_off
S 0.05 a 0.125i 255/255/255 0.5p 0.2 Ship is stationary
EOF
	psbasemap $R -JM -O -K -B1mf10cWSne -Lfx0.95/$HS/$lat/0.5 >> $PS
	gmtset ANOT_FONT_SIZE +12p
done
psxy -R -JM -O /dev/null >> $PS

if [ $multi -eq 0 ]; then
	ps2eps $PS 2> /dev/null
	mv -f $EPS ../drft06rr/dredges/
	convert -density 100x100 -crop 0x0 $PS $GIF
	ps2pdf $PS $PDF
else
	gv $PS &
fi
rm -f tmp mb.grd
cd ..