#pragma once

#include <mars_calc.h>

#include "geomorph.h"
#include "gwutil.h"
#include "linalgadapt.h"
#include "gmregistry.h"
#include "fluidyn.h"
#include "geoworld.h"

namespace geoworld
{
	class VolcanoGM : public HeightMapGeoMorph, public BoundedGeoMorph
	{
	private:
		typedef mars::GaussianFn <float> ConeFn;
		typedef mars::ParabolicFn <float, +1> CalderaFn;
		typedef mars::MidpointDisplacementGrid< float > FlowPerturbation;

		typedef mars::SolutionEquation< float, ConeFn, CalderaFn > ConeCalderaSolnEq;
		typedef mars::BisectMethod <float, ConeCalderaSolnEq, 10> ConeCalderaSolver;
		typedef VectorTag< float >::V3 Vector3;
		typedef VectorTag< float >::V2 Vector2;

		typedef mars::ParabolicFn<float, -1> EaseInFn;								// Used for tapering the flow line
		typedef mars::ParabolicFn<float, -1> EaseOutFn;								// Used for tapering the flow line

		class ConeZone : public geoworld::CylindricalZone <ConeZone>
		{
		private:
			class Intersect
			{
			public:
				const mars::Magnitudinal< float > val;
				const geoworld::TaperFn taper;
				const float outside;

				Intersect ()
					: val(0), outside(0) {}
				Intersect (const mars::Magnitudinal< float > & val, const float fScale, const float fFalloff)
					: val(val), taper(geoworld::TaperFn::createTaperFn(fFalloff * fScale, val)), outside(fFalloff * fScale + val) {}

				inline Intersect & operator = (const Intersect & copy)
				{
					const_cast< mars::Magnitudinal< float > & > (val) = copy.val;
					const_cast< geoworld::TaperFn & > (taper) = copy.taper;
					const_cast< float & > (outside) = copy.outside;
					return *this;
				}
			} _intersect;

			inline float area (const float fUntil) const 
				{ return cone.f2(fUntil) - cone.f2(_intersect.val) + caldera.f2(_intersect.val); }

			inline static Intersect createIntersect (const ConeFn & cone, const CalderaFn & caldera, const float fFalloffFactor)
				{ return Intersect(ConeCalderaSolver ::solve(ConeCalderaSolnEq(cone, caldera), 0, cone.scale *2), cone.scale, fFalloffFactor); }

		public:
			const ConeFn cone;
			const CalderaFn caldera;
			const Vector2 delta;
			const float falloff;

			ConeZone () : falloff(0) {}
			ConeZone (const ConeFn & cone, const CalderaFn & caldera, const float fFalloffFactor, const Vector2 & delta = Vector2())
				: cone(cone), caldera(caldera), delta(delta),  falloff(fFalloffFactor),
				_intersect(createIntersect(cone, caldera, fFalloffFactor))
			{}

			inline ConeZone & operator = (const ConeZone & copy)
			{
				const_cast< ConeFn & > (cone) = copy.cone;
				const_cast< CalderaFn & > (caldera) = copy.caldera;
				const_cast< Vector2 & > (delta) = copy.delta;
				const_cast< float & > (falloff) = copy.falloff;
				_intersect = copy._intersect;
				return *this;
			}
			inline mars::ObjectStream & operator << (mars::ObjectStream & ins)
			{
				float fScale, fConeAmplitude, fConeBroadness, fCalderaBroadness, fCalderaDepth;

				ins >> fScale >> fConeAmplitude >> fConeBroadness >> fCalderaBroadness >> fCalderaDepth 
					>> const_cast< float & > (falloff) 
					>> const_cast< Vector2 & > (delta);
				const_cast< ConeFn & > (cone) = ConeFn(fScale, fConeAmplitude, fConeBroadness);
				const_cast< CalderaFn & > (caldera) = CalderaFn(fScale, fCalderaBroadness, fCalderaDepth);
				_intersect = createIntersect(cone, caldera, falloff);
				return ins;
			}
			inline mars::ObjectStream & operator >> (mars::ObjectStream & outs) const
			{
				outs << cone.scale << cone.amplitude << cone.broadness << caldera.broadness << caldera.depth << falloff << delta;
				return outs;
			}

			inline float taper (const float x, const float y) const 
				{ return _intersect.taper.f(x, y); }

			inline const mars::Magnitudinal< float > & getIntersect() const
				{ return _intersect.val; }

