EarthModelBinned.cxx

Go to the documentation of this file.

/*************************************************************************
 * EarthModelBinned.cxx
 * Implements a 3D earth model (binned in depth, latitude, and longitude)
 * for the programming library OscProb
 * Created by Rebekah Pestes (last modified: 4/4/23)
 *************************************************************************/
 
#include <cmath>
#include <fstream>
#include <iostream>
 
#include "EarthModelBinned.h"
#include "prem_default.hpp"
 
using namespace std;
 
using namespace OscProb;
 
//.............................................................................
void TrajConstants::UpdateNuAngles(double cosTheta, double phi)
{
  cosT   = cosTheta;
  sinSqT = 1 - cosT * cosT;
  sinT   = sqrt(sinSqT);
 
  cosA = cos(phi);
  sinA = sin(phi);
 
  sinTsinA = sinT * sinA;
  sinTcosA = sinT * cosA;
}
 
//.............................................................................
void TrajConstants::UpdateDetPos(double rDet, double DetLat, double DetLon)
{
  sinDetLat = sin(DetLat);
  cosDetLat = cos(DetLat);
  cosDetLon = cos(DetLon);
  sinDetLon = sin(DetLon);
  DetRadius = rDet;
 
  rDetSinDetLat = DetRadius * sinDetLat;
  rDetCosDetLat = DetRadius * cosDetLat;
}
 
//.............................................................................
void TrajConstants::Recalculate()
{
  cosTcosDetLat      = cosT * cosDetLat;
  sinTsinAsinDetLon  = sinTsinA * sinDetLon;
  sinTsinAcosDetLon  = sinTsinA * cosDetLon;
  rDetCosT           = DetRadius * cosT;
  rDetSinT           = DetRadius * sinT;
  alpha              = sinTcosA * sinDetLat - cosTcosDetLat;
  beta               = sinT * sinDetLat - cosA * cosTcosDetLat;
  gamma              = sinTcosA * cosDetLat + cosT * sinDetLat;
  gammaSq            = gamma * gamma;
  rDetGammaSinDetLat = rDetSinDetLat * gamma;
  maxSinSqLat        = 1 - pow(sinA * cosDetLat, 2);
  if (beta == 0) { xlatextreme = 0; }
  else {
    xlatextreme = cosA * rDetCosDetLat / beta;
  }
}
 
//.............................................................................
EarthModelBinned::EarthModelBinned(string filename)
{
  SetDetPos(6368, 0, 0);
  LoadModel(filename);
  SetRemoveSmallPaths(false);
  fLonError = 0;
}
 
//.............................................................................
EarthModelBinned::~EarthModelBinned() {}
 
//.............................................................................
vector<EarthBin> EarthModelBinned::GetEarthBins() { return fEarthBins; }
 
//.............................................................................
void EarthModelBinned::ClearModel()
{
  fnDepthBins = 0;
  fnLonBins   = 0;
  fnLatBins   = 0;
  fEarthBins.clear();
}
 
//.............................................................................
void EarthModelBinned::SetDetPos(double rad, double lat, double lon)
{
  SetDetectorCoordinates(rad, lat, lon);
  fC.UpdateDetPos(fDetRadius, fDetLat, fDetLon);
}
 
//.............................................................................
void EarthModelBinned::AddBin(double radius_out, double radius_in,
                              double latitude, double longitude, double density,
                              double zoa, double index)
{
  fEarthBins.push_back(EarthBin(radius_out, radius_in, latitude, longitude,
                                density, zoa, index));
}
 
