#!/bin/bash

###################################################
# 
# Shell script designed to use a trench distance grids and trench-geometry files
# to step along a margin and create closed ploygon files of the small, intermediate
# and large masks required for modelling of subduction zone outer-rise and 
# seaward trench-slope gradient.
# 
# 
# 
# 
###################################################
##### Paramaters #####
###################################################


#	Switch GMT version and parameters manually
	gmtswitch /opt/local/lib/gmt4
	gmtset ELLIPSOID = WGS-84

margin=$1 # margin you are analysing

trench_distance_grid="/radar3/flex/trench_distance_grids/${margin}.grd" # Path to trench-distance grids
trench="/radar3/flex/trench_geometries/${margin}.xy" # Path to trench- geometry files 

echo $trench
echo $trench_distance_grid

# Specify the size of a-box
side_a="1024" # specify in # of minutes

# Specify the size of b-box in km (-ve is up-dip (seaward) of the trench-axis, +ve is down-dip (landward)
up_b="-800"
down_b="0"
width_b="800"

# Specify the size of c-box in km (-ve is up-dip (seaward) of the trench-axis, +ve is down-dip (landward)
up_c="-600" # distance up-dip of bathymetrically defined trench-axis
down_c="-10" # distance down-dip of bathymetrically defined trench-axis
width_c="300" # along strike extent

###################################################
###################################################
###################################################
# Calculate distance along the trench-axis

points=$(wc $trench | awk '{print $1}')
awk '{if(NR==1){print $1,$2,0}}' $trench > trench_distance.txt
for (( p=2 ; p <= $points ; p++ )); do
start=$(awk -v p=$p '{if(NR==(p-1)){print $1"/"$2}}' $trench)
stop=$(awk -v p=$p '{if(NR==p){print $1"/"$2}}' $trench)
dist=$(project -C$start -E$stop -G100 -Q | tail -n1 | awk -v d=$d2 '{print $3+d}')
d2=$dist
echo $stop | awk -F"/" -v d=$dist '{print $1,$2,d}' >> trench_distance.txt
done
sample1d trench_distance.txt -Fl -T2 -I1 > trench_distance_fine.txt
sample1d trench_distance.txt -Fl -T2 -I10k > td_10km.txt

###################################################
# Extract contours for down-dip limit of c-mask

###### C-mask ######
if [ $down_c -ne 0 ]; then
cpt_command=$(echo $down_c | awk '{print "-T"$1-1"/"$1+1"/"1}')
makecpt -Chaxby $cpt_command > contour.cpt
grdcontour $trench_distance_grid -Ccontour.cpt -Jx0.5 -A- -m -D > junk.ps
awk -v dn=$down_c '{if($3==dn){print $1,$2,NR}}' contour | sample1d -Fl -T2 -I0.2 | awk '{print $1,$2}' > c_down.txt
else
cp td_10km.txt c_down.txt
fi

###################################################

# Determine location of c-mask mid-points
pi=$(echo "scale=10; 4*a(1)" | bc -l)

# set "marching increment"
inc=$(echo $width_c | awk '{print $1/2}')
sample1d trench_distance.txt -Fl -T2 -I$inc > junk
sample1d trench_distance.txt -Fl -T2 -I1 > td_km.txt
points=$(wc junk | awk '{print $1}')


# LOOP OVER THE TRENCH AXIS STARTS HERE

for (( p=3 ; p < $points ; p++ )); do
# for (( p=20 ; p < 25 ; p++ )); do
mid=$(awk -v p=$p '{if(NR==p){print $1"/"$2}}' junk)
start=$(awk -v p=$p '{if(NR==(p-1)){print $1"/"$2}}' junk)
stop=$(awk -v p=$p '{if(NR==(p+1)){print $1"/"$2}}' junk)
start_td=$(awk -v p=$p '{if(NR==(p-1)){print $3}}' junk)
mid_td=$(awk -v p=$p '{if(NR==p){printf "%.0f\n", $3}}' junk)
stop_td=$(awk -v p=$p '{if(NR==(p+1)){print $3}}' junk)
X_dev=$(echo $start $stop | awk -F"/" '{print $1,$2,$3,$4}' | awk '{print $3-$1}')
Y_dev=$(echo $start $stop | awk -F"/" '{print $1,$2,$3,$4}' | awk '{print $4-$2}')


