/************************************************************************* * 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 #include #include #include "EarthModelBinned.h" #include "prem_default.hpp" using namespace std; using namespace OscProb; //............................................................................. /// /// Update values of zenith angle and azimuthal angle for neutrino /// trajectory. /// /// @param cosTheta - Cosine of the zenith angle /// @param phi - The azimuthal angle in radians /// 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; } //............................................................................. /// /// Update detector position for neutrino trajectory calculations. /// /// @param rDet - The distance from the center of the Earth to the detector /// (in km) /// @param DetLat - The latitude of the detector (in rad) /// @param DetLon - The longitude of the detector (in rad) /// void TrajConstants::UpdateDetPos(double rDet, double DetLat, double DetLon) { sinDetLat = sin(DetLat); cosDetLon = cos(DetLon); sinDetLon = sin(DetLon); DetRadius = rDet; cosDetLat = 1 - sinDetLat * sinDetLat; rDetSinDetLat = DetRadius * sinDetLat; rDetCosDetLat = DetRadius * cosDetLat; } //............................................................................. /// /// Recalculate constants that use combinations of detector position /// variables and neutrino direction angles. /// void TrajConstants::Recalculate() { cosTcosDetLat = cosT * cosDetLat; sinTsinAsinDetLon = sinTsinA * sinDetLon; sinTsinAcosDetLon = sinTsinA * cosDetLon; rDetCosT = DetRadius * cosT; rDetSinT = DetRadius * sinT; rDetCosAcosDetLat = cosA * rDetCosDetLat; 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); } //............................................................................. /// /// Constructor. /// /// By default this implements the model stored in /// EarthTables/earth_binned_default.txt /// with the detector 6368 km from the center of the Earth, having a /// latitude of 0 deg N and a longitude of 0 deg E. /// /// @param filename - The txt file containing a table of earth layers /// EarthModelBinned::EarthModelBinned(string filename) { SetDetPos(6368, 0, 0); LoadModel(filename); SetRemoveSmallPaths(false); fLonError = 0; } //............................................................................. /// /// Nothing to clean. /// EarthModelBinned::~EarthModelBinned() {} //............................................................................. /// /// Get the set of Earth bins /// /// This returns the set of bins for this Earth model. /// vector EarthModelBinned::GetEarthBins() { return fEarthBins; } //............................................................................. /// /// Clear the earth model information. /// void EarthModelBinned::ClearModel() { fnDepthBins = 0; fnLonBins = 0; fnLatBins = 0; fEarthBins.clear(); } //............................................................................. /// /// Set the coordinates of the detector: /// radius in km, latitude in degrees, longitude in degrees /// Also, updates the detector-coordinate-dependent parts of the trajectory /// constants. /// /// @param rad - The distance from the detector to the Earth's center in km /// @param lat - The latitude of the detector in deg N (between -90 and 90) /// @param lon - The longitude of the detector in deg E (between 0 and 360) /// void EarthModelBinned::SetDetPos(double rad, double lat, double lon) { SetDetectorCoordinates(rad, lat, lon); fC.UpdateDetPos(fDetRadius, fDetLat, fDetLon); } //............................................................................. /// /// Add a bin to the earth model. /// /// @param radius_out - The outer depth of the bin in km /// @param radius_in - The inner depth of the bin in km /// @param latitude - The latitude of the bin center in degrees /// @param longitude - The longitude of the bin center in degrees /// @param density - The density of the matter in the bin in g/cm^3 /// @param zoa - The effective Z/A value of the matter in the bin /// @param index - Region index (all bins with same index have same Z/A) /// 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)); } //............................................................................. /// /// Load an earth model from a file. /// /// By default it loads the model stored in /// EarthTables/earth_binned_default.txt /// /// The row format for the model table needs to be: /// longitude (deg) latitude (deg) outer depth (km) density /// (g/cm^3?) Z/A region_index The data in the table must be in order of /// decreasing depth, followed by increasing longitude, and then increasing /// latitude. Latitude bin widths and longitude bin widths must each be /// constant. In each bin of one coordinate, the bins centers for the other /// coordinates must be the same. /// This assumes that the center of the Earth is at a depth of 6371 km. /// /// @param filename - The txt file containing a table of earth layers /// 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.1) { // 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; } //............................................................................. /// /// Set the effective Z/A value for all bins with given region index. /// /// Use this to change the Z/A of indexed bin, /// e.g. all outer-core layers /// /// @param index - The region index /// @param zoa - The effective Z/A value to use /// 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; } } //............................................................................. /// /// Get the effective Z/A value for all bins in region specified by /// index, e.g. all outer-core layers. /// (Assumes that all bins of given index have same Z/A value) /// /// @param index - The region index /// @return Z/A corresponding to index /// 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; } //............................................................................. /// /// Scale density by scaling factor for all bins in region specified by /// given index. /// /// @param index - The region index /// @param factor - The scaling factor for changing the density /// 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; } } //............................................................................. /// /// Add a path segment to the sequence. /// /// For a given EarthBin, adds a path of a given length in that material /// /// @param length - The length of the path segment in km /// @param bin - The bin we are crossing /// void EarthModelBinned::AddPath(double length, EarthBin bin) { AddPathSegment(length, bin.density, bin.zoa, bin.index); } //............................................................................. /// /// Calculate the distance to the detector along a neutrino trajectory, /// specified by cosT and the azimuthal angle, at the edge of /// the current latitude bin in the direction specified by dir. /// /// In an attempt to make this calculation more robust, this function turns /// the latitude bin iteration around when sin^2(lat) gets within 10^(-14) /// of its extreme limit. /// /// @param cur_index - Index of current Earth bin /// @param L - Information for the lat bin crossings /// @return Calculated detector distance /// double EarthModelBinned::DetDistForNextLatBin(int cur_index, LatBinInfo& L) { // Find next lat bin boundary double NewLat = fEarthBins[cur_index].latitude + L.sign * L.dLat; // Turn around if lat reaches +/- 90 degrees if (abs(NewLat) - 0.5 * M_PI > -1.0e-14) { if (L.maxreached > 0) return -1; // Only turn around once L.maxreached = 1; 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; // Turn around if the radicand is 0 (close enough) or less if (radicand < 1.0e-14) { if (L.maxreached > 0) return -1; // Only turn around once L.maxreached = 1; 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 lon bin } double denominator = fC.gammaSq + NsinSqNewLat; // Calculate distance from the detector at next lat bin boundary double answer = 0; if (denominator == 0) { answer = 0.5 * (fC.rDetCosAcosDetLat / fC.beta - fC.rDetSinDetLat / fC.gamma); } else { answer = -(fC.rDetCosT * NsinSqNewLat + fC.rDetGammaSinDetLat + L.sign * fC.rDetSinT * sinNewLat * sqrt(radicand)) / denominator; } // Indicates that going to the edge of the bin jumped to an invalid part of // the function if (L.maxreached > 0 && answer >= fC.rDetCosAcosDetLat / fC.beta) { return -1; // The rest of the track is contained in this lon bin } return answer; } //............................................................................. /// /// Calculate the distance to the detector at the edge of the current /// longitude bin along a neutrino trajectory specified by cosT and the /// azimuthal angle, setting the index for the next longitude bin in the /// process. The maximum/minimum longitude are expressed around the 0-2pi /// boundary. So, if the track spans that boundary, the actual allowed /// range for lon is [min_lon,2pi)U[0,max_lon]. /// /// @param prev_lon - Value at center of initial longitude bin /// @param L - Information about longitude bin crossings /// @return Calculated detector distance /// 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; } //............................................................................. /// /// Record all the path segments as the neutrino crosses into new latitude /// and longitude bins until it reaches the next depth bin. /// /// All the variables not listed below are the same as the arguments for /// the DetDistForNextLonBin and DetDistForNextLatBin functions. /// /// @param detDist_nextDbin - Distance from the detector at which neutrino /// passes into next depth bin /// @param DetDist - Current distance from the detector /// @param index - Index of current Earth bin /// @param latI - Information about lat bin crossings /// @param lonI - Information about lon bin crossings /// void EarthModelBinned::RecordLatLonBinCrossings(double detDist_nextDbin, double& DetDist, int& index, LatBinInfo& latI, LonBinInfo& lonI) { while (detDist_nextDbin < latI.detDist_nextBin || detDist_nextDbin < lonI.detDist_nextBin) { if (latI.detDist_nextBin > lonI.detDist_nextBin) { // Add segment to next lat bin to path AddPath(DetDist - latI.detDist_nextBin, fEarthBins[index]); // Move lat bins DetDist = latI.detDist_nextBin; index += latI.nextBin - latI.bin; latI.bin = latI.nextBin; // Find next lat bin if (latI.dLat == 0) { latI.detDist_nextBin = -1; } else { latI.detDist_nextBin = DetDistForNextLatBin(index, latI); latI.nextBin = latI.bin + latI.sign; } } else { // 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; // Find next lon bin if (lonI.dLon == 0) { lonI.detDist_nextBin = -1; } else { lonI.detDist_nextBin = DetDistForNextLonBin(fEarthBins[index].longitude, lonI); } } } } //............................................................................. /// /// Find and return longitude bin index that contains longitude. /// /// @param longitude - longitude for which bin is found /// 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); } //............................................................................. /// /// Fill the path sequence in a vector. /// /// This will start at the upper-most layer (where the neutrino's path /// started) and find path lengths to the next bin crossed. When reaching /// the inner-most depth in the given direction, the paths move back to /// outer layers until hitting the detector. The neutrino direction is /// specified as the zenith angle and azimuthal angle of the vector pointing /// from the detector to where the neutrino entered the Earth. /// /// The path sequence is stored as an attribute and can be retrieved with /// the function GetNuPath. /// /// @param cosT - The cosine of the zenith angle of the neutrino direction /// @param phi - The azimuthal angle (in degrees from North) of the neutrino /// direction /// @return The number of path segments in the sequence /// 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 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) * fInvLatBinWidth); int index = fEarthBins.size() - binsPerDepth + init_lonBin * fnLatBins + init_latBin; latinfo.bin = init_latBin; loninfo.bin = init_lonBin; /* 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 = 0; // this changes to 1 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.rDetCosAcosDetLat / fC.beta) { latinfo.sign = -(latinfo.sign); latinfo.maxreached++; } // 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; } /* 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; } /* 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(); }