//.............................................................................
void EarthModelBinned::LoadModel(string filename)
{
  // Clear the current model
  ClearModel();
  double EarthRadius = 6371.0;
 
  // Use default if no file provided
  if (filename == "") { filename = PREM3D_DEFAULT; }
  cout << "Loading Earth model table from " << filename << "..." << endl;
 
  // Open the file
  ifstream fin;
  fin.open(filename.c_str());
  if (!fin) {
    cerr << "ERROR: File " << filename << " not found!" << endl;
    return;
  }
 
  // Variables for storing table rows
  float depth, latitude, longitude, density, zoa, index;
 
  // Keep track of previous depth/latitude/longitude
  double radius_prev = 0;
  double lat_prev    = -90.0;
  double lon_prev    = 0;
 
  // Hold onto previous different depth, so we have depth maximum for bin
  double radius_min = radius_prev;
 
  // Initialize max region index
  fmaxRegIndex = -1;
 
  // Loop over table rows
  while (fin >> longitude >> latitude >> depth >> density >> zoa >> index) {
    double outer_radius = EarthRadius - depth;
 
    // Move minimum depth (and reset latitude/longitude) if depth has changed
    // from previous one
    if (outer_radius > radius_prev) {
      radius_min = radius_prev;
      lat_prev   = -90.0;
      lon_prev   = 0.0;
      fnDepthBins++;
    }
    else if (outer_radius <
             radius_prev) { // Depths must be ordered in model file
      cerr << "ERROR: Depths are not sorted in decreasing order in the model "
              "file "
           << "(or depth is greater than " << EarthRadius << " km)." << endl;
      ClearModel();
      return;
    }
    else if (longitude > lon_prev) { // Reset latitude if longitude has changed
                                     // from previous one
      lat_prev = -90.0;
      fnLonBins++;
    }
    else if (longitude < lon_prev) { // Longitudes must be ordered in model file
                                     // (after depths)
      cerr << "ERROR: Longitudes are not sorted in increasing order (after "
              "depths) in the model file."
           << endl;
      ClearModel();
      return;
    }
    else if (latitude < lat_prev) { // Latitudes must be ordered in model file
                                    // (after longitudes)
      cerr << "ERROR: Latitudes are not sorted in increasing order (after "
              "longitudes) in the model file."
           << endl;
      ClearModel();
      return;
    }
 
    // Add this bin to the model
    AddBin(outer_radius, radius_min, latitude, longitude, density, zoa, index);
 
    // Check region index
    if (index > fmaxRegIndex) { fmaxRegIndex = index; }
 
    // Set previous coordinates for next bin
    radius_prev = outer_radius;
    lat_prev    = latitude;
    lon_prev    = longitude;
  }
 
  // Set Outermost Radius
  fRadiusMax = radius_prev;
 
  // Modify number of longitude bins to how many there actually are (instead of
  // how many times we changed the longitude bin without changing the depth bin
  // = nLonBins-1 times per depth bin)
  fnLonBins        = fnLonBins / fnDepthBins + 1;
  fInvLonBinWidth  = 0.5 * fnLonBins / M_PI; // in radians^(-1)
  fHalfLonBinWidth = M_PI / fnLonBins;       // in radians
 
  // Calculate number of latitude bins
  fnLatBins        = fEarthBins.size() / (fnLonBins * fnDepthBins);
  fInvLatBinWidth  = fnLatBins / M_PI;       // in radians^(-1)
  fHalfLatBinWidth = 0.5 * M_PI / fnLatBins; // in radians
 
  // Error if first latitude bin center deviates from about half of the
  // calculated latitude bin width
  if (abs(fHalfLatBinWidth - (fEarthBins[0].latitude + 0.5 * M_PI)) >
      fHalfLatBinWidth * 0.01) { // Warn if first latitude bin center isn't
                                 // about half of the latitude bin width
    cerr << "ERROR: Latitude of 1st bin center (" << fEarthBins[0].latitude
         << ") is not in the middle of the first bin, whose size should be "
         << 2 * fHalfLatBinWidth << " rad, as calculated from having "
         << fnLatBins << " latitude bins." << endl;
    ClearModel();
    return;
  }
 
  cout << "\t...done (" << fEarthBins.size() << " bins: " << fnDepthBins
       << " depth, " << fnLatBins << " latitude, and " << fnLonBins
       << " longitude)." << endl;
}
 
