Removing missed merge conflict indicator.

CMakeLists.txt

@@ -7,11 +7,12 @@
 IF (NOT TRIBITS_PROCESSING_PACKAGE)
   # This CMakeLists.txt file is being processed as the radix TriBITS projects's base
   # CMakeLists.txt file!  (See comments at bottom of this file.)
-  CMAKE_MINIMUM_REQUIRED(VERSION 2.8.11 FATAL_ERROR)
+  CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12 FATAL_ERROR)
   INCLUDE("${CMAKE_CURRENT_SOURCE_DIR}/ProjectName.cmake")
   PROJECT(${PROJECT_NAME} NONE)
   # enable C++11 by default
   SET(radix_ENABLE_CXX11 ON CACHE BOOL "Compile using the C++11 standard" FORCE)
+
   # disable generating HTML dependencies webpage and xml files
   SET(${PROJECT_NAME}_DEPS_XML_OUTPUT_FILE OFF CACHE BOOL "")
   # disable TriBITS export system to save time configuring

radixams/CMakeLists.txt

 TRIBITS_SUBPACKAGE(ams)
 
 SET(SOURCE
+main.cc
 )
 SET(HEADERS
 )
 
 INSTALL(FILES ${HEADERS} DESTINATION "include/radixams/")
+
+TRIBITS_ADD_EXECUTABLE(amscli
+    NOEXEPREFIX
+    NOEXESUFFIX
+    SOURCES ${SOURCE} ${HEADERS}
+)
+
 TRIBITS_SUBPACKAGE_POSTPROCESS()

radixams/cmake/Dependencies.cmake

+IF(NOT USE_QT5)
+   SET(USE_QT4 ON)
+ENDIF()
+
+IF(USE_QT4)
+   SET(QT_PACKAGES QT)
+ELSE()
+   SET(QT_PACKAGES Qt5Widgets)
+ENDIF()
+
 TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
-LIB_REQUIRED_PACKAGES radixmath radixbug
-LIB_OPTIONAL_PACKAGES
-TEST_REQUIRED_PACKAGES
-TEST_OPTIONAL_PACKAGES
-LIB_REQUIRED_TPLS
-LIB_OPTIONAL_TPLS
-TEST_REQUIRED_TPLS
-TEST_OPTIONAL_TPLS
+   LIB_REQUIRED_PACKAGES radixio radixmath radixbug
+   LIB_OPTIONAL_TPLS ${QT_PACKAGES}
 )
+

radixams/main.cc

+#include "radixmath/util.hh"
+
+#include <iomanip>
+
+int main(int argc, char **argv)
+{
+    // expecting three arguments
+    if(argc != 4)
+    {
+        // print usage info
+        std::cout << "Usage: amscli <options>" << std::endl << std::endl
+                  << "Options:" << std::endl
+                  << "  altitude - Altitude at which to calculate dose rate (0 for sea level)" << std::endl
+                  << "  activity - Activity of the area in Ci/m^2" << std::endl
+                  << "  energy   - Gamma ray energy in MeV" << std::endl;
+        return 0;
+    }
+
+    // altitude above sea level - default to 0
+    auto z = 0.0;
+    z = std::stod(argv[1]);
+
+    // activity in Ci/m^2 of the area of interest - default to 1.767e3
+    auto a = 1.767e3;
+    a = std::stod(argv[2]);
+
+    // gamma-ray energy in MeV - default to 7e-1
+    auto e = 7e-1;
+    e = std::stod(argv[3]);
+
+    // radius of area of interest (assumed to be altitude), distance
+    auto r = z;
+    auto s = std::sqrt(z*z + r*r);
+
+    // sigma scalar; constant at sea level
+    auto sigma = 1.225;
+
+    // not at sea level, update sigma
+    if(z != 0)
+    {
+        sigma = 1.911515e-16 * (z * z * z * z)
+              - 9.857582e-13 * (z * z * z)
+              + 5.3399174e-9 * (z * z)
+              - 1.180686e-4  * z
+              + 1.22507;
+    }
+
+    // absorption/attenuation of gamma ray as it passes through air
+    auto mua = radix::gammaRayAbsorptionInAir(e);
+    auto mut = radix::gammaRayAttenuationInAir(e);
+
+    // mu scalars
+    auto mutg = -mut * 1.225;
+    auto muag = -mua * 1.225;
+    auto mutz = -mut * sigma;
+    auto muaz = -mua * sigma;
+
+    // exponential integral values and their ratio
+    auto EIg = radix::exponentialIntegral(mutg);
+    auto EIz = radix::exponentialIntegral(mutz * z);
+    auto EIRatio = EIz / EIg;
+
+    // dose rate at sea level/altitude
+    auto DRg = 0.62 * a;
+    auto DRz = DRg * EIRatio;
+
+    // report dose rate at sea level
+    if(z == 0)
+    {
+        std::cout << std::scientific << DRg << " is the Dose Rate at Ground Sea Level in R/hr" << std::endl;
+    }
+    // report dose rate at altitude above sea level
+    else
+    {
+        std::cout << std::scientific << DRz << " is the Dose Rate at your Altitude above Ground Sea Level in R/hr" << std::endl;
+    }
+
+    return 0;
+}

