JGeometry.hh

Go to the documentation of this file.

#ifndef __JACOUSTICS__JGEOMETRY__
#define __JACOUSTICS__JGEOMETRY__

#include <ostream>
#include <iomanip>
#include <vector>
#include <map>
#include <algorithm>
#include <cmath>

#include "JLang/JException.hh"
#include "JLang/JPredicate.hh"
#include "JLang/JComparator.hh"
#include "JLang/JManip.hh"

#include "JGeometry3D/JPosition3D.hh"

#include "JTools/JHashMap.hh"

#include "JDetector/JLocation.hh"
#include "JDetector/JHydrophone.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JModuleSupportkit.hh"

#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JAcousticsSupportkit.hh"
#include "JAcoustics/JMechanics.hh"
#include "JAcoustics/JModel.hh"


/**
 * \file
 *
 * Acoustic geometries.
 * \author mdejong
 */
namespace JACOUSTICS {}
namespace JPP { using namespace JACOUSTICS; }

namespace JACOUSTICS {

  using JLANG::JValueOutOfRange;
  using JLANG::JNoValue;
  using JGEOMETRY3D::JVector3D;
  using JGEOMETRY3D::JPosition3D;
  using JTOOLS::JHashMap;
  using JDETECTOR::JLocation;
  using JDETECTOR::JHydrophone;
 

  /**
   * Auxiliary namespace to encapsulate different geometries.\n
   */
  namespace JGEOMETRY {

    /**
     * Floor geometry.
     */
    struct JFloor {
      /**
       * Default constructor.
       */
      JFloor() :
        height(0.0)
      {}


      /**
       * Constructor.\n
       * The height refers to the maximal distance from the top of the T-bar of the parent string to the actual sensor.
       *
       * \param  height                   height
       */
      JFloor(const double height) :
        height(height)
      {}


      /**
       * Get height of this floor.\n
       * The height refers to the maximal distance from
       * the fixed reference point of the parent string.
       *
       * \return                          height
       */
      double getHeight() const
      {
        return height;
      }


      /**
       * Write floor parameters to output stream.
       *
       * \param  out                      output stream
       * \param  floor                    floor
       * \return                          output stream
       */
      friend inline std::ostream& operator<<(std::ostream& out, const JFloor& floor)
      {
        return out << FIXED(7,2) << floor.getHeight();
      }


    protected:
      double height;
    };


