#pragma once

#include <mars_util.h>

#include "gwutil.h"
#include "gmregistry.h"
#include "geomorph.h"

namespace geoworld
{
	class CraterGM : public HeightMapGeoMorph, public BoundedGeoMorph, public PersistentGeoMorph
	{
	private:
		typedef mars::GaussianFn <float> CentralUpliftEq;
		typedef mars::ParabolicFn <float, +1> BowlFn;
		typedef mars::InvertedSquareFn <float, +1> PerimeterUpsideEq;
		typedef mars::InvertedSquareFn <float, -1> PerimeterDownsideEq;

		class CraterZone : public CylindricalZone <CraterZone>
		{
		public:
			const BowlFn bowl;
			const float radius;

			inline CraterZone (const BowlFn & bowl)
				: bowl(bowl), radius(mars::NewtonsMethod <float, BowlFn, 10>::solve(bowl, bowl.scale))
			{}

			inline float outside () const { return radius; }
			inline float inside () const { return 0; }
			inline float downside (const float x, const float y) const { return bowl.f(x, y); }
			inline float upside (const float x, const float y) const { return 0; }
			inline float area () const { return bowl.f2(radius); }
		};
	
		class SurfaceDensityZoneFn : public CylindricalZone <SurfaceDensityZoneFn>
		{
		private:
			mars::GaussianFn <float> _fn;
			float _fOutside;

		public:
			SurfaceDensityZoneFn (const BowlFn & bowl) // HACK: Export constants, 25 makes it high, B * 11/32 makes it wide enough for the crater-rim
				: _fn(bowl.scale, 25, bowl.broadness / 32 * 11, 0), _fOutside(bowl.scale * bowl.broadness)
			{}

			inline float inside() const { return 0; }
			inline float outside() const { return _fOutside; }
			inline float upside(const float x, const float y) const { return _fn.f(x, y); }
			inline float upsideGrad (const float x, const float y) const { return _fn.df(x); }
			inline float downside(const float x, const float y) const { return 0; }
			inline float area () const { return _fn.f2(_fOutside); }
		};

		class DensityFieldBase
		{
		protected:
			const GWCoords & _gcOrigin;
			const SurfaceDensityZoneFn & _zone;

		public:
			DensityFieldBase (const SurfaceDensityZoneFn & zone, const GWCoords & gcOrigin)
				: _zone(zone), _gcOrigin(gcOrigin) {}
		};

		class DensityGradientField : public GMGradientField, public DensityFieldBase
		{
		public:
			DensityGradientField(const SurfaceDensityZoneFn & zone, const GWCoords & gcOrigin)
				: DensityFieldBase(zone, gcOrigin) {}

			ResultType operator () (const IndexType x, const IndexType y, const IndexType z) const;
		};

		class DensityScalarField : public GMScalarField, public DensityFieldBase
		{
		public:
			DensityScalarField(const SurfaceDensityZoneFn & zone, const GWCoords & gcOrigin)
				: DensityFieldBase(zone, gcOrigin) {}

			const GMGradientField * queryGradient() const { return new DensityGradientField(_zone, _gcOrigin); }
			void releaseGradient (const GMGradientField * pGF) const { delete pGF; }

			ResultType operator () (const IndexType x, const IndexType y, const IndexType z) const;
		};

		class CentralUpliftZone : public CylindricalZone <CentralUpliftZone>
		{
		private:
			typedef mars::SolutionEquation <float, CentralUpliftEq, BowlFn > SolutionEq;
			CentralUpliftEq _up;
			BowlFn _bowl;
			float _fOutside;

		public:
			// FIXME: Gaussians are not configured right yet
			inline CentralUpliftZone (const CentralUpliftEq & central, const BowlFn & bowl)
				: _up (central), _bowl(bowl), 
				_fOutside (
					mars::BisectMethod <float, SolutionEq >
					::solve(
						SolutionEq(central, bowl), 
						0, 
						mars::NewtonsMethod <float, BowlFn>
						::solve(
							bowl, 
							bowl.scale
						)
					)
				)
			{}

			inline float outside () const { return _fOutside; }
			inline float inside () const { return 0; }
			inline float upside (const float x, const float y) const { return _up.f(x, y); }
			inline float downside (const float x, const float y) const { return _bowl.f(x, y); }
			inline float area () const { return _bowl.f2(_fOutside) - _up.f2(_fOutside); }
		};

		class PerimiterUpliftZone : public CylindricalZone <PerimiterUpliftZone>
		{
		public:
			const class OutsideSlope 
			{
			public:
				const mars::InvertedSquareFn< float, +1 >	fn;
				const mars::Magnitudinal< float >			crest;
				const float									outside, falloff;

				const TaperFn taper;

				inline OutsideSlope (const mars::InvertedSquareFn< float, +1 > & fn, const mars::Magnitudinal< float > & crest, const float fSpallFalloffFactor)
					: fn(fn), crest(crest), taper(TaperFn::createTaperFn(crest * fSpallFalloffFactor, crest)), outside(crest * fSpallFalloffFactor + crest), falloff(fSpallFalloffFactor) {}
			} outslope;