//.............................................................................
void EarthModelBinned::SetRegionZoA(int index, double zoa)
{
  if (index > fmaxRegIndex) {
    cerr << "WARNING: Index " << index << " does not exist in model "
         << "(max index = " << fmaxRegIndex << ")." << endl
         << "Doing nothing." << endl;
    return;
  }
 
  int nbins = fEarthBins.size();
 
  // Loop over all bins and change the ones
  // with the given index
  for (int i = 0; i < nbins; i++) {
    if (fEarthBins[i].index != index) continue;
 
    fEarthBins[i].zoa = zoa;
  }
}
 
//.............................................................................
double EarthModelBinned::GetRegionZoA(int index)
{
  if (index > fmaxRegIndex) {
    cerr << "ERROR: Index " << index << " does not exist in model "
         << "(max index = " << fmaxRegIndex << "). Returning 0" << endl;
    return 0;
  }
 
  int nbins = fEarthBins.size();
 
  for (int i = 0; i < nbins; i++) {
    if (fEarthBins[i].index != index)
      continue;
 
    else
      return fEarthBins[i].zoa;
  }
 
  // End of vector reached without finding index
  cerr << "ERROR: Index " << index << " not found! Returning 0" << endl;
  return 0;
}
 
//.............................................................................
void EarthModelBinned::ScaleRegionDensity(int index, double factor)
{
  if (index > fmaxRegIndex) {
    cerr << "WARNING: Index " << index << " does not exist in model "
         << "(max index = " << fmaxRegIndex << ")." << endl
         << "Doing nothing." << endl;
    return;
  }
 
  int nbins = fEarthBins.size();
 
  // Loop over all bins and change the ones
  // with the given index
  for (int i = 0; i < nbins; i++) {
    if (fEarthBins[i].index != index) continue;
 
    double density_init   = fEarthBins[i].density;
    fEarthBins[i].density = factor * density_init;
  }
}
 
//.............................................................................
void EarthModelBinned::AddPath(double length, EarthBin bin)
{
  AddPathSegment(length, bin.density, bin.zoa, bin.index);
}
 
//.............................................................................
double EarthModelBinned::DetDistForNextLatBin(int cur_index, LatBinInfo& L)
{
  // Find next lat bin boundary (bin center +/- 1/2-binwidth)
  double NewLat = fEarthBins[cur_index].latitude + L.sign * L.dLat;
  // Turn around if lat reaches +/- 90 degrees (or close enough)
  if (abs(NewLat) - 0.5 * M_PI > -1.0e-14) {
    if (L.maxreached) return -1; // Only turn around once
    L.maxreached = true;
    //    cout << "Max latitude reached..." << endl;
    L.sign = -L.sign;
    NewLat = fEarthBins[cur_index].latitude + L.sign * L.dLat;
    if (abs(NewLat) - 0.5 * M_PI > -1.0e-14)
      return -1; // This would happen if there's only one latitude bin...
  }
 
  // Calculate parts of answer
  double sinNewLat    = sin(NewLat);
  double NsinSqNewLat = -sinNewLat * sinNewLat;
  double radicand     = fC.maxSinSqLat +
                    NsinSqNewLat; // 1-[sin(NewLat)]^2-[sin(phi)cos(DetLat)]^2
  // Turn around if the radicand is 0 (close enough) or less
  if (radicand < 1.0e-14) {
    if (L.maxreached) return -1; // Only turn around once
    L.maxreached = true;
    //    cout << "Max latitude reached..." << endl;
    L.sign       = -L.sign;
    NewLat       = fEarthBins[cur_index].latitude + L.sign * L.dLat;
    sinNewLat    = sin(NewLat);
    NsinSqNewLat = -sinNewLat * sinNewLat;
    radicand     = fC.maxSinSqLat + NsinSqNewLat;
    if (radicand < 1.0e-14)
      return -1; // The rest of the track is contained in this lat bin
  }
  double denominator = fC.gammaSq + NsinSqNewLat;
 
  // Calculate distance from the detector at next lat bin boundary
  double answer;
  if (denominator == 0) {
    answer = 0.5 * (fC.xlatextreme - fC.rDetSinDetLat / fC.gamma);
    // There is a "possible" case where answer would be 0:
    //   if the detector is at the center of the Earth... just don't...
  }
  else {
    answer = -(fC.rDetCosT * NsinSqNewLat + fC.rDetGammaSinDetLat +
               L.sign * fC.rDetSinT * sinNewLat * sqrt(radicand)) /
             denominator;
  }
 
  // Check to see if going to the edge of the bin jumped to an invalid
  // part of the function (cosA*rDetCosDetLat/beta = x at max lat)
  if (L.maxreached && answer >= fC.xlatextreme) {
    return -1; // The rest of the track is contained in this lat bin
  }
 
  return answer;
}
 
