#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/JDetectorSupportkit.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;
 

  /**
   * 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 so-called 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 mounted in 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();

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


      /**
       * 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 ts = 0.5 * 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);

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

	  z -= (ls - length) / (1.0 + ts * (1.0 - mechanics.a*mechanics.b / (1.0 - mechanics.a*z)));
	}

	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 point of the string from
       * which the positions of the piezo sensors and hydrophone are calculated.
       *
       * \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 calculated.
       *
       * 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 given floor.
       *
       * \param  floor                    floor
       * \return                          true if floor present; else false
       */
      bool hasFloor(size_t floor) const
      {
	if (floor == 0)
	  return has_hydrophone;
	else
	  return (floor < this->size());
      }


      /**
       * Get height of given floor.
       *
       * \param  floor                    floor
       * \return                          height
       */
      double getHeight(const size_t floor) const
      {
	if (floor == 0) {

	  if (has_hydrophone) {
	    return hydrophone.getZ();
	  }

	} else if (floor < this->size()) {
	
	  return (*this)[floor].getHeight();
	}
 
	THROW(JValueOutOfRange, "Invalid floor " << floor);
      }


      /**
       * Get position of 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 (floor == 0) {

	  if (has_hydrophone) {
	    return this->getPosition() + hydrophone.getPosition();
	  }

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

	  return this->getPosition() + getPosition(parameters, mechanics, (*this)[floor].getHeight());
	}

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


      /**
       * Get position of given floor according to default string model parameters.
       *
       * \param  floor                    floor
       * \return                          position
       */
      JPosition3D getPosition(const size_t floor) const
      {
	return getPosition(JMODEL::JString(), floor);
      }


      /**
       * 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     = this->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;

    protected:
      /**
       * 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<JDETECTOR::JHydrophone>& hydrophones = std::vector<JDETECTOR::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