#include <u.h>
#include <libc.h>

#include "project.h"

static point_xy
tmxy(point_λφ λφ, elliptic *ell, projection *p) {
	point_xy xy;
	double N, Sφ, Tφ, Cφ, A, A⁲, M/*, k*/; /* point dependent */

	λφ.λ -= p->o.λ;
	λφ.φ -= p->o.φ;

	Sφ = sin(λφ.φ);
	Cφ = cos(λφ.φ);
	Tφ = Sφ/Cφ;
	Tφ *= Tφ;
	A = Cφ * λφ.λ;
	A⁲ = A*A;
	A /= sqrt(1. - ell->e⁲*Sφ*Sφ); /* what is this about? */
	M = lenMarc(λφ.φ, Sφ, Cφ, ell);
	N = ell->E⁲ * Cφ*Cφ;

	/* TODO: investigate more terms, and perhaps generalize to arcPell? */
	xy.x = 1./42 * A⁲ * (61. + Tφ * (Tφ * (179. - Tφ) - 479.));
	xy.x = 1./20 * A⁲ * (5. + Tφ * (Tφ - 18.) + N * (14. - 58. * Tφ) + xy.x);
	xy.x = 1./6 * A⁲ * (1. - Tφ + N + xy.x);
	xy.x = 1./1 * A * (1. + xy.x);
	xy.x = p->k₀ * xy.x;

	xy.y = 1./56 *	A⁲ * (1385. + Tφ * (Tφ * (543. - Tφ) - 3111.));
	xy.y = 1./30 *	A⁲ * (61. + Tφ * (Tφ - 58.) + N * (270. - 330 * Tφ) + xy.y);
	xy.y = 1./12 *	A⁲ * (5. - Tφ + N * (9. + 4. * N) + xy.y);
	xy.y = 1./2 *	A * (1. + xy.y);
	xy.y = p->k₀ * (M - p->M₀ + Sφ * λφ.λ * xy.y);

	xy.x = (ell->a * xy.x + p->false.x);
	xy.y = (ell->a * xy.y + p->false.y);

	return xy;
};

static point_λφ
tmλφ(point_xy xy, elliptic *ell, projection *p) {
	double n, con, Cφ, d, D⁲, Sφ, Tφ, φ;
	point_λφ λφ;

	xy.x = (xy.x - p->false.x) / ell->a;
	xy.y = (xy.y - p->false.y) / ell->a;

	φ = arcMlen(p->M₀ + xy.y / p->k₀, ell);

	Sφ = sin(φ);
	Cφ = cos(φ);
	Tφ = tan(φ);
	n = ell->E⁲ * Cφ * Cφ;
	con = 1. - ell->e⁲ * Sφ*Sφ;
	d = xy.x * sqrt(con) / p->k₀;
	con *= Tφ;
	Tφ *= Tφ;
	D⁲ = d * d;

	λφ.φ = 1./56 * D⁲ * (1385. + Tφ * (3633. + Tφ * (4095. + 1574. * Tφ)));
	λφ.φ = 1./30 * D⁲ * (61. + Tφ * (90. - 252. * n + 45. * Tφ) + 46. * n - λφ.φ);
	λφ.φ = 1./12 * D⁲ * (5. + Tφ * (3. - 9. *  n) + n * (1. - 4 * n) - λφ.φ);
	λφ.φ = 1./2 * D⁲ * (1. - λφ.φ);
	λφ.φ = φ - con / (1. - ell->e⁲) * λφ.φ;

	λφ.λ = 1./42 * D⁲ * (61. + Tφ * (662. + Tφ * (1320. + 720. * Tφ)));
	λφ.λ = 1./20 * D⁲ * (5. + Tφ*(28. + 24.*Tφ + 8.*n) + 6.*n - λφ.λ);
	λφ.λ = 1./6 * D⁲ * (1. + 2.*Tφ + n - λφ.λ);
	λφ.λ = d/Cφ * (1 - λφ.λ);

	λφ.λ += p->o.λ;
	λφ.φ += p->o.φ;

	return λφ;
}

int
tmin(projection *p, elliptic *ell, char *){
	p->λφ=tmλφ;
	p->xy=tmxy;
	p->M₀ = lenMarc(p->o.φ, sin(p->o.φ), cos(p->o.φ), ell);
	return 1;
};

int
utmin(projection *p, elliptic *ell, char *arg){
	int zone, south;
	char *c;

	zone = -1; south = 0;
	for(c=arg;c-1!=nil;){
		switch(*c){
		case 's':	c++;
			south=1;
			break;
		case 'z':	c++;
			zone = (strtol(c, &c, 10) - 1) % 60;
			break;
		default:	c++;
		}
		c=strchr(c, ',')+1;
	}

	p->k₀ = 0.9996;
	p->false.x = 500000.;
	p->false.y = south?10000000:0.;
	p->o.λ = ((float)zone + 0.5)*PI/30 - PI;
	p->o.φ = 0;
	return tmin(p, ell, arg);
}