#!/bin/sh

#####
#	this is a program to compute the geometrical spreading from the
#	velocity model for the specified source depth
#
#	the geometrical spreading of from Bullen.
#	I assume that MODEL is a spherical model
#
#	We follow Bullen and the Herrmann notes http://www.eas.slu.edu/People/RBHerrmann/Courses/EASA462/
#		for Ass14.pdf
#
#	The geometrical spreading is dimensionless and gives the decrease in amplitude from
#	a distance of 1 km from the source to the receiver
#                 2                            2
#            (   a sin Is Vs                  d T      )
#       sqrt |  ------------------------     -----     |
#            |                              2      2   |
#            (  sin DEG  cos Ir Rs Cos Is Rr  d DEL    )
#
#	a = radius of sphere about source - we use 1.0 km
#	Is= incident angle at source
#       Ir= incident angle at receiver
#       Rs= distance from center of Earth to source
#       Rr= distance from center of Earth to receiver
#       DEG=epicental distance in degrees
#       DEL=distance in radians
#
#	note that the second derivative divided by Rr^2 is sec/km^2 - which we can get
#       from dp/dx  where dp is in sec/km and dx is the distance in km, so we compute
#                 2                            2
#            (   a sin Is Vs                  d T      )
#       sqrt |  ------------------------     -----     |
#            |                                  2      |
#            (  sin DEG  cos Ir Rs Cos Is     dx       )
#
#
#	The only imperfection is the numerical computation of dp/dD
#
#	TO DO:
#		modify time96 to give the velocity at the source
#		modify time96 to give the velocity at the receiver
#		modify time96 to give the T* for the path
#		modify time96 to give everything for pP and sP
#####
HS=10

MODEL=tak135sph.mod
#####
#	compute the geometrical spreading which is dimensionless
#	it propagates the amplitude from a distance of 1,0 km from the source to
#	the receiver
#####
for PHASE in P pP sP
do
echo " "
echo "--------------------------------------------------"
echo "          Phase ${PHASE}              "
echo " DEG    HS    T(s)    p(s/deg)    Geom     T*(sec)"
echo "--------------------------------------------------"
for DEG in 30. 32.5 35. 37.5 40. 42.5 45. 50. 55. 60. 65. 70.  75. 80. 85. 90.
do
	#####
	# VS  P velocity at the source in km/sec
	# VR  P velocity at the receiver in km/sec
	# RR  Radius of the earth sinc ethe receiver is assumed to be at the surface
	# RS  Radius of the source which is RR - HS (computed)
	# DDEG increment in degrees for computing DP/DD by centered difference
	#####
	VS=6.2
	VR=6.2
	RR=6371
	RS=`echo $RR $HS | awk '{ print $1 - $2 }' `
	#####
	# compute the travel time and ray parameter and T*
	#####
	DDEG=5
	DEGP=`echo $DEG $DDEG | awk '{print $1 + $2  }' `
	DEGM=`echo $DEG $DDEG | awk '{print $1 - $2  }' `
	A=`time96    -${PHASE} -T    -M ${MODEL} -GCARC ${DEG} -EVDP $HS`
	RAYP=`time96 -${PHASE} -RAYP -M ${MODEL} -GCARC ${DEG} -EVDP $HS`
	TS=`time96 -${PHASE} -TS -M ${MODEL} -GCARC ${DEG} -EVDP $HS`
	#####
	# get the ray parameter in sec/deg for output
	#####
	RAYPDEG=`echo $RAYP | awk '{print $1 * 111.195 }' `
	#####
	# get ray parameter in sec/radian for computation of angles
	#####
	RAYPRAD=`echo $RAYPDEG | awk '{print $1 * 180.0/3.1415927 }' `
	#####
	#	 RAYP = sec/km
	#		sec/km * 111.195 km/deg * 180 deg/pi radians
	#####
	RAYPP=`time96 -${PHASE} -RAYP -M ${MODEL} -GCARC ${DEGP} -EVDP $HS`
	RAYPM=`time96 -${PHASE} -RAYP -M ${MODEL} -GCARC ${DEGM} -EVDP $HS`
	#####
	# get DP/Dkm by centered difference which is obtained from the two valeus of
	# p(sec/km) divided by the distance in km, hence the 111.195
	#####
	DPDKM=`echo $RAYPM $RAYPP $DDEG | awk '{print ( $1 - $2)/(2. * $3 * 111.195 ) }' `
	SINIS=`echo $VS $RS $RAYPRAD | awk '{print $3 * $1 / $2 }' `
	COSIS=`echo $SINIS | awk '{print sqrt ( 1.0 - $1 * $1 ) }' `
	SINIR=`echo $VS $RR $RAYPRAD | awk '{print $3 * $1 / $2 }' `
	COSIR=`echo $SINIR | awk '{print sqrt ( 1.0 - $1 * $1 ) }' `
	SIND=`echo $DEG | awk '{print sin ( $1 * 3.1415927 / 180.0 ) }' `
	
	GEOM2=`echo $SINIS $VS $DPDKM $SIND $COSIR $RS $COSIS | awk '{ print ($1 * $2 * $3 )/($4 * $5 *$6 *$7) }' `
	GEOM=`echo $GEOM2 | awk '{print sqrt ( $1 ) }' `

	echo $DEG $HS $A $RAYPDEG $GEOM $TS| awk '{printf "%5.1f %5.1f %8.2f %7.3f %12.4g %5.2f\n",$1, $2, $3, $4, $5, $6 }' 
done
done