Azimuthal Gap Downweighting option for local grids

glasscore/glasslib/include/Hypo.h

@@ -81,7 +81,8 @@ class CHypo {
 			std::string web, double bayes, double thresh, int cut,
 			std::shared_ptr<traveltime::CTravelTime> firstTrav,
 			std::shared_ptr<traveltime::CTravelTime> secondTrav,
-			std::shared_ptr<traveltime::CTTT> ttt, double resolution = 100);
+			std::shared_ptr<traveltime::CTTT> ttt, double resolution = 100,
+			double aziTaper = 360.);
 
 	/**
 	 * \brief CHypo alternate constructor
@@ -148,7 +149,7 @@ class CHypo {
 					std::shared_ptr<traveltime::CTravelTime> firstTrav,
 					std::shared_ptr<traveltime::CTravelTime> secondTrav,
 					std::shared_ptr<traveltime::CTTT> ttt, double resolution =
-							100);
+					100, double aziTap = 360.);
 
 	/**
 	 * \brief Add pick reference to this hypo
@@ -478,6 +479,15 @@ class CHypo {
 									double tStart, double tStop, int nucleate =
 											0);
 
+	/**
+	 * Calculates azimuthal gap for a proposed location
+	 *
+	 * \param lat - latitude of test location
+	 * \param lon - longitude of test location
+	 * \param z - depth of test location
+	 */
+	double gap(double lat, double lon, double z);
+
 	/**
 	 * Gets the stack of associated arrivals at location
 	 *
@@ -549,6 +559,12 @@ class CHypo {
 	 */
 	void trap();
 
+	/**
+	 * \brief aziTaper
+	 * \return the aziTaper
+	 */
+	double getAziTaper() const;
+
 	/**
 	 * \brief Latitude getter
 	 * \return the latitude
@@ -853,6 +869,11 @@ class CHypo {
 	 */
 	double searchVals[5];
 
+	/**
+	 * \brief where to taper bayesVal from azi gap
+	 */
+	double aziTaper;
+
 	/**
 	 * \brief An integer value containing this hypo's processing cycle count
 	 */

glasscore/glasslib/include/Web.h

@@ -92,7 +92,7 @@ class CWeb {
 			std::shared_ptr<traveltime::CTravelTime> firstTrav,
 			std::shared_ptr<traveltime::CTravelTime> secondTrav,
 			int numThreads = 0, int sleepTime = 100,
-			int checkInterval = 60);
+			int checkInterval = 60, double aziTaper = 360.);
 
 	/**
 	 * \brief CWeb destructor
@@ -157,7 +157,8 @@ class CWeb {
 					int numNucleate, int resolution, int numRows, int numCols,
 					int numZ, bool update,
 					std::shared_ptr<traveltime::CTravelTime> firstTrav,
-					std::shared_ptr<traveltime::CTravelTime> secondTrav);
+					std::shared_ptr<traveltime::CTravelTime> secondTrav,
+					double aziTap);
 
 	/**
 	 * \brief Generate a local detection grid
@@ -338,6 +339,13 @@ class CWeb {
 	 */
 	void addJob(std::function<void()> newjob);
 
+	/**
+	 * \brief aziTapre Getter
+	 * \return double with the azi taper start
+	 */
+	double getAziTaper() const;
+
+
 	/**
 	 * \brief CGlass getter
 	 * \return the CGlass pointer
@@ -495,6 +503,12 @@ class CWeb {
 	 */
 	double dResolution;
 