			inline float full (const float frDepth) const
				{ return caldera.f(0) + (cone.f(0) - caldera.f(0)) * frDepth; }
			inline float outside () const { return _intersect.outside; }
			inline float inside () const { return 0; }
			inline float downside (const float x, const float y) const { return 0; }
			inline float upside (const float x, const float y) const { return std::min(cone.f(x, y), caldera.f(x, y)); }
			inline float area () const { return area(_intersect.outside); }
		} _czone;

		GWSurfacePos
			_ptOrigin;
		float _frCalderaFull;
		mars::RangeX< float > 
			_mmfInnerConeRelativeScale,
			_mmfInnerConeBroadness,
			_mmfInnerConeAmplitude, 
			_mmfInnerCalderaBroadness, 
			_mmfrInnerCalderaDepthRatio,
			_mmfrInnerConeDisplaceRatio;

		mars::BBox< long > _bbox;

		typedef std::vector< ConeZone > InnersList;
		InnersList _inners;

		class FluiDynQueryBridge : public fldyn::IQueryBridge <FluiDynQueryBridge, FieldReal, GeoHeightMap::Precision>
		{
		private:
			const IGeoWorldAccessor * _pAccessor;
			const GMScalarField * _pfsDensity;
			const GWSurfacePos _ptOffset;

		public:
			inline FluiDynQueryBridge (const IGeoWorldAccessor * pAccessor, const GWSurfacePos & ptOffset)
				: _pAccessor (pAccessor), _pfsDensity(pAccessor->queryField(FT_Density)), _ptOffset(ptOffset) {}

			inline DensityPrecision density (const GeoHeightMap::Precision x, const GeoHeightMap::Precision y, const GeoHeightMap::Precision z) const
				{ return (*_pfsDensity) (x + _ptOffset.x, y + _ptOffset.y, z); }
			inline DensityPrecision density (const GeoHeightMap::Precision x, const GeoHeightMap::Precision y) const
				{ return (*_pfsDensity) (x + _ptOffset.x, y + _ptOffset.y, (*_pAccessor)(x + _ptOffset.x, y + _ptOffset.y)); }
			inline GeoHeightMap::Precision elevation (const unsigned int x, const unsigned int y) const
				{ return _pAccessor->mean(x + _ptOffset.x, y + _ptOffset.y); }
		};

		class LavaFlowData
		{
		public:
			GeoHeightMap::View & view;
			GeoHeightMap hmbb, hmbbb;

			LavaFlowData (GeoHeightMap::View * pView)
				: view(*pView), hmbb(*pView), hmbbb(*pView) {}
			~LavaFlowData()
				{ view << hmbbb; }
		};

		typedef fldyn::ChannelSet <FluiDynQueryBridge, float, typename GeoHeightMap::Precision, LavaFlowData> ChannelSetType;

		class MyChannelSet : public ChannelSetType
		{
		private:
			class HumpyBevelFn
			{
			private:
				const float _fBL2;

				inline static float createBL2 () { return 1-(1/exp(1.0f)); }

			public:
				const float amplitude;

				HumpyBevelFn ()
					: _fBL2(createBL2()), amplitude(0) {}
				HumpyBevelFn(const float amplitude)
					: _fBL2(createBL2()), amplitude(amplitude) {}

				inline const float f (const float x) const
				{
					using namespace mars;
					return GaussianFn< float >::f(x, 1.0f, 0.5f * GaussianFn< float >::f(x, 1.0f, _fBL2, amplitude, 0.0f), amplitude, 0.0f);
				}
			};

			class FlowBevelFn : public mars::CacheFnProxy< float, HumpyBevelFn >
			{
			private:
				static inline CacheFnProxy createCacheFnProxy (const float fAmplitude, const unsigned short nResolution)
					{ return CacheFnProxy(HumpyBevelFn(fAmplitude), 0, 1, nResolution); }

			public:
				const float amplitude;

				FlowBevelFn (const unsigned short nResolution, const float fAmplitude)
					: CacheFnProxy(createCacheFnProxy(fAmplitude, nResolution)), amplitude(fAmplitude) {}
				FlowBevelFn (const float fMaxRadius, const float fAmplitude)
					: CacheFnProxy(createCacheFnProxy(fAmplitude, static_cast< unsigned short > (ceil(fMaxRadius)))), amplitude(fAmplitude) {}
			} * _pfnBevel;

			class WaveFn : private mars::GaussianFn< float >
			{
			private:
				float _fScale;

				inline static GaussianFn< float > createGaussian (const float fViscosity)
					{ return GaussianFn(1.0f, 1.0f / fViscosity, 1.0f/3.0f, 0.0f); }

			public:
				const float waves, wavesize, perturbance, waveamplitude, viscosity;

