1#!/usr/bin/env bash
2#
3# Compute FAA and VGG for case of variable density
4ps=variablerho.ps
5gmt set GMT_FFT kiss
6order=4
7dx=0.5
8# 2 panels of topo and grav, with top profile of admittance & coherence
9# NOT FINISHED
10# 1. Create a bathymetry data set with one circular plateau
11#    with R_base = 50 km, height = 3500 m, depth = -5000 m,
12#    and let density contrast to from 1300 at flanks and reach
13#    1700 at a radius of 25 km and stay at that through the center.
14#    The equivalent average density is 1533.33333
15rhoc=1533.333333
16echo "0	0	50	3500" | gmt grdseamount -R-256/256/-256/256 -I$dx -r -Cd -Gplateau.nc -Z-5000
17gmt grdmath -Rplateau.nc 50 0 0 CDIST SUB 25 DIV 0 MAX 400 MUL 1300 ADD 1700 MIN 0 0 CDIST 50 LT MUL = rho.nc
18# BL Plot the bathymetry
19gmt grdtrack -Gplateau.nc -ELM/RM -o0,2 -nn > plateau.trk
20gmt psxy -R-256/256/-5000/0 -JX5.5i/2.75i -P -K -Ggray -L+yb -BWSn -Baf -Bx+ukm -By+l"Depth (m)" plateau.trk -X1.5i > $ps
21echo "-256 0 BATHYMETRY" | gmt pstext -R -J -O -K -F+jTL+f14p,gray -Dj0.1i -Gwhite -C+tO >> $ps
22# Plot the variable density profile in red and constant in blue
23gmt grdtrack -Grho.nc -ELM/RM -o0,2 -nn > rho.trk
24gmt psxy -R-256/256/0/1800 -J -O -K -W1p,red rho.trk -BE -Baf -By+l"Density (kg/m@+3@+)" >> $ps
25gmt psxy -R -J -O -K -W1p,blue << EOF >> $ps
26-256	0
27-50	0
28-50	$rhoc
2950	$rhoc
3050	0
31256	0
32EOF
33echo "+256 1800 CONSTANT DENSITY" | gmt pstext -R -J -O -K -F+jTR+f14p,blue -Dj0.1i -Gwhite -C+tO >> $ps
34echo "+256 1800 VARIABLE DENSITY" | gmt pstext -R -J -O -K -F+jTR+f14p,red -Dj0.1i/0.4i -Gwhite -C+tO >> $ps
35# 2a. Compute the VGG anomaly for variable density
36gmt gravfft plateau.nc+uk -Drho.nc+uk -Nf+a -Fv -E$order -Gvgg.nc
37# BR plot the VGG anomaly
38gmt grdtrack -Gvgg.nc -ELM/RM -o0,2 -nn > vgg.trk
39gmt psxy -R-256/256/-100/250 -J -O -K -W1p,red vgg.trk -BWSne -Bafg1000 -Bx+ukm -By+l"VGG (E\371tv\371s)" -Y3.2i >> $ps
40# 2b. Compute the VGG anomaly for constant density
41gmt gravfft plateau.nc+uk -D$rhoc -Nf+a -Fv -E$order -Gvgg.nc
42# BR plot the VGG anomaly
43gmt grdtrack -Gvgg.nc -ELM/RM -o0,2 -nn > vgg.trk
44gmt psxy -R -J -O -K -W1p,blue vgg.trk >> $ps
45echo "-256 250 VGG" | gmt pstext -R -J -O -K -F+jTL+f14p -Dj0.1i -Gwhite -C+tO >> $ps
46echo "+256 250 CONSTANT DENSITY" | gmt pstext -R -J -O -K -F+jTR+f14p,blue -Dj0.1i -Gwhite -C+tO >> $ps
47echo "+256 250 VARIABLE DENSITY" | gmt pstext -R -J -O -K -F+jTR+f14p,red -Dj0.1i/0.4i -Gwhite -C+tO >> $ps
48# 3a. Compute the FAA anomaly for variable density
49gmt gravfft plateau.nc+uk -Drho.nc+uk -Nf+a -Ff -E$order -Gfaa.nc
50# ML plot the FAA anomaly
51gmt grdtrack -Gfaa.nc -ELM/RM -o0,2 -nn > faa.trk
52gmt psxy -R-256/256/-25/250 -J -O -K -W1p,red faa.trk -BWSne -Bafg1000 -Bx+ukm -By+l"FAA (mGal)" -Y3.2i >> $ps
53# 3a. Compute the FAA anomaly for constant density
54gmt gravfft plateau.nc+uk -D$rhoc -Nf+a -Ff -E$order -Gfaa.nc
55gmt grdmath faa.nc DUP LOWER SUB = faa.nc
56# ML plot the FAA anomaly
57gmt grdtrack -Gfaa.nc -ELM/RM -o0,2 -nn > faa.trk
58gmt psxy -R -J -O -K -W1p,blue faa.trk >> $ps
59echo "-256 250 FAA" | gmt pstext -R -J -O -K -F+jTL+f14p -Dj0.1i -Gwhite -C+tO >> $ps
60echo "+256 250 CONSTANT DENSITY" | gmt pstext -R -J -O -K -F+jTR+f14p,blue -Dj0.1i -Gwhite -C+tO >> $ps
61echo "+256 250 VARIABLE DENSITY" | gmt pstext -R -J -O -F+jTR+f14p,red -Dj0.1i/0.4i -Gwhite -C+tO >> $ps
62rm -f vgg.trk faa.trk plateau.trk rho.trk rho.nc plateau.nc vgg.nc faa.nc
63