#!/bin/sh

#####
#	THIS IS A SCRIPT TO CROSS CORRELATE SLU BB STATIONS
#	THESE STATIONS HAVE THE SAME SENSORS AND DATALOGGERS
#	SO I DO NOT DECONVOLVE THE INSTRUMENT
#####

#####
#	PROCESSING STEPS
#
#####

#####
#	Master parameters
#	for final display, just show cross-correlation at +- 300 sec
#####
WINL=-600
WINH=+600
#####
#	in correlating with the one hour data segments, divide into this number of
#	sub-segments. This means that the final time series will be about +-900 seconds
#	instead of the +- 3600 seconds for the cross-correlation. This will speed processing
#	Note however than no additional whitening is done on these subsegments
#####
SEGMENT=2
#####
#	define the freqlimitf for the whitening - this should be large enough
#	to include the expected signal
#####
FREQLIMITS="0.01 0.02 5. 8.0"
#####
#	High pass corner - this is to get rid of very low frequency terms. This
#	corresponding period should be less than the time series length and
#	greater than periods of interest
#####
HPCORNER="0.01"

#####
#	MYPWD: location of the directory in which this script is run.
#		this directory will have all of the stacked results
#####
MYPWD=`pwd`
#####
#	BASE:	location of the high level at which the raw SAC files
#		are stored. For this data set, the SAC files will
#		be in ${BASE}/2004/350, for example
#####
BASE=/rrd0/home/rbh/DATA/FLOR
#####
#	TEMP:	location of a temporary storage place
#	this contains a low-pass filtered/decimated version of the two traces 
#	in the trace pair to be processed
#
#	it will also contain the daily cross-corrleations and auto-correlations
#	for each hour segment in the form ${F1}.${HS}.xx ${F2}.${HS}.xy
#
#	it will also contain the final stacked cross-correlation for the day in the
#	form
#	2005.035.seoBHZsesbgz.rev
#	2005.035.seoBHZsesbgz.cor
#####
TEMP=${BASE}/TEMP
#####
#	create TEMP if it does not exist, clean it up if it does exist
#####
if [ ! -d ${TEMP} ]
then
	echo creating ${TEMP}
	mkdir ${TEMP}