//.............................................................................
double EarthModelBinned::DetDistForNextLonBin(double prev_lon, LonBinInfo& L)
{
  double twoPi = 2 * M_PI;
  // Increment to next bin in direction dLon
  double lon = prev_lon + L.dLon;
  if (lon < 0 || lon >= twoPi) { lon -= floor(lon / twoPi) * twoPi; }
  if (L.dLon > 0) {
    L.nextBin = L.bin + 1;
    if (L.nextBin >= fnLonBins) {
      // cross from 2pi to 0
      L.nextBin = 0;
    }
  }
  else {
    L.nextBin = L.bin - 1;
    if (L.nextBin < 0) {
      // cross from 0 to 2pi
      L.nextBin = fnLonBins - 1;
    }
  }
 
  // Check if lon is outside of range covered by neutrino's trajectory
  if (L.max < L.min) {
    if (0 > lon - L.min && L.max - lon < 0) return -1;
  }
  else {
    if (0 > lon - L.min || L.max - lon < 0) return -1;
  }
 
  // Parts of the answer
  double cosNewLon   = cos(lon);
  double sinNewLon   = sin(lon);
  double expression  = fC.sinDetLon * cosNewLon - sinNewLon * fC.cosDetLon;
  double denominator = fC.alpha * expression +
                       fC.sinTsinAcosDetLon * cosNewLon +
                       fC.sinTsinAsinDetLon * sinNewLon;
 
  // Calculate & return detector distance
  if (denominator == 0) { // Is there ever a case when this is true???
    // Error message will be printed at the end of the FillPath function
    L.err_message.assign("Longitude along neutrino trajectory = " +
                         std::to_string(lon));
    L.error = -1;
    return -1;
  }
  return fC.rDetCosDetLat * expression / denominator;
}
 