    /**
     * String geometry.
     *
     * This data structure provides for the implementation of the dynamical geometry of a detector string.\n
     * In this,
     * the position of floor > 0 corresponds to the piezo sensor which is located inside the optical module and \n 
     * the position of floor = 0 to the hydrophone which is optionally mounted on the anchor.\n
     * The reference position of a string corresponds to the top of the T-bar located on the anchor.\n
     * The position of a piezo sensor depends on the tilt of the string and the mechanical model.\n
     * The position of the hydrophone is fixed.
     * Its value is relative to the reference position of the string.
     */
    struct JString :
      public JPosition3D,
      public std::vector<JFloor>
    {
      using JPosition3D::getPosition;
      using JPosition3D::getDistance;
      using std::vector<JFloor>::operator[];

      static constexpr double PRECISION_M = 1.0e-4;      //!< precision of height evaluation [m]
      
      /**
       * Get approximate length of string.
       *
       * \param  parameters               parameters
       * \param  mechanics                mechanics
       * \param  height                   height
       * \return                          length
       */
      static inline double getLength(const JMODEL::JString& parameters,
                                     const JMechanics&      mechanics,
                                     const double           height)
      {
        const double T2 = parameters.getLengthSquared();
        
        const double x  = 1.0 - mechanics.a*height;

        return sqrt(1.0 + T2) * height  +  T2*mechanics.b * log(x)  +  0.5*T2*mechanics.a*mechanics.b*mechanics.b * (1.0/x - 1.0);
      }


      /**
       * Get approximate height of string.
       *
       * \param  parameters               parameters
       * \param  mechanics                mechanics
       * \param  length                   length
       * \param  precision                precision
       * \return                          height
       */
      static inline double getHeight(const JMODEL::JString& parameters,
                                     const JMechanics&      mechanics,
                                     const double           length,
                                     const double           precision = PRECISION_M)
      {
        const size_t MAXIMUM_NUMBER_OF_ITERATIONS = 10;
        
        const double T2 = parameters.getLengthSquared();

        double z = length;    // start value

        for (size_t i = 0; i != MAXIMUM_NUMBER_OF_ITERATIONS; ++i) {

          const double ls = getLength(parameters, mechanics, z)  -  length;

          if (fabs(ls) <= precision) {
            break;
          }

          const double vs = 1.0 - mechanics.a*mechanics.b / (1.0 - mechanics.a*z);

          z -= ls / (1.0 + 0.5*T2 * vs*vs);
        }

        return z;
      }


      /**
       * Get position at given height according to given string model parameters and mechanics.
       *
       * \param  parameters               parameters
       * \param  mechanics                mechanics
       * \param  height                   height
       * \return                          position
       */
      static inline JPosition3D getPosition(const JMODEL::JString& parameters,
                                            const JMechanics&      mechanics,
                                            const double           height)
      {
        const double h1 = height * (1.0 + parameters.vs);
        const double z1 = mechanics.getHeight(h1);

        return JPosition3D(parameters.tx * z1  +  parameters.tx2 * h1*h1,
                           parameters.ty * z1  +  parameters.ty2 * h1*h1,
                           getHeight(parameters, mechanics, h1));
      }


      /**
       * Default constructor.
       */
      JString() :
        has_hydrophone(false)
      {}


      /**
       * Constructor.
       *
       * The given position corresponds to the reference position of the string from
       * which the positions of the piezo sensors and hydrophone are determined.
       *
       * \param  position                 position
       * \param  mechanics                mechanical model parameters
       */
      JString(const JVector3D&  position,
              const JMechanics& mechanics) :
        JPosition3D(position),
        has_hydrophone(false),
        mechanics(mechanics)
      {}


      /**
       * Constructor.
       *
       * The given position corresponds to the reference position of the string from
       * which the positions of the piezo sensors and hydrophone are determined.
       *
       * The template parameter should correspond to a data type which implements the following policy methods.
       * <pre>
       *     int                       %getFloor();
       *     JGEOMETRY3d::JPosition3D  %getPosition();
       * </pre>
       * In this, the position should correspond to the center of the optical module. 
       *
       * Note that the position of the piezo is offset by JDETECTOR::getPiezoPosition with respect to the center of the optical module.\n
       * The position of the hydrophone should separately be set.
       *
       * \param  position                 position
       * \param  mechanics                mechanical model parameters
       * \param  __begin                  begin of optical modules
       * \param  __end                    end   of optical modules
       */
      template<class T>
      JString(const JVector3D&  position,
              const JMechanics& mechanics,
              T __begin, 
              T __end) :
        JPosition3D(position),
        has_hydrophone(false),
        mechanics(mechanics)
      {
        for (T i = __begin; i != __end; ++i) {
          (*this)[i->getFloor()] = this->getDistance(i->getPosition() + JDETECTOR::getPiezoPosition());
        }
      }


      /**
       * Get mechanical model parameters.
       *
       * \return                          mechanical model parameters
       */
      const JMechanics& getMechanics() const
      {
        return mechanics;
      }


      /**
       * Get floor data.
       *
       * \param  floor                    floor number
       * \return                          floor data
       */
      JFloor& operator[](size_t floor)
      {
        if (floor >= this->size()) {
          this->resize(floor + 1);
        }

        return static_cast<std::vector<JFloor>&>(*this)[floor];
      }


      /**
       * Check if this string has receiver at given floor.
       *
       * \param  floor                    floor
       * \return                          true if receiver present; else false
       */
      bool hasFloor(size_t floor) const
      {
        if (floor == 0)
          return has_hydrophone;
        else
          return (floor < this->size());
      }


      /**
       * Get height of receiver at given floor with respect to reference position.
       *
       * \param  floor                    floor
       * \return                          height
       */
      double getHeight(const size_t floor) const
      {
        if (!hasFloor(floor)) {
          THROW(JValueOutOfRange, "Invalid floor " << floor);
        }

        if (floor == 0)
          return hydrophone.getZ();
        else
          return (*this)[floor].getHeight();
      }


      /**
       * Get position of receiver at given floor according to given string model parameters.
       *
       * \param  parameters               parameters
       * \param  floor                    floor
       * \return                          position
       */
      JPosition3D getPosition(const JMODEL::JString& parameters,
                              const size_t           floor) const
      {
        if (!hasFloor(floor)) {
          THROW(JValueOutOfRange, "Invalid floor " << floor);
        }

        if (floor == 0)
          return this->getPosition() + hydrophone.getPosition();
        else
          return this->getPosition() + getPosition(parameters, this->mechanics, (*this)[floor].getHeight());
      }


      /**
       * Get position of receiver at given floor.
       *
       * \param  floor                    floor
       * \return                          position
       */
      JPosition3D getPosition(const size_t floor) const
      {
        if (!hasFloor(floor)) {
          THROW(JValueOutOfRange, "Invalid floor " << floor);
        }

        if (floor == 0)
          return this->getPosition() + hydrophone.getPosition();
        else
          return this->getPosition() + JPosition3D(0.0, 0.0, (*this)[floor].getHeight());
      }


      /**
       * Get distance between given position and floor according to given string model parameters.
       *
       * \param  parameters               parameters
       * \param  position                 position
       * \param  floor                    floor
       * \return                          distance
       */
      double getDistance(const JMODEL::JString& parameters,
                         const JVector3D&       position,
                         const size_t           floor) const
      {
        return this->getPosition(parameters, floor).getDistance(position);
      }


      /**
       * Get model gradient of distance between given position and floor according to given string model parameters.
       *
       * \param  parameters               parameters
       * \param  position                 position
       * \param  floor                    floor
       * \return                          gradient
       */
      JMODEL::JString getGradient(const JMODEL::JString& parameters,
                                  const JVector3D&       position,
                                  const size_t           floor) const
      {
        if (floor == 0) {

          return JMODEL::JString();

        } else if (floor < this->size()) {

          const JPosition3D pos    = this->getPosition(parameters, floor);
          const double      height = (*this)[floor].getHeight();
          const double      h1     = height * (1.0 + parameters.vs);
          const double      z1     = mechanics.getHeight(h1);

          const double tx = parameters.tx;
          const double ty = parameters.ty;
          const double tz = sqrt(1.0 - tx*tx - ty*ty);

          const double dx = pos.getX() - position.getX();
          const double dy = pos.getY() - position.getY();
          const double dz = pos.getZ() - position.getZ();

          const double D  = sqrt(dx*dx + dy*dy + dz*dz);
          const double vw = 1.0  -  mechanics.a * mechanics.b / (1.0 - mechanics.a*h1);
          const double vs = 1.0  +  0.5 * parameters.getLengthSquared() * vw;
          
          return JMODEL::JString(z1 * dx / D  -  height * (tx / tz) * dz / D,
                                 z1 * dy / D  -  height * (ty / tz) * dz / D,
                                 h1*h1 * dx / D,
                                 h1*h1 * dy / D,
                                 height * vw * (tx * dx  +  ty * dy) / D  +  h1 * (dz / vs) / D);

        } else {

          THROW(JValueOutOfRange, "Invalid floor " << floor);
        }       
      }


      /**
       * Write string parameters to output stream.
       *
       * \param  out                      output stream
       * \param  string                   string
       * \return                          output stream
       */
      friend inline std::ostream& operator<<(std::ostream& out, const JString& string)
      {
        using namespace std;

        for (size_t i = 0; i != string.size(); ++i) {
          if (string.hasFloor(i)) {
            out << setw(2)    << i              << ' ' 
                << FIXED(7,3) << string[i]      << " | "
                << string.getPosition(i)        << ' ' 
                << string.mechanics             << endl;
          }
        }

        return out;
      }


      /**
       * Hydrophone.
       *
       * The position of the hydrophone is relative to the reference position of the string.
       */
      JPosition3D hydrophone;
      bool        has_hydrophone;

    private:
      /**
       * Mechanical data.
       */
      JMechanics mechanics;
    };