else
	echo cleaning #{TEMP}
	rm -f  ${TEMP}/*
fi

#####
#	define the component loop values 
#	Note even though we can do a BHZ - bge or bgz - bge  we will only do
#	BHZ - bgz, bge - bge, bgn-bgn
#
#	CMPZ defines the vertical component so that we can define the interstation
#	distances and azimuths
#
#	COMP defines the list of components ot process
#####

CMPZ="BHZ"
COMP="BHZ BHR BHT"
OCOMP="BHZ BHN BHE"

#####
#	define functions - the purpose of these functions is to
#	make the processing script much easier to read and to understand
#####

#####
#	determine if the files exist and if they have 1728000 data points
#
#	fileok  fileok ${YEAR}/${DOY} ${STA1} ${STA2}
#####
fileok () {
	DOPROCESS="NO"
#	echo fileok $1 $2 $3
	DIR=$1
	STA1=$2
	STA2=$3
	for  i in ${STA1} ${STA2}
	do
		for j in ${OCOMP}
		do
			if [ ! -f ${DIR}/${i}${j} ]
			then
				return 0
			fi
			#####
			#####
			# get the file name - in case of more than one
			# use the last 
			#####
			F=`ls ${DIR}/${i}${j} | tail -1`
			if [ -z ${F} ]
			then
				return 0
			fi
			#####
			# file exists - now try to get number of points
			#####
			NPTS=`saclhdr -NPTS ${F}`
echo ${NPTS} ${F}
			DEPMAX=`saclhdr -DEPMAX ${F}`
			DEPMIN=`saclhdr -DEPMIN ${F}`
			DEPMEN=`saclhdr -DEPMEN ${F}`
#			if [ ${NPTS} -ne "1728000" ]
#			then
#				echo wrong number of points ${F} ${NPTS}
#				return 0
#			fi
			#####
			#	 this is also the place to check for maximum amplitude
			#####
			MAX1=`echo $DEPMAX $DEPMIN | awk '{printf "%d", $1 - $2}' `
			if [ ${MAX1} -gt "1000000" ]
			then
				echo amplitudes too large ${F} ${MAX1}
				return 0
			fi
		done
	done
	#####
	#	all is OK
	#####
	DOPROCESS="YES"
	echo $1 $2 $3 is OK
	return 0
}

doprep () {
	echo doprep $1 $2 $3
	DIR=$1
	STA1=$2
	STA2=$3
	#####
	#	for the station pair rotate horizontals to the great circle 
	#	 doprep ${BASE}/${YEAR}/${DOY} ${STA1} ${STA2}
	#####
	SLAT1=`saclhdr -STLA ${DIR}/${STA1}BHZ `
	SLON1=`saclhdr -STLO ${DIR}/${STA1}BHZ `
	SLAT2=`saclhdr -STLA ${DIR}/${STA2}BHZ `
	SLON2=`saclhdr -STLO ${DIR}/${STA2}BHZ `
	AZ12=`udelaz -BAZ -ELAT ${SLAT2} -ELON ${SLON2} -SLAT ${SLAT1} -SLON ${SLON1} | awk '{ printf "%d", $1 }' | awk '{print ($1 )%360}' `
	AZ21=`udelaz -BAZ -ELAT ${SLAT1} -ELON ${SLON1} -SLAT ${SLAT2} -SLON ${SLON2} | awk '{ printf "%d", $1 }' |awk '{print ($1 + 180)%360}' `
	#####
	# now rotate the traces
	#####
	F1Z=${DIR}/${STA1}BHZ
	F1N=${DIR}/${STA1}BHN
	F1E=${DIR}/${STA1}BHE
	F2Z=${DIR}/${STA2}BHZ
	F2N=${DIR}/${STA2}BHN
	F2E=${DIR}/${STA2}BHE
echo ${F1Z} ${F1N} ${F1E}
#####
#	BUG FIX FOR CSS - PUT IN CMPINC CMPAZ
#	This should not be donw here since it repeats a lot
#####
gsac << EOF
r ${F1Z} ${F2Z}
ch cmpinc 0 cmpaz 0
wh
r ${F1N} ${F2N}
ch cmpinc 90 cmpaz 0
wh
r ${F1E} ${F2E}
ch cmpinc 90 cmpaz 90
wh
EOF


gsac  << EOF
r ${F1Z} ${F1N} ${F1E}
rtr
hp c ${HPCORNER} n 2 p 2
taper w 0.02
rot3 to ${AZ12}
ch EVLA ${SLAT2} EVLO ${SLON2}
cd ${TEMP}
w ${2}BHR ${2}BHT ${2}BHZ
quit
EOF

echo ${F2Z} ${F2N} ${F2E}
gsac  << EOF
r ${F2Z} ${F2N} ${F2E}
rtr
hp c ${HPCORNER} n 2 p 2
taper w 0.02
rot3 to ${AZ21}
ch EVLA ${SLAT1} EVLO ${SLON1}
cd ${TEMP}
w ${3}BHR ${3}BHT ${3}BHZ
quit
EOF

}

docorr () {
	echo docorr $1 $2 $3 $4
	#####
	#	for the station pair create the cross correlations for each day
	#		  1       2       3        4
	#	 docorr ${YEAR} ${DOY} ${STA1}  ${STA2} 
	#####
	YEAR=$1
	DOY=$2
	STA1=$3
	STA2=$4
	for C in ${COMP}
	do
		for HS in 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 
		do
		#####
		# define end time for this segment
		#####
		case ${HS} in
			00) HE=01 ;;
			01) HE=02 ;;
			02) HE=03 ;;
			03) HE=04 ;;
			04) HE=05 ;;
			05) HE=06 ;;
			06) HE=07 ;;
			07) HE=08 ;;
			08) HE=09 ;;
			09) HE=10 ;;
			10) HE=11 ;;
			11) HE=12 ;;
			12) HE=13 ;;
			13) HE=14 ;;
			14) HE=15 ;;
			15) HE=16 ;;
			16) HE=17 ;;
			17) HE=18 ;;
			18) HE=19 ;;
			19) HE=20 ;;
			20) HE=21 ;;
			21) HE=22 ;;
			22) HE=23 ;;
		esac
		F1=${STA1}${C}
		F2=${STA2}${C}
		echo crosscorrelation $1 $2 $F1 $F2 hours ${HS}:30 - ${HE}:30
gsac > /dev/null << EOF
cd ${TEMP}
cut GMT ${YEAR} ${DOY} $HS 31 00 000 GMT ${YEAR} ${DOY} $HE 29 00 000
r $F1 $F2
rtr
hp c ${HPCORNER} n 2 p 2
whiten freqlimits ${FREQLIMITS} absolute
rtr
correlate 2 N ${SEGMENT} 
w ${HS}.xx ${HS}.xy
quit
EOF
		#####
		# end HS HE
		#####
		done
	#####
	# now do the stack for the day
	#####
gsac > /dev/null << EOF
cd ${TEMP}
r ??.xy 
stack
w ${YEAR}.${DOY}.${STA1}${C}${STA2}${C}.cor
reverse
w ${YEAR}.${DOY}.${STA1}${C}${STA2}${C}.rev
q
EOF
	#####
	# clean up
	#####
	echo clean up
	rm -f ${TEMP}/??.xy
	rm -f ${TEMP}/??.xx
	#####
	# end C
	#####
 done
	rm -f ${TEMP}/${STA1}BH?
	rm -f ${TEMP}/${STA2}BH?
}

dostack () {
	#####
	# do the complete stack for all days
	#            1       2       3      4
	# dostack ${STA1} ${CMP} ${STA2} ${CMP}
	#####
	STA1=$1
	CMP1=$2
	STA2=$3
	CMP2=$4
	##########
	#	NOW THAT THE INDIVIDUAL CROSSCORRELATIONS BY DAY ARE
	#	COMPLETED, PERFORM THE STACK. CREATE TWO STACKED FILES
	#	ONE FOR THE ENTIRE 1 hour WINDOW and another for a +- 300 sec WINDOW
	#####

#####
#	FINAL STACK
#	NOTE ONLY ZZ RR TT
#####
cd ${MYPWD}
gsac > /dev/null << EOF
r ${TEMP}/*.*.${STA1}${CMP1}${STA2}${CMP2}.rev ${TEMP}/*.*.${STA1}${CMP1}${STA2}${CMP2}.cor
lh dist
rtr
stack
w cstack.stk
cut o ${WINL} o ${WINH}
r cstack.stk
w ${STA1}${CMP1}${STA2}${CMP2}.WSTK
quit
EOF
rm -f ${TEMP}/*.*.${STA1}${CMP1}${STA2}${CMP2}.rev 
rm -f ${TEMP}/*.*.${STA1}${CMP1}${STA2}${CMP2}.cor
rm -f cstack.stk
}


#####################################################################################
#		BEGIN PROCESSING LOOP
#####################################################################################
 
for STA1 in  ACLR ASCB 
do
	DOSTK="NO"
case ${STA1} in
ACLR) STR="ASCB VELZ" ;;
ASCB) STR="VELZ" ;;
esac
	for STA2 in $STR
	do
	
	#####
	#	Now that the station pair is defined, get the data
	#	a) first determine if the SAC file exists
	#	b) determine if the SAC file has 1728000 points
	
	#####
	
		for YEAR in 2005 
		do
			for DOY in     001 002 003 004 005 006 007 008 009 \
	 			010 011 012 013 014 015 016 017 018 019 \
	 			020 021 022 023 024 025 026 027 028 029 \
	 			030 031 032 033 034 035 036 037 038 039 \
	 			040 041 042 043 044 045 046 047 048 049 \
	 			050 051 052 053 054 055 056 057 058 059 \
	 			060 061 062 063 064 065 066 067 068 069 \
	 			070 071 072 073 074 075 076 077 078 079 \
	 			080 081 082 083 084 085 086 087 088 089 \
	 			090 091 092 093 094 095 096 097 098 099 \
	 			100 101 102 103 104 105 106 107 108 109 \
	 			110 111 112 113 114 115 116 117 118 119 \
	 			120 121 122 123 124 125 126 127 128 129 \
	 			130 131 132 133 134 135 136 137 138 139 \
	 			140 141 142 143 144 145 146 147 148 149 \
	 			150 151 152 153 154 155 156 157 158 159 \
	 			160 161 162 163 164 165 166 167 168 169 \
	 			170 171 172 173 174 175 176 177 178 179 \
	 			180 181 182 183 184 185 186 187 188 189 \
	 			190 191 192 193 194 195 196 197 198 199 \
	 			200 201 202 203 204 205 206 207 208 209 \
	 			210 211 212 213 214 215 216 217 218 219 \
	 			220 221 222 223 224 225 226 227 228 229 \
	 			230 231 232 233 234 235 236 237 238 239 \
	 			240 241 242 243 244 245 246 247 248 249 \
	 			250 251 252 253 255 255 256 257 258 259 \
	 			260 261 262 263 265 265 266 267 268 269 \
	 			270 271 272 273 275 275 276 277 278 279 \
	 			280 281 282 283 285 285 286 287 288 289 \
	 			290 291 292 293 295 295 296 297 298 299 \
	 			300 301 302 303 305 305 306 307 308 309 \
	 			310 311 312 313 315 315 316 317 318 319 \
	 			320 321 322 323 325 325 326 327 328 329 \
	 			330 331 332 333 335 335 336 337 338 339 \
	 			340 341 342 343 345 345 346 347 348 349 \
	 			350 351 352 353 355 355 356 357 358 359 \
	 			360 361 362 363 365 366
			do

#				echo Processing ${STA1} ${STA2} DOY ${DOY} YEAR ${YEAR}
				#####
				#	before we process the directory must exist 
				#	then the files must exist
				#####
				if [ -d ${BASE}/${YEAR}/${DOY} ]
				then
	
					
					fileok ${BASE}/${YEAR}/${DOY} ${STA1} ${STA2}
					if [ $DOPROCESS = "YES" ]
					then
					echo  ${BASE}/${YEAR}/${DOY}    exist
						doprep  ${BASE}/${YEAR}/${DOY} ${STA1} ${STA2}
						docorr  ${YEAR} ${DOY} ${STA1}  ${STA2} 
						DOSTK="YES"
					else
						echo reject ${BASE}/${YEAR}/${DOY} ${STA1} ${STA2}
					fi
				else
					echo  ${BASE}/${YEAR}/${DOY} no exist
				fi
			#####
			#	end of DOY loop
			#####
			done

		#####
		#	end of YEAR loop
		#####
		done

		#####
		#	do final stacks for this station pair for the year
		#####
			if [ $DOSTK = "YES" ]
			then
				for CMP in ${COMP}
				do
					echo ${STA1} ${CMP} ${STA2} ${CMP}
					dostack ${STA1} ${CMP} ${STA2} ${CMP}
				done
			fi
	
	#####
	#	end of STA2 loop
	#####
	done
#####
#	end of STA1 loop
#####
done

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