#!/bin/csh
# First run fircircle and get L2 average position and N hemisphere pole
#
fitcircle sat.xyg -L2 >! report
set cpos = `grep "L2 Average Position" report | awk '{ print $(NF-1), $NF}'` 
set ppos = `grep "L2 N Hemisphere" report | awk '{print $(NF-1), $NF}'` 
#
# Now you will need to project the data to find out what their ranges are in the projected 
# coordinate:
#
project  sat.xyg -C$cpos[1]/$cpos[2] -T$ppos[1]/$ppos[2] -S -V -Fpz -M >! sat.pg
project ship.xyg -C$cpos[1]/$cpos[2] -T$ppos[1]/$ppos[2] -S -V -Fpz -M >! ship.pg
#
# The minmax utility will report the minimum and maximum values for all columns. 
# Use this information to find the range of the projected distance coordinate they 
# have in common:
#
minmax sat.pg ship.pg
echo "We need the minimum and maximum distance the data sets have in common."
echo -n "Enter min max to nearest km (here -1167 1169): "
set line = $<
set limits = ($line)
#
# Make a file of equidistant sampling points spaced 1 km apart over the common interval.
# You will call this file samp.x  Here is one way to do it:
#
awk 'BEGIN { for (i = '$limits[1]'; i <= '$limits[2]'; i++) print i }' /dev/null >! samp.x
#
# Now you can resample the projected data:
#
sample1d sat.pg -Nsamp.x >! samp_sat.pg
sample1d ship.pg -Nsamp.x >! samp_ship.pg
#
# Now to do the cross-spectra, assuming that the ship is the input and the sat is the output 
# data, do this:
# 
paste samp_ship.pg samp_sat.pg | cut -f2,4 | spectrum1d -S256 -D1 -V -W -C
# 
# Now you want to plot the spectra.  The following commands will plot the ship and sat 
# power in one diagram and the coherency on another diagram,  both on the same page.  
# Note the extended use of pstext and psxy to put labels and legends directly on the plots.  
# For that purpose we often use -Jx1 and specify positions in inches directly:

psxy spectrum.coh -Ba1f3p:"Wavelength (km)":/a0.25f0.05:"Coherency@+2@+":WeSn -JX-4l/3.75 -R1/1000/0/1 -U/-2.25/-1.25/"Example 3 in Cookbook" -P -K -X2.5 -Sc0.07 -G0 -Ey/2 -Y1.5 >! example_03.ps
echo "3.85 3.6 18 0.0 1 11 Coherency@+2@+" | pstext -R0/4/0/3.75 -Jx1 -O -K >> example_03.ps
cat << END >! box.d
2.375   3.75
2.375   3.25
4       3.25
END
psxy -R -Jx1 -O -K -W5 box.d >> example_03.ps
psxy -St0.07 -O -Ba1f3p/a1f3p:"Power (mGal@+2@+km)"::."Ship and Satellite Gravity":WeSn spectrum.xpower -R1/1000/0.1/10000 -JX-4l/3.75l -Y4.2 -K -Ey/2 >> example_03.ps
psxy spectrum.ypower -R -JX -O -K -G0 -Sc0.07 -Ey/2 >> example_03.ps
echo "3.9 3.6 18 0.0 1 11 Input Power" | pstext -R0/4/0/3.75 -Jx1 -O -K >> example_03.ps
psxy -R -Jx -O -K -W5 box.d >> example_03.ps
psxy -R -Jx -O -K -G240 -L -W5 << END >> example_03.ps
0.25    0.25
1.4     0.25
1.4     0.9
0.25    0.9
END
echo "0.4 0.7" | psxy -R -Jx -O -K -St0.07 -G0 >> example_03.ps
echo "0.5 0.7 14 0.0 1 5 Ship" | pstext -R -Jx -O -K >> example_03.ps
echo "0.4 0.4" | psxy -R -Jx -O -K -Sc0.07 -G0 >> example_03.ps
echo "0.5 0.4 14 0.0 1 5 Satellite" | pstext -R -Jx -O >> example_03.ps
\rm -f box.d report samp.x *.pg spectrum.* .gmtcommands