				WaveFn ()
					: waves(0), wavesize(0), perturbance(0), waveamplitude(0), viscosity(0), _fScale(0) {}
				WaveFn (
					const float frWaveSize, 
					const float fWaveAmplitude,
					const float fViscosity,
					const float frPerturbance,
					const float fScale = 1.0f 
				)
					: GaussianFn(createGaussian(fViscosity)), viscosity(fViscosity), _fScale(fScale), waves(1.0f / frWaveSize), wavesize(frWaveSize),
					perturbance (frPerturbance), waveamplitude(fWaveAmplitude)
				{}

				// r : Distance from spine (radi)
				// t : The total distance traveled thus-far in the current spline
				// rT : The total radius of the flow
				inline float f (const float r, const float t, const float rT, const FlowPerturbation & perturb) const
				{
					// TODO: Externalize a wave amplitude parameter, right now depends on pre-configured scale
					// TODO: Scale parameter isn't used, always 1.0, causes perturbance to become dependent on actual scale to have any real effect
					return cosf(
						(
							(
								t
								- GaussianFn::f(r / rT) * rT	// Cancels broadness feature of Gaussian
								+ perturb(
									mars::RNDi(r),	// TODO: Verify use of RNDi
									mars::RNDi(t)
								) * perturbance
							) / static_cast< float > (mars::PI) * waves
						) / _fScale
					) * waveamplitude / 2.0f * _fScale;
				}

				mars::ObjectStream & operator >> (mars::ObjectStream & outs) const
					{ return outs << _fScale << waves << wavesize << perturbance << waveamplitude << viscosity; }

				mars::ObjectStream & operator << (mars::ObjectStream & ins)
				{
					ins
						>> _fScale
						>> const_cast< float & > (waves) 
						>> const_cast< float & > (wavesize) 
						>> const_cast< float & > (perturbance) 
						>> const_cast< float & > (waveamplitude) 
						>> const_cast< float & > (viscosity);

					GaussianFn::operator = (createGaussian(viscosity));
					return ins;
				}
			} _fnWave;

			const EaseInFn _fnEaseIn;
			const EaseOutFn _fnEaseOut;
			geoworld::NoiseGenerator< float > _ngenPert;

			inline static EaseInFn createEaseInFn () { return EaseInFn(1, 1.1f, -1.001f); }
			inline static EaseOutFn createEaseOutFn () { return EaseOutFn(1, 1.1f, -1.001f); }

		protected:
			virtual void applyFlowLine (const FlowLineType & fl, LavaFlowData * pData) const;

		public:
			MyChannelSet();
			MyChannelSet(
				const mars::RangeX< float > & mmfBendZenith,
				const unsigned int nBendZenithSteps,
				const unsigned int nCandidateSampleAmt,
				const mars::RangeX< float > & mmfSplineSeg,
				const float fGravitySlopeGrace,
				const float frSnapToGuide,
				const mars::RangeX< float > & mmfFlowLineWidth,
				const mars::RangeX< unsigned int > & mmnKidsPerLevel,
				const float fwStraightness,
				const float fwElevationInfluence,
				const float fwDensityInfluence,
				const float fwChaosInfluence,
				const float fLavaFlowHeight,
				const float fWaveScale,
				const float fPhase,
				const float fAmplitude,
				const float fViscosity,
				const float fPerturbance,
				const float fPertNoise
			);

			~MyChannelSet();

			void init ( const mars::ptr< FluiDynQueryBridge > & pQuery, const unsigned int nLevels, const mars::BBox< float > & bboxSearch );
			void run ( LockedGWSection::View * pView ) const;

			mars::ObjectStream & operator >> (mars::ObjectStream & outs) const;
			mars::ObjectStream & operator << (mars::ObjectStream & ins);

		} _chans;

		unsigned int	_nNumLavaFlows;
		geoworld::NoiseGenerator< float > _noise;
		mars::RangeX< unsigned int > _mmnLevels;

		inline void addInner (
			const float frdScale, 
			const float fAmplitude, 
			const float fConeBroadness, 
			const float fCaldBroadness, 
			const float frCaldDepthRatio, 
			const float fFalloffFactor,
			VectorTag< float >::PC dpc
		)
		{ 
			const float fScale = _czone.caldera.scale * frdScale;
			dpc.p *= _czone.caldera.scale;
			const Vector2 dpos = dpc;

			_inners.push_back(
				ConeZone(
					ConeFn(
						fScale, 
						fAmplitude, 
						fConeBroadness,
						-ConeFn(
							fScale, 
							fAmplitude, 
							fConeBroadness,
							0
						).f(_czone.getIntersect() + mars::MAG(dpos)) / fScale
					), 
					CalderaFn(
						fScale,
						fCaldBroadness,
						fAmplitude * frCaldDepthRatio
					),
					fFalloffFactor,
					dpos
				)
			);
		}

