/* -*- C++ -*-
 * Copyright ©2006 Hugo Mills <hugo@carfax.org.uk>
 *
 * This software is distributed under the terms of the GNU GPL v3
 * For more information on the GPL, see the file COPYING or 
 * visit http://www.gnu.org/
 *
 * This software is distributed without warranty
 */

#include "xfmap_mercator.h"

extern "C" PluginInfo* getmodules(int i)
{
	if(i==0)
	{
		PluginInfo_Map* pi = new PluginInfo_Map();

		pi->name = "Mercator Projection";
		pi->author = "Hugo Mills";
		pi->copyright = "©2006 Hugo Mills";
		pi->version = 0x10000;
		pi->parser = XfMap_Mercator::parser;
		return pi;
	}
	else
		return NULL;
}

XfMap* XfMap_Mercator::parser(
	ConfigLexer* lex,
	const std::string& type,
	const std::string& subtype
	)
{
	if(type == "mercator")
	{
		while(!lex->eos()) {
			throw(UnknownConfigToken(lex));
		}

		return new XfMap_Mercator();
	}
	else
		return NULL;
}

bool XfMap_Mercator::f(dvec2& result, const dvec3& pos) const
{
	// Derivation: consider a loxodrome of slope k. This is a straight
	// line on the projected map (dq/dp = k). It is a line of constant
	// bearing on the sphere. This line goes a distance k in \phi
	// direction for every unit in \theta direction. At (\theta, \phi)
	// therefore, it goes k d\phi units for each d\theta \cos \phi, so
	// d\phi/d\theta = k \cos \phi = dq/dp \cos \phi. 

	// Now, we can state that \theta = p (given a scaling factor), so
	// we just need to solve q = \integral dq / d\theta d\theta
	//                         = \integral 1/\cos \phi d\phi / d\theta d\theta
	//                         = \integral d\phi/\cos \phi

	// Solve this by substituting x = tan \phi, then x = sinh y

	// p = \theta
	// q = \sinh^{-1} \tan \phi
	// \phi = tan^{-1} \sinh q

	// \theta = atan(z, x)
	// \phi = atan(y, |x, z|)

	double r = sqrt(pos[0]*pos[0] + pos[2]*pos[2]);
	if(r == 0.0) { return false; }
	result[0] = atan2(pos[0], pos[2]) / M_PI;
	double tanphi = pos[1] / r;
	result[1] = asinh(tanphi) / M_PI; // Same scaling as result[0]

	return true;
}

bool XfMap_Mercator::g(dvec3& result, const dvec2& pos) const
{
	// x = \cos \theta \cos \phi
	// y = \sin \phi
	// z = \sin \theta \cos \phi
	double ang = sinh(pos[1]*M_PI);
	result[1] = ang/sqrt(1+ang*ang); // = sin(atan(ang))

	if(pos[0] < -1 || pos[0] > 1)
		return false;

	double m = sqrt(1-result[1]*result[1]);
	double a = pos[0] * M_PI;

	result[0] = m * sin(a);
	result[2] = m * cos(a);

	return true;
}