//.............................................................................
void EarthModelBinned::RecordLatLonBinCrossings(double  detDist_nextDbin,
                                                double& DetDist, int& index,
                                                LatBinInfo& latI,
                                                LonBinInfo& lonI)
{
  // Keep going as long as there is a lat or lon bin crossing before
  // the next depth bin crossing
  while (detDist_nextDbin < latI.detDist_nextBin ||
         detDist_nextDbin < lonI.detDist_nextBin) {
    if (latI.detDist_nextBin >
        lonI.detDist_nextBin) { // if next lat bin crossing is before next lon
                                // bin crossing
      // Add segment to next lat bin to path
      AddPath(DetDist - latI.detDist_nextBin, fEarthBins[index]);
      double latitude_crossing =
          fEarthBins[index].latitude + latI.sign * latI.dLat;
 
      // Move lat bins
      DetDist = latI.detDist_nextBin;
      index += latI.nextBin - latI.bin;
      latI.bin = latI.nextBin;
      //      cout << "Latitude bin crossing at x=" << DetDist << "km,lat=" <<
      //      latitude_crossing << "rad...now in Earth bin " << index << " = (r"
      //      << floor(index/(fnLonBins*fnLatBins)) << ", lon" << lonI.bin << ",
      //      lat" << latI.bin << ")" << endl;
 
      // Find next lat bin
      if (latI.dLat ==
          0) { // Only 1 lat bin crossing along entire neutrino trajectory
        latI.detDist_nextBin = -1;
      }
      else {
        latI.detDist_nextBin = DetDistForNextLatBin(index, latI);
        latI.nextBin         = latI.bin + latI.sign;
      }
    }
    else { // if next lon bin crossing is before (or simultaneous with) next lat
           // bin corssing
      // Add segment to next lon bin to path
      AddPath(DetDist - lonI.detDist_nextBin, fEarthBins[index]);
 
      // Move lon bins
      DetDist = lonI.detDist_nextBin;
      index += (lonI.nextBin - lonI.bin) * fnLatBins;
      lonI.bin = lonI.nextBin;
      //      cout << "Longitude bin crossing at x=" << DetDist << "km...now in
      //      Earth bin " << index << " = (r" <<
      //      floor(index/(fnLonBins*fnLatBins)) << ", lon" << lonI.bin << ",
      //      lat" << latI.bin << ")" << endl;
 
      // Find next lon bin
      if (lonI.dLon ==
          0) { // Only 1 lon bin crossing along entire neutrino trajectory
        lonI.detDist_nextBin = -1;
      }
      else {
        lonI.detDist_nextBin =
            DetDistForNextLonBin(fEarthBins[index].longitude, lonI);
      }
    }
  }
}
 
//.............................................................................
int EarthModelBinned::LonBinIndex(double longitude)
{
  // Find difference between longitude and beginning of first lon bin
  double londiff = longitude - (fEarthBins[0].longitude - fHalfLonBinWidth);
 
  // Shift londiff to be between 0 and 2pi
  if (londiff < 0 || londiff >= 2 * M_PI)
    londiff -= floor(londiff / (2 * M_PI));
 
  // Number of binwidths after from first lon bin
  return floor(londiff * fInvLonBinWidth);
}
 