		private:
			const mars::InvertedSquareFn< float, -1 >	_downside;
			const float _fInside, _fBowlRadius;
			BowlFn _bowl;

			inline static mars::Magnitudinal<float> solveBowl2Upside (const BowlFn & bowl, const PerimeterUpsideEq & upside)
			{
				using namespace mars;

				return NewtonsMethod <float, SolutionEquation <float, BowlFn, InvertedSquareFn <float, +1> > > 
					::solve (
						SolutionEquation <float, BowlFn, InvertedSquareFn <float, +1> > (
							bowl, 
							upside
						), 
						bowl.scale * 2
					);
			}
			inline static float solveBowl2Downside( const PerimeterDownsideEq & downside, const CraterZone & czone)
			{
				return mars::NewtonsMethod <float, mars::SolutionEquation <float, BowlFn, mars::InvertedSquareFn <float, -1> > > 
					::solve (
						mars::SolutionEquation <float, BowlFn, mars::InvertedSquareFn <float, -1> > (
							czone.bowl, 
							downside
						), 
						czone.radius
					);
			}

		public:
			inline PerimiterUpliftZone (const CraterZone & czone, const PerimeterUpsideEq & upside, const PerimeterDownsideEq & downside, const float fSpallFalloffFactor)
				: _bowl(czone.bowl), 
				outslope(
					upside, 
					solveBowl2Upside(czone.bowl, upside),
					fSpallFalloffFactor
				), 
				_downside(downside),
				_fInside ( solveBowl2Downside(downside, czone) ),
				_fBowlRadius(czone.radius)
			{
			}

			inline float compareToCrest (const float x, const float y) const { return outslope.crest.compareTo(x, y); }

			inline float outside () const { return outslope.outside; }
			inline float inside () const { return _fInside; }
			inline float crest () const { return outslope.crest; }
			inline float downside (const float x, const float y) const { return _downside.f(x, y); }
			inline float upside (const float x, const float y) const { return outslope.crest.compareTo (x, y) > 0 ? outslope.fn.f(x, y) : _bowl.f(x, y); }

			// The most Newton-method-trigger-happy-intersection-routine currently on the planet
			const float area () const
			{
				return 
					// Above the x-axis
					((outslope.fn.f2(outslope.outside) - outslope.fn.f2(outslope.crest)) + (_bowl.f2(outslope.crest) - _bowl.f2(_fBowlRadius)))
					+
					// Below the x-axis
					_downside.f2(outslope.outside) - _downside.f2(_fInside) - (_bowl.f2(_fBowlRadius) - _bowl.f2(_fInside));
			}
		};

		class Zones
		{
		public:
			const PerimiterUpliftZone puz;
			const CentralUpliftZone central;
			const CraterZone czone;
			SurfaceDensityZoneFn density;

			inline Zones (const CraterZone & czone, const PerimeterUpsideEq & perimUp, const PerimeterDownsideEq & perimDown, const CentralUpliftEq & central, const float fSpallFalloffFactor)
				: czone(czone), puz(czone, perimUp, perimDown, fSpallFalloffFactor), central(central, czone.bowl), density(czone.bowl) {}
		} _zones;

		const GWSurfacePos _ptOrigin;
		GWCoords _gcOrigin;
		
		const float 
			_fCentralUpliftDepthScale,
			_fCentralUpliftBroadnessScale;

		const NoiseGenerator< float > _noise;

	public:
		class CraterGMState
		{
		public:
			float frScale, frBowlBroadness, fBowlDepth, fCentralUpliftBroadnessScale, fCentralUpliftDepthScale, fSpallFalloffFactor, frNoise;
			GWSurfacePos ptOrigin;
			GWCoords gcOrigin;

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

		CraterGM (const CraterGMState & state);

		CraterGM(
			const VectorTag< GeoHeightMap::Precision >::V2 & ptOrigin, 
			const float frSample, 
			const float frScale, 
			const float frBowlBroadness, 
			const float fBowlDepth, 
			const float fSpallFalloffFactor, 
			const float frNoise, 
			const float fCentralUpliftDepthScale, 
			const float fCentralUpliftBroadnessScale
		);

		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;
		virtual GMScalarField * queryField (const FieldType enft, const GeoHeightMap & hm);
		virtual void releaseField (const GMScalarField * pField);

		virtual mars::ptr< State > createState () const;
	};

	class CraterGMFactory : public GeoMorphFactory, public ISurfaceBoundedGeoMorphFactory, public IPersistentGeoMorphFactory
	{
	private:
		mars::RangeX< float >
			_rngfrSample, 
			_rngfrScale, 
			_rngfrBowlBroadness, 
			_rngfBowlDepth;

		float
			_fSpallFalloffFactor, 
			_frNoise, 
			_fCentralUpliftDepthScale, 
			_fCentralUpliftBroadnessScale;

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