radixio/CMakeLists.txt

@@ -2,6 +2,7 @@ TRIBITS_SUBPACKAGE(io)
 
 
 SET(SOURCE
+   commandline.cc
    csvfile.cc
    eafstream.cc
    endian.cc
@@ -10,6 +11,8 @@ SET(SOURCE
 )
 
 SET(HEADERS
+   commandline.hh
+   commandline.i.hh
    csvfile.hh
    eafstream.hh
    endian.hh
@@ -21,6 +24,7 @@ SET(HEADERS
 # Add library
 TRIBITS_ADD_LIBRARY(radixio
   SOURCES ${SOURCE}
+  NOINSTALLHEADERS ${HEADERS}
 )
 #
 # Add testing directory

radixio/commandline.cc

+/*
+ * @file:   commandline.cc
+ * @author: Jordan P. Lefebvre, lefebvrejp@ornl.gov
+ *
+ * Created on July 3, 2013, 2:49 PM
+ */
+#include <iostream>
+#include <cstdlib>
+
+#include "radixio/commandline.hh"
+namespace radix {
+
+CommandLine::CommandLine(int argc, char ** argv)
+    : mArgc(argc), mArgv(argv){
+    this->mExecutable = std::string(argv[0]);
+    for (int i = 1; i < argc; i++) {
+        std::string temp = std::string(argv[i]);
+        //
+        // Check for '-' as first character
+        //
+        if (temp.size() > 0) {
+            if (temp.at(0) == '-') {
+                //
+                // Check for flag=value
+                //
+                size_t index = temp.find('=');
+                if (index == std::string::npos) {//we didn't find equal sign
+                    //
+                    // This is just a flag with empty value
+                    //
+                    this->mFlagToValues[temp] = std::string("");
+                } else {
+                    //
+                    // We found equal sign, so we must divide into flag and value
+                    //
+                    std::string flag = temp.substr(0, index);
+                    //
+                    // Check to make sure empty string wasn't provided
+                    // flag=''
+
+                    std::string value = "";
+                    if (temp.size() - 1 != index) {
+                        value = temp.substr(index + 1, temp.size() - 1);
+                    }
+                    this->mFlagToValues.insert(std::make_pair(flag, value));
+                }//else
+            }//if = '-'
+            else {
+                mArgs.push_back(temp);
+            }
+        }//if size > 0
+    }//loop
+}//constructor
+
+CommandLine::CommandLine(const CommandLine& orig)
+    : mArgc(orig.mArgc),
+      mArgv(orig.mArgv),
+      mFlagToValues(orig.mFlagToValues),
+      mArgs(orig.mArgs) {
+}
+
+CommandLine::~CommandLine() {
+}
+
+std::string CommandLine::toString() const {
+    std::stringstream ss;
+    ss << "executable: " << executable() << std::endl;
+    ss << "options: " << std::endl;
+    std::map<std::string, std::string>::const_iterator iter;
+    for(iter=mFlagToValues.begin(); iter!=mFlagToValues.end(); ++iter){
+        ss << iter->first << "=" << iter->second << std::endl;
+    }
+    ss << "arguments: " << std::endl;
+    const std::vector<std::string>& args=arguments();
+    for(size_t i=0; i<args.size(); ++i){
+        ss << (i+1) << ". " << args[i] << std::endl;
+    }
+    return ss.str();
+}
+
+bool CommandLine::exists(std::string flag) const {
+    std::map<std::string, std::string>::const_iterator iter;
+    iter = mFlagToValues.find(flag);
+    if (iter == mFlagToValues.end()) return false;
+    return true;
+}//exists
+
+const std::vector<std::string>& CommandLine::arguments() const {
+    return mArgs;
+}//arguments
+
+int CommandLine::valueAsInt(std::string flag) const {
+    std::map<std::string, std::string>::const_iterator iter;
+    int result = 0;
+    iter = this->mFlagToValues.find(flag);
+    if (iter == this->mFlagToValues.end()) return -1;
+    try {
+        result = std::atoi(iter->second.c_str());
+    } catch (...) {
+        return -1;
+    }
+    return result;
+}//value
+
+float CommandLine::valueAsFloat(std::string flag) const {
+    std::map<std::string, std::string>::const_iterator iter;
+    float result = 0.0f;
+    iter = this->mFlagToValues.find(flag);
+    if (iter == this->mFlagToValues.end()) return -1;
+    try {
+        result = std::atof(iter->second.c_str());
+    } catch (...) {
+        return -1.0f;
+    }
+    return result;
+}//valueAsFloat
+
+std::string CommandLine::value(std::string flag) const {
+    std::map<std::string, std::string>::const_iterator iter;
+    iter = this->mFlagToValues.find(flag);
+    if (iter == this->mFlagToValues.end()) return "";
+    return iter->second;
+}//value
+
+std::string CommandLine::value(std::string flag, std::string default_value) const {
+    std::map<std::string, std::string>::const_iterator iter;
+    iter = this->mFlagToValues.find(flag);
+    if (iter == this->mFlagToValues.end()){
+        return default_value;
+    } else {
+        return iter->second;
+    }
+}//value
+
+std::string CommandLine::executable() const {
+    return this->mExecutable;
+}
+
+//specialize boolean so we can handle yes/true/on and no/false/off
+template<>
+bool CommandLine::valueAs(std::string flag, bool default_value) const{
+
+    if( exists(flag) ){
+        std::string sv=value(flag);
+        if( sv.compare("yes")==0 ||
+                sv.compare("true")==0 ||
+                sv.compare("on")==0 ||
+                sv.compare("1")==0 )return true;
+        if( sv.compare("no")==0 ||
+                sv.compare("false")==0 ||
+                sv.compare("off")==0 ||
+                sv.compare("0")==0 )return false;
+    }
+    return default_value;
+}
+
+}//namespace radix
+

radixio/commandline.hh

+/*
+ * @file:   commandline.hh
+ * @author: Jordan P. Lefebvre, lefebvrejp@ornl.gov
+ *
+ * Created on July 3, 2013, 2:49 PM
+ */
+
+#ifndef RADIX_RADIXIO_COMMANDLINE_HH_
+#define RADIX_RADIXIO_COMMANDLINE_HH_
+#include <vector>
+#include <map>
+#include <string>
+#include <memory>
+#include <sstream>
+
+namespace radix {
+
+    /**
+     * @class CommandLine
+     * @brief Simple class designed to decipher commandline arguments
+     * Provides ability to check flag existences.
+     * Provides ability to get flag=value options
+     * Flags are designated with '-' prefix
+     * All non flag prefixed values are stored as arguments
+     */
+    class CommandLine {
+    public:
+        CommandLine(int argc, char ** argv);
+        CommandLine(const CommandLine& orig);
+        virtual ~CommandLine();
+        /**
+         * Does the flag exist
+         * @param flag
+         * @return
+         */
+        virtual bool exists(std::string flag) const;
+        /**
+         * @brief Get the value of flag as int
+         * @param flag std::string
+         * @return int
+         */
+        int valueAsInt(std::string flag) const;
+        /**
+         * @brief Get the value of flag as float
+         * @param flag std::string
+         * @return float
+         */
+        float valueAsFloat(std::string flag) const;
+        /**
+         * @brief Get the value of flag as std::string
+         * @param flag std::string
+         * @return std::string
+         */
+        std::string value(std::string flag) const;
+        std::string value(std::string flag, std::string default_value) const;
+        /**
+         * @brief Templated get the value of flag with default
+         * @param flag std::string
+         * @param default_value T
+         * @return T
+         */
+        template<class T>
+            T valueAs(std::string flag, T default_value) const;
+
+        /**
+         * @brief Get arguments on the commandline which were not flag=value
+         * @return const std::vector<std::string>&
+         */
+        virtual const std::vector<std::string>& arguments() const;
+        /**
+         * @brief Get the executable path
+         * @return std::string
+         */
+        std::string executable() const;
+        /**
+         * @brief Get the contents as a string
+         * @return std::string
+         */
+        std::string toString() const;
+
+        /**
+         * @brief argc inline argc()
+         * @return int number of arguments
+         */
+        int argc() const { return mArgc; }
+        /**
+         * @brief argv inline argv()
+         * @return char array of commandline arguments
+         */
+        char** argv() const { return mArgv; }
+
+    private:
+
+        int mArgc;
+        char ** mArgv;
+        /** Executable path*/
+        std::string mExecutable;
+        /** Flag to values map*/
+        std::map<std::string, std::string> mFlagToValues;
+        /** List of arguments which were not flag=value*/
+        std::vector<std::string> mArgs;
+
+    };
+
+}//namespace radix
+
+#include "radixio/commandline.i.hh"
+
+#endif	/* RADIX_RADIXIO_COMMANDLINE_HH_ */
+
+

radixio/commandline.i.hh

+#include<sstream>
+
+namespace radix {
+
+    template<class T>
+    T CommandLine::valueAs(std::string flag, T default_value) const{
+        T tvalue;
+        if( exists(flag) && std::istringstream(value(flag)) >> tvalue){
+            return tvalue;
+        } else {
+            return default_value;
+        }
+    }
+
+}//namespace radix
+
+

radixmath/util.cc

@@ -177,4 +177,145 @@ Real greatCircleArea(Real lat1, Real lon1, Real lat2, Real lon2, Real r1)
 
 }
 
+double interpolate(const std::vector<double> &xvals, const std::vector<double> &yvals, double x, bool linear)
+{
+    // iterator to upper bound
+    auto ub = std::upper_bound(xvals.begin(), xvals.end(), x);
+
+    // below range, return first item
+    if(ub == xvals.begin())
+    {
+        return yvals.front();
+    }
+    // above range, return last item
+    else if(ub == xvals.end())
+    {
+        return yvals.back();
+    }
+
+    // indices to right/left side
+    auto r = std::distance(xvals.begin(), ub);
+    auto l = r - 1;
+
+    // min/max vals
+    auto xmin = xvals[l];
+    auto xmax = xvals[r];
+    auto ymin = yvals[l];
+    auto ymax = yvals[r];
+
+    // adjust for log interpolation
+    if(!linear)
+    {
+        xmin = std::log(xmin);
+        xmax = std::log(xmax);
+        ymin = std::log(ymin);
+        ymax = std::log(ymax);
+        x = std::log(x);
+    }
+
+    // scaling ratio
+    auto ratio = (x - xmin) / (xmax - xmin);
+
+    // interpolation result
+    auto result = yvals[l] + (ymax - ymin) * ratio;
+    return result;
+}
+
+double gammaRayAbsorptionInAir(double energy, bool linear)
+{
+    // supported energies in MeV
+    static const std::vector<double> energies =
+    {
+        0.001, 0.0015, 0.002, 0.003,
+        0.004, 0.005,  0.006, 0.008,
+        0.01,  0.015,  0.02,  0.03,
+        0.04,  0.05,   0.06,  0.08,
+        0.1,   0.15,   0.2,   0.3,
+        0.4,   0.5,    0.6,   0.8,
+        1,     1.25,   1.5,   2,
+        3,     4,      5,     6,
+        8,     10,     15,    20
+    };
+
+    // absorption coefficients
+    const std::vector<double> absorption =
+    {
+        3599,    1188,    526.2,   161.4,
+        76.36,   39.31,   22.7,    9.446,
+        4.742,   1.334,   0.5389,  0.1537,
+        0.6833,  0.04098, 0.03041, 0.02407,
+        0.02325, 0.02496, 0.02672, 0.02872,
+        0.02949, 0.02966, 0.02953, 0.02882,
+        0.02789, 0.02666, 0.02547, 0.02345,
+        0.02057, 0.0187,  0.0174,  0.01647,
+        0.01525, 0.0145,  0.01353, 0.01311
+    };
+
+    // interpolate to determine the absorption for the given energy
+    return interpolate(energies, absorption, energy, linear);
+}
+
+double gammaRayAttenuationInAir(double energy, bool linear)
+{
+    // supported energies in MeV
+    static const std::vector<double> energies =
+    {
+        0.001, 0.0015, 0.002, 0.003,
+        0.004, 0.005,  0.006, 0.008,
+        0.01,  0.015,  0.02,  0.03,
+        0.04,  0.05,   0.06,  0.08,
+        0.1,   0.15,   0.2,   0.3,
+        0.4,   0.5,    0.6,   0.8,
+        1,     1.25,   1.5,   2,
+        3,     4,      5,     6,
+        8,     10,     15,    20
+    };
+
+    // attenuation coefficients
+    static const std::vector<double> attenuation =
+    {
+        3606,    1191,    527.9,   162.5,
+        77.88,   40.27,   23.41,   9.921,
+        5.12,    1.614,   0.7779,  0.3538,
+        0.2485,  0.208,   0.1875,  0.1662,
+        0.1541,  0.1356,  0.1233,  0.1067,
+        0.09549, 0.08712, 0.08055, 0.07074,
+        0.06358, 0.05687, 0.05175, 0.04447,
+        0.03581, 0.03079, 0.02751, 0.02522,
+        0.02225, 0.02045, 0.0181,  0.01705
+    };
+
+    // interpolate to determine the attenuation for the given energy
+    return interpolate(energies, attenuation, energy, linear);
+}
+
+double exponentialIntegral(double d)
+{
+    // based on polynomial approximation in document provided by Vince
+    // specialization for d <= 1
+    if(d <= 1)
+    {
+        return 0.00107857 * (d * d * d * d * d)
+             - 0.00976004 * (d * d * d *d)
+             + 0.05519968 * (d * d * d)
+             - 0.24991055 * (d * d)
+             + 0.99999193 * d
+             - 0.57721566
+             - std::log(d);
+    }
+
+    // specialization for d > 1
+    auto numer = 1             * (d * d * d * d)
+               + 8.5733287401  * (d * d * d)
+               + 18.0590169730 * (d * d)
+               + 8.6347608925  * d
+               + 0.2677737343;
+    auto denom = 1             * (d * d * d * d)
+               + 9.5733223454  * (d * d * d)
+               + 25.6329561486 * (d * d)
+               + 21.0996530827 * d
+               + 3.9584969228;
+
+    return numer / denom / (d * std::exp(d));
+}
 } // namespace radix

radixmath/util.hh

@@ -122,6 +122,38 @@ inline Real toTurns( Real degrees )
  *          t == maxReal if no intersection occurs
  */
 Real lineIntersect(const class Point3D& p, const class Vector3D& v, const class Point3D& sp,const class Point3D& ep);
+
+/**
+ * @brief interpolate interpolates the given data vectors
+ * @param xvals - a list of bin boundaries in the x-dimension
+ * @param yvals - a list of y-dimension values at each x-value
+ * @param x - an x-value at which to interpolate its associated y-value
+ * @param linear - whether to interpolate in a linear fashion (false for log)
+ */
+double interpolate(const std::vector<double> &xvals, const std::vector<double> &yvals, double x, bool linear = false);
+
+/**
+ * @brief gammaRayAbsorptionInAir calculates the absorption of gamma rays as they pass through air
+ * @param energy - energy of the gamma ray in MeV
+ * @param linear - whether to interpolate across energies in a linear fashion (false for log)
+ * return double - the absorption of gamma rays of a given energy as they pass through air
+ */
+double gammaRayAbsorptionInAir(double energy, bool linear = false);
+