    /**
     * Detector geometry.
     */
    struct JDetector :
      public JHashMap<int, JString>
    {
      /**
       * Default constructor.
       */
      JDetector()
      {}


      /**
       * Constructor.
       *
       * Note that the positions of the base modules correspond to the reference position of the string.\n
       * As a consequence, a base module (i.e.\ floor = 0) is required for each string.\n
       * Missing base modules should therefore be added beforehand (e.g.\ using application JDetectorDB.cc).
       *
       * Note that if the position of a hydrophone is not available,
       * it is assumed that there is no hydrophone on that string.\n
       * If the position of the hydrophone is manually set,
       * the corresponding parameter JString::has_hydrophone should be set to <tt>true</tt>.
       *
       * \param  detector                 detector
       * \param  hydrophones              container with data of hydrophones
       */
      JDetector(const JDETECTOR::JDetector&     detector,
                const std::vector<JHydrophone>& hydrophones = std::vector<JHydrophone>())
      {
        using namespace std;
        using namespace JPP;

        map<int, vector<module_type> > buffer;                            // string -> modules

        for (JDETECTOR::JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
          buffer[module->getString()].push_back(module_type(module->getLocation(), module->getPosition()));
        }

        for (map<int, vector<module_type> >::iterator i = buffer.begin(); i != buffer.end(); ++i) {

          sort(i->second.begin(), i->second.end(), make_comparator(&module_type::getFloor));

          if (i->second[0].getFloor() == 0) {

            vector<module_type>::iterator p = i->second.begin();

            ++p;

            (*this)[i->first] = JGEOMETRY::JString(i->second[0].getPosition(), getMechanics(i->first), p, i->second.end());

            try {

              (*this)[i->first].hydrophone     = getPosition(hydrophones.begin(),
                                                             hydrophones.end(),
                                                             make_predicate(&JHydrophone::getString, i->first));

              (*this)[i->first].has_hydrophone = true;
            }
            catch(const exception&) {
              (*this)[i->first].has_hydrophone = false;
            }

          } else {

            THROW(JNoValue, "No floor 0 in string " << i->first << "; use e.g. JDetectorDB -W.");
          }
        }
      }


      /**
       * Check if this detector has given string.
       *
       * \param  string                   string
       * \return                          true if string present; else false
       */
      bool hasString(int string) const
      {
        return this->has(string);
      }


      /**
       * Check if this detector has given location.
       *
       * \param  location                 location
       * \return                          true if location present; else false
       */
      bool hasLocation(const JLocation& location) const
      {
        return this->hasString(location.getString()) && (*this)[location.getString()].hasFloor(location.getFloor());
      }


      /**
       * Write detector parameters to output stream.
       *
       * \param  out                      output stream
       * \param  detector                 detector
       * \return                          output stream
       */
      friend inline std::ostream& operator<<(std::ostream& out, const JDetector& detector)
      {
        using namespace std;

        for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
          out << setw(4) << i->first << endl << i->second;
        }

        return out;
      }


      /**
       * Auxiliary data structure for module location and position.
       */
      struct module_type :
        public JLocation,
        public JPosition3D
      {
        /**
         * Constructor.
         *
         * \param  location           module location
         * \param  position           module position
         */
        module_type(const JLocation&   location,
                    const JPosition3D& position) :
          JLocation  (location),
          JPosition3D(position)
        {}


        /**
         * Less-than operator.
         *
         * \param  first              first  module
         * \param  second             second module
         * \return                    true if floor of first module less than that of second; else false
         */
        friend inline bool operator<(const module_type& first, const module_type& second)
        {
          return first.getFloor() < second.getFloor();
        }
      };
    };
  }


  /**
   * Type definition of detector geometry.
   */
  typedef JGEOMETRY::JDetector  JGeometry;
}

#endif