<?php
class Sourcemap_Proj_Transform_Somerc extends Sourcemap_Proj_Transform {
/*******************************************************************************
NAME                       SWISS OBLIQUE MERCATOR

PURPOSE:	Swiss projection.
WARNING:  X and Y are inverted (weird) in the swiss coordinate system. Not
   here, since we want X to be horizontal and Y vertical.

ALGORITHM REFERENCES
1. "Formules et constantes pour le Calcul pour la
 projection cylindrique conforme à axe oblique et pour la transformation entre
 des systèmes de référence".
 http://www.swisstopo.admin.ch/internet/swisstopo/fr/home/topics/survey/sys/refsys/switzerland.parsysrelated1.31216.downloadList.77004.DownloadFile.tmp/swissprojectionfr.pdf

*******************************************************************************/

    public function init() {
        $phy0 = $this->lat0;
        $this->lambda0 = $this->long0;
        $sinPhy0 = sin($phy0);
        $semiMajorAxis = $this->a;
        $invF = $this->rf;
        $flattening = 1 / $invF;
        $e2 = 2 * $flattening - pow($flattening, 2);
        $e = $this->e = sqrt($e2);
        $this->R = $semiMajorAxis * sqrt(1 - $e2) / (1 - $e2 * pow($sinPhy0, 2.0));
        $this->alpha = sqrt(1 + $e2 / (1 - $e2) * pow(cos($phy0), 4.0));
        $this->b0 = asin($sinPhy0 / $this->alpha);
        $this->K = log(tan(Sourcemap_Proj::PI / 4.0 + $this->b0 / 2.0))
            - $this->alpha
            * log(tan(Sourcemap_Proj::PI / 4.0 + $phy0 / 2.0))
            + $this->alpha
            * $e / 2
            * log((1 + $e * $sinPhy0)
                    / (1 - $e * $sinPhy0));
    }

    public function forward($p) {
        $Sa1 = log(tan(Sourcemap_Proj::PI / 4.0 - $p->y / 2.0));
        $Sa2 = $this->e / 2.0
            * log((1 + $this->e * sin($p->y))
                    / (1 - $this->e * sin($p->y)));
        $S = -$this->alpha * ($Sa1 + $Sa2) + $this->K;

        // spheric latitude
        $b = 2.0 * (atan(exp($S)) - Sourcemap_Proj::PI / 4.0);

        // spheric longitude
        $I = $this->alpha * ($p->x - $this->lambda0);

        // psoeudo equatorial rotation
        $rotI = atan(sin($I)
                / (sin($this->b0) * tan($b) +
                    cos($this->b0) * cos($I)));

        $rotB = asin(cos($this->b0) * sin($b) -
                sin($this->b0) * cos($b) * cos($I));

        $p->y = $this->R / 2.0
            * log((1 + sin($rotB)) / (1 - sin($rotB)))
            + $this->y0;
        $p->x = $this->R * $rotI + $this->x0;
        return $p;
    }

    public function inverse($p) {
        $Y = $p->x - $this->x0;
        $X = $p->y - $this->y0;

        $rotI = $Y / $this->R;
        $rotB = 2 * (atan(exp(X / $this->R)) - Sourcemap_Proj::PI / 4.0);

        $b = asin(cos($this->b0) * sin($rotB)
                + sin($this->b0) * cos($rotB) * cos($rotI));
        $I = atan(sin($rotI)
                / (cos($this->b0) * cos($rotI) - sin($this->b0)
                    * tan($rotB)));

        $lambda = $this->lambda0 + $I / $this->alpha;

        $S = 0.0;
        $phy = $b;
        $prevPhy = -1000.0;
        $iteration = 0;
        while (abs($phy - $prevPhy) > 0.0000001) {
            if (++$iteration > 20) {
                throw new Exception("Infinity...");
            }
            //S = log(tan(Sourcemap_Proj::PI / 4.0 + $phy / 2.0));
            $S = 1.0 / $this->alpha * (log(tan(Sourcemap_Proj::PI / 4.0 + $b / 2.0)) - $this->K) + $this->e * log(tan(Sourcemap_Proj::PI / 4.0 + asin($this->e * sin($phy)) / 2.0));
            $prevPhy = $phy;
            $phy = 2.0 * atan(exp($S)) - Sourcemap_Proj::PI / 2.0;
        }

        $p->x = $lambda;
        $p->y = $phy;
        return $p;
    }
}