echo 25 > temp
Trig_ratio=$(gawk -v Xd=$X_dev -v Yd=$Y_dev '{print Xd/Yd}' temp)
Angle_rads=$(echo "scale=5; a($Trig_ratio)" | bc -l)
Strike=$(gawk -v AR=$Angle_rads -v PI=$pi -v YD=$Y_dev -v XD=$X_dev '{if(XD>=0 && YD<0) {print 180+(AR*(180/PI))} else if(XD>=0 && YD>=0){print AR*(180/PI)}else if(XD<0 && YD>0){print 360+(AR*(180/PI))}else if (XD<0 && YD<0){print (360-(AR*(180/PI)))}}' temp)
dip_direction=$(echo $Strike | awk '{if($1>270){print $1-90}else{print $1+90}}')

# echo X_dev $X_dev Y_dev $Y_dev Strike $Strike dip_direction $dip_direction
# echo start $start stop $stop start_td $start_td  $stop_td

###################################################
##### Determine down-dip direction #####
flip=$(project -C$mid -A$dip_direction -L-20/20 -G1 -Q | awk '{if($1<0){print 360+$1,$2,$3}else{print $1,$2,$3}}' | grdtrack -G$trench_distance_grid | awk '{if(NR==1 && $4>0){print 1}else if(NR==1 && $4<0){print 2}}')

if [ $flip -eq 1 ]; then
dip_direction=$(echo $dip_direction | awk '{if($1<=180){print $1+180}else{print $1-180}}')
fi

strike=$( echo $dip_direction | awk '{if($1>=45 && $1<=135){print 2}else if($1>=215 && $1<305){print 2}else{print 1}}')

###################################################
##### Construct C-Box #####

project -C$start -A$dip_direction -L$up_c/500 -G1 -Q | awk '{if($1<0){print 360+$1,$2,$3}else{print $0}}' > c_box_junk
# Determine intersection with trench-distance contour
grdtrack c_box_junk -G$trench_distance_grid | awk -v dw=$down_c '{print $1,$2,$3,(($4-dw)^2)^0.5}' > junky
tx=$(sort -n -k4 junky | head -n1 | awk '{print $1}')
ty=$(sort -n -k4 junky | head -n1 | awk '{print $2}')
td=$(sort -n -k4 junky | head -n1 | awk '{print $3}')
awk -v td=$td '{if($3<=td){print $0}}' c_box_junk > c_box


project -C$stop -A$dip_direction -L$up_c/500 -G1 -Q | awk '{if($1<0){print 360+$1,$2,$3}else{print $0}}' > c_box_junk
# Determine intersection with trench-distance contour
grdtrack c_box_junk -G$trench_distance_grid | awk -v dw=$down_c '{print $1,$2,$3,(($4-dw)^2)^0.5}' > junky
px=$(sort -n -k4 junky | head -n1 | awk '{print $1}')
py=$(sort -n -k4 junky | head -n1 | awk '{print $2}')
pd=$(sort -n -k4 junky | head -n1 | awk '{print $3}')



rm tsegs_int.dat
touch tsegs_int.dat

if [ $strike -eq "1" ]; then
awk -v tx=$tx -v px=$px '{if($1>tx && $1<px){print $1,$2}}' c_down.txt | sort -n >> c_box
awk -v tx=$tx -v px=$px '{if($1>tx && $1<px){print $1,$2}}' td_10km.txt | sort -n >> tsegs_int.dat
else
awk -v ty=$ty -v py=$py '{if($2>ty && $2<py){print $1,$2}}' c_down.txt  >> c_box
awk -v ty=$ty -v py=$py '{if($2>ty && $2<py){print $1,$2}}' td_10km.txt >> tsegs_int.dat
fi



if [ $flip -eq 1 ]; then
awk -v md=$pd '{if($3<=md){print $0}}' c_box_junk | sort -n -k2  -r >> c_box
else
awk -v md=$pd '{if($3<=md){print $0}}' c_box_junk | sort -n -k2 >> c_box
fi

head -n1 c_box >> c_box


###################################################
##### Construct B-Box #####
wid=$(echo $width_b | awk '{printf "%.0f\n", $1/2}')
b_start=$(awk -v wd=$wid -v p=$mid_td '{if($3==(p-wd)){print $1"/"$2}}' td_km.txt)
b_stop=$(awk -v wd=$wid -v p=$mid_td '{if($3==(p+wd)){print $1"/"$2}}' td_km.txt)
start_td=$(awk -v wd=$wid -v p=$mid_td '{if($3==(p-wd)){print $3}}' td_km.txt)
stop_td=$(awk -v wd=$wid -v p=$mid_td '{if($3==(p+wd)){print $3}}' td_km.txt)

echo b_start $b_start b_stop $b_stop
echo mid_td $mid_td start_td $start_td stop_td $stop_td

project -C$b_start -A$dip_direction -L$up_b/$down_b -G1 -Q -Dg > b_box


rm tsegs.dat
touch tsegs.dat