//.............................................................................
int EarthModelBinned::FillPath(double cosT, double phi)
{
  // Clear current path sequence
  fNuPath.clear();
 
  // Do nothing if cosine is unphysical
  if (fabs(cosT) > 1) return 0;
 
  // Calculate useful combinations of variables
  double Az = phi / 180.0 * M_PI; // convert from degrees to radians
  fC.UpdateNuAngles(cosT, Az);
  fC.Recalculate();
  double rDetSqSinSqT = pow(fC.rDetSinT, 2);
  double distEdgeToMin =
      sqrt(fRadiusMax * fRadiusMax -
           rDetSqSinSqT); // distance from point where neutrino enters to
                          // minimum depth (neglecting stopping at the detector)
  double baseline = distEdgeToMin - fC.rDetCosT; // total baseline
 
  // Define the maximum depth (minimum radius) along the path (detector
  // depth/radius, if cosT is non-negative)
  double minRadius = fDetRadius;
  if (cosT < 0) { minRadius = fC.rDetSinT; }
 
  // Starting latitude
  /*  cout << "Pieces of init_lat calculation:" << endl
         << "\trDetSinT = " << fC.rDetSinT << endl
         << "\tbeta = " << fC.beta << endl
         << "\tgamma = " << fC.gamma << endl
         << "\tdistEdgeToMin = " << distEdgeToMin << endl
         << "\tRadiusMax = " << fRadiusMax << endl; */
  double init_lat = asin((fC.rDetSinT * fC.beta + fC.gamma * distEdgeToMin) /
                         fRadiusMax); // in radians
 
  // Starting longitude
  double init_lon = atan2(
      fC.sinDetLon * fC.rDetCosDetLat -
          baseline * (fC.alpha * fC.sinDetLon + fC.sinTsinAcosDetLon),
      fC.cosDetLon * fC.rDetCosDetLat +
          baseline *
              (fC.sinTsinAsinDetLon -
               fC.alpha * fC.cosDetLon)); // in radians (between -M_PI and M_PI)
  if (init_lon < 0) { init_lon += 2 * M_PI; }
 
  // Get initial Earth bin and # of bins per depth layer
  LatBinInfo latinfo;
  LonBinInfo loninfo;
  int        binsPerDepth = fEarthBins.size() / fnDepthBins;
  int        init_lonBin  = LonBinIndex(init_lon);
  int        init_latBin  = floor((init_lat + M_PI / 2.0) * fInvLatBinWidth);
  int        index =
      fEarthBins.size() - binsPerDepth + init_lonBin * fnLatBins + init_latBin;
 
  latinfo.bin = init_latBin;
  loninfo.bin = init_lonBin;
  //  cout << "Neutrino starting at (lat,lon)=(" << init_lat << "rad," <<
  //  init_lon << "rad) in Earth bin " << index << " = (r" <<
  //  floor(index/binsPerDepth) << ", lon" << loninfo.bin << ", lat" <<
  //  latinfo.bin << ")" << endl;
 
  /* Find detector distances for 1st latitude/longitude bin changes */
  latinfo.dLat = fHalfLatBinWidth; // change in lat from bin center to next bin
                                   // (0 if no more than 1 lat bin change)
  // Condition calculation for starting sign in x(lat) equation
  latinfo.sign = 1; //+ => lat incr., - => lat decr. (with decr. x)
  latinfo.maxreached =
      false; // this changes to true if past the lat func transition
  if (fC.beta < 0) latinfo.sign = -1;
  if (fC.beta == 0) {
    if (fC.cosA > 0) latinfo.sign = -1;
  }
  else if (baseline < fC.xlatextreme) {
    latinfo.sign       = -(latinfo.sign);
    latinfo.maxreached = true;
  }
  // Condition calculation for incr./decr. Longitude (with decr. x)
  loninfo.dLon =
      fHalfLonBinWidth;   // change in lon from bin center to next bin (includes
                          // direction; 0 if no more than 1 lon bin change)
  loninfo.min = init_lon; // initial lon if lon is inc.
  loninfo.max = fDetLon;  // detector lon if lon is inc.
  if (fC.sinA < 0) {
    loninfo.dLon = -(loninfo.dLon);
    loninfo.min  = loninfo.max;
    loninfo.max  = init_lon;
    // Switch max/min if dec., instead of inc.
    loninfo.max = init_lon;
    loninfo.min = fDetLon;
  }
  // Remainder of detDist (and nextBin) calculations
  if (fC.sinT == 0) {
    // both latitude and logitude flip at center of the Earth
    loninfo.detDist_nextBin =
        -fC.rDetCosT; // formula is divided by cosT, but cosT = +-1
    loninfo.nextBin         = LonBinIndex(fDetLon);
    loninfo.dLon            = 0;
    latinfo.detDist_nextBin = loninfo.detDist_nextBin;
    latinfo.nextBin         = floor((fDetLat + 0.5 * M_PI) * fInvLatBinWidth);
    latinfo.dLat            = 0;
  }
  else {
    if (fC.sinA == 0) {
      if (fC.alpha == 0) { // crosses "at" infinity
        loninfo.detDist_nextBin = -1;
      }
      else { // longitude flips at center of the Earth
        loninfo.detDist_nextBin = fC.rDetCosDetLat / fC.alpha;
        loninfo.nextBin         = LonBinIndex(fDetLon);
      }
      loninfo.dLon = 0;
    }
    else if (fC.cosDetLat == 0) {
      loninfo.detDist_nextBin = -1;
    }
    else {
      // x(lon) calculation
      loninfo.detDist_nextBin = DetDistForNextLonBin(
          fEarthBins[index].longitude, loninfo); // also sets loninfo.nextBin
    }
    // x(lat) calculation
    latinfo.detDist_nextBin = DetDistForNextLatBin(index, latinfo);
    latinfo.nextBin         = init_latBin + latinfo.sign;
  }
 
  // Get initial Detector Distance
  double prev_DetDistance = baseline;
 
  /* Start at the top layer of Earth model and go down to minimum (not including
   * minimum), looping over depth bins */
  int dBin = 0; // depth bin index
  for (dBin = 0; dBin < fnDepthBins - 1; dBin++) {
    // Check if minimum radius is contained in depth bin
    double r_in = fEarthBins[index].radius_in; // inner-most radius
    if (r_in < minRadius) break;
 
    // Distance from detector at inner-most radius
    double DetDistance = -fC.rDetCosT + sqrt(r_in * r_in - rDetSqSinSqT);
 
    // Check for latitude/longitude bin crossings
    RecordLatLonBinCrossings(DetDistance, prev_DetDistance, index, latinfo,
                             loninfo);
 
    // Add Segment to next depth bin to path
    AddPath(prev_DetDistance - DetDistance, fEarthBins[index]);
 
    // Reset variables for next depth bin
    prev_DetDistance = DetDistance;
    index -= binsPerDepth;
    //    cout << "Depth bin crossing at x=" << DetDistance << "km...now in
    //    Earth bin " << index << " = (r" << floor(index/binsPerDepth) << ",
    //    lon" << loninfo.bin << ", lat" << latinfo.bin << ")" << endl;
  }
 
  /* Do minimum depth bin and then go up to the detector */
  int index2 = index;
  for (index2 = index; fEarthBins[index2].radius_out < fDetRadius;
       index2 += binsPerDepth) {
    double r_out = fEarthBins[index2].radius_out; // outer-most radius
 
    // Distance from detector at outer-most radius
    double DetDistance = -fC.rDetCosT - sqrt(r_out * r_out - rDetSqSinSqT);
 
    // Check for latitude/longitude bin crossings
    RecordLatLonBinCrossings(DetDistance, prev_DetDistance, index2, latinfo,
                             loninfo);
 
    // Add Segment in depth bin to path
    AddPath(prev_DetDistance - DetDistance, fEarthBins[index2]);
 
    // Reset variables for "previous" depth bin
    prev_DetDistance = DetDistance;
    //    cout << "Depth bin crossing at x=" << DetDistance << "km...now in
    //    Earth bin " << index2 << " = (r" <<
    //    floor((index2+binsPerDepth)/binsPerDepth) << ", lon" << loninfo.bin <<
    //    ", lat" << latinfo.bin << ")" << endl;
  }
 
  /* Do path to detector */
  // Check for latitude/longitude bin crossings
  RecordLatLonBinCrossings(0, prev_DetDistance, index2, latinfo, loninfo);
  // Add the path segment within the detector bin
  AddPath(prev_DetDistance, fEarthBins[index2]);
 
  // Longitude calculation error message, if needed
  if (loninfo.error < 0) {
    cerr
        << "ERROR:  Oops... It turns out that I was wrong about the denominator"
        << " of the x(long) equation never being 0 apart from the cases where "
        << "sin(Zenith) = 0, sin(Azimuthal) = 0, or cos(Detector Longitude) = "
        << "0."
        << "  Send the following message to rpestes@apc.in2p3.fr:" << endl
        << endl
        << "You were wrong about the denominator of x(long)... "
        << "Here are the values of the variables that broke it:" << endl
        << "\tDetector Coordinates (r,lat,lon) = (" << fDetRadius << ", "
        << fDetLat << ", " << fDetLon << ")" << endl
        << "\tcos(Zenith Angle) = " << cosT << endl
        << "\tAzimuthal Angle = " << phi << " deg" << endl
        << "\t" << loninfo.err_message << endl
        << "Please fix this ASAP!" << endl;
    return -1;
  }
 
  // Return the number of path segments
  return fNuPath.size();
}