		VolcanoGM ();

	public:
		VolcanoGM (
			const GWSurfacePos & ptOrigin,
			const float frNoise,
			const float fScale, 
			const float fConeBroadness, 
			const float fConeAmplitude, 
			const float fCalderaBroadness, 
			const float frCalderaDepthRatio, 
			const float frCalderaFull,
			const float fConeFalloffFactor,
			const unsigned int nNumInnerCones,
			const mars::RangeX< float > & mmfInnerConeRelativeScale,
			const mars::RangeX< float > & mmfInnerConeBroadness,
			const mars::RangeX< float > & mmfInnerConeAmplitude, 
			const mars::RangeX< float > & mmfInnerCalderaBroadness, 
			const mars::RangeX< float > & mmfrInnerCalderaDepthRatio,
			const mars::RangeX< float > & mmfrInnerConeDisplaceRatio,
			// Lava flow configs
			// TODO: Sort out this mess of parameters, way too many!
			const unsigned int nNumLavaFlows,
			const float fLavaWaveScale,
			const float fLavaWavePhase,
			const float fLavaWaveAmplitude,
			const float fLavaViscosity,
			const float frLavaWavePerturbance,
			const float fLavaFlowHeight,
			// FluiDyn configs
			const mars::RangeX< float > & mmfBendZenith,
			const unsigned short nBendZenithSteps,
			const unsigned short nCandidateSampleAmt,
			const mars::RangeX< float > & mmfSplineSeg,
			const float fGravitySlopeGrace,
			const float frSnapToGuide,
			const mars::RangeX< float > & mmfFlowLineWidth,
			const mars::RangeX< unsigned int > & mmnLevels,
			const mars::RangeX< unsigned int > & mmnKidsPerLevel,
			const float fwStraightness,
			const float fwElevationInfluence,
			const float fwDensityInfluence,
			const float fwChaosInfluence
		);

		virtual void init (SynthesisSession & manager, const IGeoWorldAccessor & accessor, const unsigned short nMaxWidth, const unsigned short nMaxHeight);
		virtual void doMorph (LockedGWSection & section) const;
		virtual mars::BBox< long > getBBox () const { return _bbox; }

		static VolcanoGM * createInstance (mars::ObjectStream & ins);
		mars::ObjectStream & operator >> (mars::ObjectStream & outs) const;
	};

	inline mars::ObjectStream & operator << (VolcanoGM *& pGM, mars::ObjectStream & ins)
	{ 
		pGM = VolcanoGM::createInstance(ins); 
		return ins;
	}

	class VolcanoGMFactory : public GeoMorphFactory, public IPersistentGeoMorphFactory, public ISurfaceBoundedGeoMorphFactory
	{
	private:
		float 
			_frNoise,
			_fConeFalloffFactor,

			_frLavaWavePerturbance,

			_fGravitySlopeGrace,
			_frSnapToGuide,

			_fwStraightness,
			_fwElevationInfluence,
			_fwDensityInfluence,
			_fwChaosInfluence;

		mars::RangeX< float > 
			_mmfScale,
			_mmfConeBroadness,
			_mmfConeAmplitude,
			_mmfCalderaBroadness, 
			_mmfrCalderaDepthRatio, 
			_mmfrCalderaFull,

			_mmfInnerConeRelativeScale,
			_mmfInnerConeBroadness,
			_mmfInnerConeAmplitude, 
			_mmfInnerCalderaBroadness, 
			_mmfrInnerCalderaDepthRatio,
			_mmfrInnerConeDisplaceRatio,

			_mmfLavaWaveScale,
			_mmfLavaWavePhase,
			_mmfLavaWaveAmplitude,
			_mmfLavaViscosity,
			_mmfLavaFlowHeight,

			_mmfBendZenith,
			_mmfSplineSeg,
			_mmfFlowLineWidth;

		mars::RangeX< unsigned int > 
			_mmnLevels,		
			_mmnKidsPerLevel,
			_mmnNumInnerCones,

			_mmnNumLavaFlows;

		unsigned short 
			_nBendZenithSteps,
			_nCandidateSampleAmt;

	public:
		virtual GeoMorph * createRandomInstance(const long x, const long y) const;
		virtual void configure( IConfigGMFactory * pFactoryConfig, const IConfigGMFactory::Settings & settings );
		virtual void save (const GeoMorph * pGM, mars::ObjectStream & outs) const;
		virtual const GeoMorph * restore (mars::ObjectStream & ins) const;
	};
}