+	/**
+	 * \brief A double which describes where the locator should start
+	 * down weighting for azimuthal gap
+	**/
+	double aziTaper = 360.;
+
 	/**
 	 * \brief A boolean flag that stores whether to update this web when a
 	 * station has changed.

glasscore/glasslib/src/Hypo.cpp

@@ -52,13 +52,15 @@ CHypo::CHypo(double lat, double lon, double z, double time, std::string pid,
 				std::string web, double bayes, double thresh, int cut,
 				std::shared_ptr<traveltime::CTravelTime> firstTrav,
 				std::shared_ptr<traveltime::CTravelTime> secondTrav,
-				std::shared_ptr<traveltime::CTTT> ttt, double resolution) {
+				std::shared_ptr<traveltime::CTTT> ttt, double resolution,
+				double aziTaper) {
+
 	// seed the random number generator
 	std::random_device randomDevice;
 	m_RandomGenerator.seed(randomDevice());
 
 	if (!initialize(lat, lon, z, time, pid, web, bayes, thresh, cut, firstTrav,
-					secondTrav, ttt, resolution)) {
+					secondTrav, ttt, resolution, aziTaper)) {
 		clear();
 	}
 }
@@ -93,7 +95,8 @@ CHypo::CHypo(std::shared_ptr<CTrigger> trigger,
 					trigger->getWeb()->getThresh(),
 					trigger->getWeb()->getNucleate(),
 					trigger->getWeb()->getTrv1(), trigger->getWeb()->getTrv2(),
-					ttt, trigger->getResolution())) {
+					ttt, trigger->getResolution(),
+					trigger->getWeb()->getAziTaper())) {
 		clear();
 	}
 }
@@ -412,11 +415,16 @@ void CHypo::annealingLocate(int nIter, double dStart, double dStop,
 		return;
 	}
 
+	// taper to lower val if large azimuthal gap
+	glassutil::CTaper taperGap;
+	taperGap = glassutil::CTaper(0., 0., aziTaper, 360.);
+
 	// these hold the values of the initial, current, and best stack location
 	double valStart = 0;
 	double valBest = 0;
 	// calculate the value of the stack at the current location
-	valStart = getBayes(dLat, dLon, dZ, tOrg, nucleate);
+	valStart = getBayes(dLat, dLon, dZ, tOrg, nucleate)
+		* taperGap.Val(gap(dLat,dLon,dZ));
 
 	char sLog[1024];
 
@@ -483,7 +491,8 @@ void CHypo::annealingLocate(int nIter, double dStart, double dStop,
 		double oT = tOrg + dt;
 
 		// get the stack value for this hypocenter
-		double val = getBayes(xlat, xlon, xz, oT, nucleate);
+		double val = getBayes(xlat, xlon, xz, oT, nucleate)
+				* taperGap.Val(gap(xlat,xlon,xz));
 
 		// if testing locator print iteration
 		if (pGlass->getTestLocator()) {
@@ -497,6 +506,7 @@ void CHypo::annealingLocate(int nIter, double dStart, double dStop,
 
 		// is this stacked bayesian value (val) better than the previous best
 		// (valBest)
+
 		if (val > valBest
 				|| (val > dThresh
 						&& (valBest - val)
@@ -1179,6 +1189,55 @@ std::shared_ptr<json::Object> CHypo::expire(bool send) {
 	return (expire);
 }
 
+// ---------------------------------------------------------Gap
+double CHypo::gap(double lat, double lon, double z) {
+	// lock mutex for this scope
+	std::lock_guard<std::recursive_mutex> guard(hypoMutex);
+
+	// set up a geographic object for this hypo
+	glassutil::CGeo geo;
+	geo.setGeographic(lat, lon, EARTHRADIUSKM - z);
+
+	// create and populate vectors containing the
+	// pick distances and azimuths
+	std::vector<double> azm;
+
+	for (auto pick : vPick) {
+		// get the site
+		std::shared_ptr<CSite> site = pick->getSite();
+
+		// compute the azimuth
+		double azimuth = geo.azimuth(&site->getGeo()) / DEG2RAD;
+
+		// add to azimuth vector
+		azm.push_back(azimuth);
+	}
+
+	int nazm = azm.size();
+
+	if(nazm <= 1)
+	{
+		return 360.;
+	}
+	// sort the azimuths
+	sort(azm.begin(), azm.end());
+
+	// add the first (smallest) azimuth to the end by adding 360
+	azm.push_back(azm.front() + 360.0);
+
+	// compute gap
+	double tempGap = 0.0;
+
+	for (int i = 0; i < nazm; i++) {
+		double gap = azm[i + 1] - azm[i];
+		if (gap > tempGap) {
+			tempGap = gap;
+		}
+	}
+
+	return tempGap;
+}
+
 // ---------------------------------------------------------Gaussian
 double CHypo::gauss(double avg, double std) {
 	// generate Gaussian pseudo-random number using the
@@ -1227,6 +1286,12 @@ double CHypo::getResidual(std::shared_ptr<CPick> pick) {
 	return tRes;
 }
 
+double CHypo::getAziTaper() const {
+	std::lock_guard<std::recursive_mutex> hypoGuard(hypoMutex);
+	return (aziTaper);
+}
+
+
 double CHypo::getBayes() const {
 	std::lock_guard<std::recursive_mutex> hypoGuard(hypoMutex);
 	return (dBayes);
@@ -1903,7 +1968,7 @@ bool CHypo::initialize(double lat, double lon, double z, double time,
 						std::shared_ptr<traveltime::CTravelTime> firstTrav,
 						std::shared_ptr<traveltime::CTravelTime> secondTrav,
 						std::shared_ptr<traveltime::CTTT> ttt,
-						double resolution) {
+						double resolution, double aziTap) {
 	// lock mutex for this scope
 	std::lock_guard<std::recursive_mutex> guard(hypoMutex);
 
@@ -1917,7 +1982,7 @@ bool CHypo::initialize(double lat, double lon, double z, double time,
 	sWebName = web;
 	dBayes = bayes;
 	dBayesInitial = bayes;
-
+	aziTaper = aziTap;
 	dThresh = thresh;
 	nCut = cut;
 	dRes = resolution;

glasscore/glasslib/src/HypoList.cpp

@@ -1133,7 +1133,7 @@ bool CHypoList::mergeCloseEvents(std::shared_ptr<CHypo> hypo) {
 							glassutil::CPid::pid(), "Merged Hypo", 0.0,
 							hypo->getThresh(), hypo->getCut(), hypo->getTrv1(),
 							hypo->getTrv2(), pGlass->getTTT(),
-							distanceCut * DEG2KM);
+							distanceCut * DEG2KM, hypo->getAziTaper());
 
 					// set hypo glass pointer and such
 					hypo3->setGlass(pGlass);

glasscore/glasslib/src/Node.cpp

@@ -509,11 +509,14 @@ std::string CNode::getSitesString() {
 	for (const auto &link : vSite) {
 		// get the site
 		std::shared_ptr<CSite> currentSite = std::get< LINK_PTR>(link);
+		double lat, lon, r;
+
+		currentSite->getGeo().getGeographic(&lat, &lon, &r);
 
 		siteString += sPid + "," + currentSite->getScnl() + ","
-				+ std::to_string(currentSite->getGeo().dLat) + ";"
-				+ std::to_string(currentSite->getGeo().dLon) + ";"
-				+ std::to_string(currentSite->getGeo().dRad) + "\n";
+				+ std::to_string(lat) + ";"
+				+ std::to_string(lon) + ";"
+				+ std::to_string(r) + "\n";
 	}
 
 	return (siteString);