if [ $down_b -eq 0 ]; then
awk -v st=$start_td -v sp=$stop_td '{if($3>st && $3<sp){print $1,$2,$3}}' td_10km.txt >> b_box
fi

awk -v st=$start_td -v sp=$stop_td '{if($3>=st && $3<=sp){print $1,$2}}' td_10km.txt >> tsegs.dat

if [ $flip -eq 1 ]; then
project -C$b_stop -A$dip_direction -L$up_b/$down_b -G1 -Q -Dg | sort -n -k2 -r >> b_box
else
project -C$b_stop -A$dip_direction -L$up_b/$down_b -G1 -Q -Dg | sort -n -k2 >> b_box
fi
project -C$b_start -A$dip_direction -L$up_b/$down_b -G1 -Q -Dg | head -n1 >> b_box

###################################################
##### Construct A-Box #####
a_x=$(echo $mid | awk -F"/" '{printf "%0.1f\n", $1}')
a_y=$(echo $mid | awk -F"/" '{printf "%0.1f\n", $2}')

min_x=$(echo $a_x $side_a | awk '{print $1-(($2/2)/60)}')
max_x=$(echo $a_x $side_a | awk '{print $1+(($2/2)/60)}')
min_y=$(echo $a_y $side_a | awk '{print $1-(($2/2)/60)}')
max_y=$(echo $a_y $side_a | awk '{print $1+(($2/2)/60)}')

echo $min_x $min_y > a_box
echo $max_x $min_y >> a_box
echo $max_x $max_y >> a_box
echo $min_x $max_y >> a_box
echo $min_x $min_y >> a_box

###################################################
####################### PLOT ######################
###################################################

map="map.ps"

page_height="21"
page_width="28"

gmtset X_ORIGIN = 0
gmtset Y_ORIGIN = 0
gmtset PAPER_MEDIA = A4
gmtset ANNOT_FONT_SIZE_PRIMARY = 5p
gmtset ANNOT_FONT_SIZE_SECONDARY = 5p
gmtset LABEL_FONT_SIZE = 5p
gmtset HEADER_FONT_SIZE = 5p
gmtset COLOR_BACKGROUND = white
gmtset COLOR_FOREGROUND = white
gmtset COLOR_NAN = white

rangex=$(echo $min_x $max_x | awk '{print ($2-$1)-3}')
rangey=$(echo $min_y $max_y | awk '{print ($2-$1)-2}')

if (( $(echo "$rangex > $rangey" | bc -l) )); then
gmtset PAGE_ORIENTATION = LANDSCAPE
pj=$(echo 1 | awk -v ph=$page_height -v pw=$page_width -v rx=$rangex -v ry=$rangey '{if(((ph-4)/rx)*ry>(pw-4)){print(pw-4)/ry}else{print (ph-4)/rx}}')   
proj="-Jx"$pj"d"
else
gmtset PAGE_ORIENTATION = PORTRAIT
pj=$(echo 1 | awk -v ph=$page_height -v pw=$page_width -v rx=$rangex -v ry=$rangey '{if(((ph-4)/ry)*rx>(pw-4)){print(pw-4)/rx}else{print (ph-4)/ry}}')   
proj="-Jx"$pj"d"
fi

region=$(echo $min_x $max_x $min_y $max_y | awk '{print "-R"$1-1"/"$2+1"/"$3-1"/"$4+1}')

pscoast -X1 -Y1 $region $proj -W1 -Df -K -B5f1::/a5f1::neWS > $map
psxy $trench $region $proj -W5t10_20:2/150/150/150 -O -K >> $map
psxy junk $region $proj -Sp0.1 -Gred -W1/255/0/0 -O -K >> $map
psxy c_down.txt $region $proj -W1t10_20:1/150/150/150 -O -K >> $map
# psxy a $region $proj -W3/0/0/0 -O -K >> $map
psxy c_box $region $proj -W1/255/0/0 -O -K >> $map
psxy b_box $region $proj -W1/255/0/0 -O -K >> $map
psxy a_box $region $proj -W1/255/0/0 -O >> $map
gv $map 
# exit

let q=p-2

segno=$(printf "%02d" $q)

mv a_box a_box_${segno}
mv b_box b_box_${segno}
mv c_box c_box_${segno}
mv tsegs.dat tsegs_${segno}.dat
mv tsegs_int.dat tsegs_int_${segno}.dat
mv map.ps map_${segno}.ps

done

##### Clean up #####

rm c_down.txt junk temp tsegs.dat td_km.txt c_box_junk junky
# rm trench_distance.txt trench_distance_fine.txt c_down.txt junk temp td_km.txt c_box_junk junky

#	Switch GMT version and parameters manually
	gmtswitch /Users/esg006/gmt5