Stereographic inverse equations--mapping x,y to lat/long

@param Point $p
@return Point
178877
public function inverse($p)
    {
        $x = ($p->x - $this->x0) / $this->a;   /* descale and de-offset */
        $y = ($p->y - $this->y0) / $this->a;

        /*
        $lon;
        $lat;
        $cosphi;
        $sinphi;
        $rho;
        $tp = 0.0;
        $phi_l = 0.0;
        $i;
        */

        $halfe = 0.0;
        $pi2 = 0.0;

        if ($this->sphere) {
            /*
            $c;
            $rh;
            $sinc;
            $cosc;
            */

            $rh = sqrt( $x * $x + $y * $y );
            $c = 2. * atan( $rh / $this->akm1 );
            $sinc = sin( $c );
            $cosc = cos( $c );
            $lon = 0.;

            switch ($this->mode) {
                case $this->EQUIT:
                    if (abs( $rh ) <= Common::EPSLN) {
                        $lat = 0.;
                    } else {
                        $lat = asin( $y * $sinc / $rh );
                    }
                    if ($cosc != 0. || $x != 0.0)
                        $lon = atan2( $x * $sinc, $cosc * $rh );
                    break;

                case $this->OBLIQ:
                    if (abs($rh) <= Common::EPSLN) {
                        $lat = $this->phi0;
                    } else {
                        $lat = asin( $cosc * $this->sinph0 + $y * $sinc * $this->cosph0 / $rh );
                    }
                    $c = $cosc - $this->sinph0 * sin( $lat );
                    if ($c != 0. || $x != 0.0) {
                        $lon = atan2( $x * $sinc * $this->cosph0, $c * $rh );
                    }
                    break;

                case $this->N_POLE:
                    $y = -$y;

                case $this->S_POLE:
                    if (abs($rh) <= Common::EPSLN) {
                        $lat = $this->phi0;
                    } else {
                        $lat = asin( $this->mode == $this->S_POLE ? -$cosc : $cosc  );
                    }
                    $lon = ($x == 0. && $y == 0.0) ? 0. : atan2($x, $y);
                    break;
            }

            $p->x = Common::adjust_lon( $lon + $this->long0 );
            $p->y = $lat;
        } else {
            $rho = sqrt($x * $x + $y * $y);

            switch ($this->mode) {
                case $this->OBLIQ:
                case $this->EQUIT:
                    $tp = 2. * atan2( $rho * $this->cosX1, $this->akm1 );
                    $cosphi = cos( $tp );
                    $sinphi = sin( $tp );
                    if ($rho == 0.0) {
                        $phi_l = asin($cosphi * $this->sinX1);
                    } else {
                        $phi_l = asin($cosphi * $this->sinX1 + ($y * $sinphi * $this->cosX1 / $rho));
                    }

                    $tp = tan( .5 * (Common::HALF_PI + $phi_l) );
                    $x *= $sinphi;
                    $y = $rho * $this->cosX1 * $cosphi - $y * $this->sinX1 * $sinphi;
                    $pi2 = Common::HALF_PI;
                    $halfe = .5 * $this->e;
                    break;

                case $this->N_POLE:
                    $y = -$y;

                case $this->S_POLE:
                    $tp = - $rho / $this->akm1;
                    $phi_l = Common::HALF_PI - 2. * atan( $tp );
                    $pi2 = -Common::HALF_PI;
                    $halfe = -.5 * $this->e;
                    break;
            }

            for ($i = $this->NITER; $i--; $phi_l = $lat) { //check this
                $sinphi = $this->e * sin($phi_l);
                $lat = 2.0 * atan( $tp * pow((1.0 + $sinphi) / (1.0 - $sinphi), $halfe)) - $pi2;

                if (abs( phi_l - lat ) < $this->CONV) {
                    if ($this->mode == $this->S_POLE) {
                        $lat = -$lat;
                    }

                    $lon = ($x == 0. && $y == 0.0) ? 0. : atan2($x, $y);
                    $p->x = Common::adjust_lon( $lon + $this->long0);
                    $p->y = $lat;

                    return $p;
                }
            }
        }
    }
218184
function(p) {
        var v = this.b.sub(this.a);
        var t = p.sub(this.a).dot(v) / v.magSquared();
        // Check to see if t is beyond the extents of the line segment
        if (t < 0.0) {
            return this.a.copy();
        } else if (t > 1.0) {
            return this.b.copy();
        }
        // Return the point between 'a' and 'b'
        return this.a.add(v.scaleSelf(t));
    }
1063281
public function inverse( $p ) {
        
        $p->x -= $this->x0;
        $p->y -= $this->y0;

        $lon = Common::adjust_lon( $this->long0 + $p->x / $this->a );
        $lat = 2.5 * (atan( exp( 0.8 * $p->y / $this->a ) ) - Common::PI / 4.0);

        $p->x = $lon;
        $p->y = $lat;
        
        return $p;
    }
492787
function inverse$12(p) {
        p.x -= this.x0;
        p.y -= this.y0;
        var x = p.x / this.a;
        var y = p.y / this.a;
        var lam, phi, cCe, sCe, q, rho, ab;
        if (this.sphere) {
            var cosz = 0,
                rh, sinz = 0;

            rh = Math.sqrt(x * x + y * y);
            phi = rh * 0.5;
            if (phi > 1) {
                return null;
            }
            phi = 2 * Math.asin(phi);
            if (this.mode === this.OBLIQ || this.mode === this.EQUIT) {
                sinz = Math.sin(phi);
                cosz = Math.cos(phi);
            }
            switch (this.mode) {
                case this.EQUIT:
                    phi = (Math.abs(rh) <= EPSLN) ? 0 : Math.asin(y * sinz / rh);
                    x *= sinz;
                    y = cosz * rh;
                    break;
                case this.OBLIQ:
                    phi = (Math.abs(rh) <= EPSLN) ? this.phi0 : Math.asin(cosz * this.sinph0 + y * sinz * this.cosph0 / rh);
                    x *= sinz * this.cosph0;
                    y = (cosz - Math.sin(phi) * this.sinph0) * rh;
                    break;
                case this.N_POLE:
                    y = -y;
                    phi = HALF_PI - phi;
                    break;
                case this.S_POLE:
                    phi -= HALF_PI;
                    break;
            }
            lam = (y === 0 && (this.mode === this.EQUIT || this.mode === this.OBLIQ)) ? 0 : Math.atan2(x, y);
        }
        else {
            ab = 0;
            if (this.mode === this.OBLIQ || this.mode === this.EQUIT) {
                x /= this.dd;
                y *= this.dd;
                rho = Math.sqrt(x * x + y * y);
                if (rho < EPSLN) {
                    p.x = 0;
                    p.y = this.phi0;
                    return p;
                }
                sCe = 2 * Math.asin(0.5 * rho / this.rq);
                cCe = Math.cos(sCe);
                x *= (sCe = Math.sin(sCe));
                if (this.mode === this.OBLIQ) {
                    ab = cCe * this.sinb1 + y * sCe * this.cosb1 / rho;
                    q = this.qp * ab;
                    y = rho * this.cosb1 * cCe - y * this.sinb1 * sCe;
                }
                else {
                    ab = y * sCe / rho;
                    q = this.qp * ab;
                    y = rho * cCe;
                }
            }
            else if (this.mode === this.N_POLE || this.mode === this.S_POLE) {
                if (this.mode === this.N_POLE) {
                    y = -y;
                }
                q = (x * x + y * y);
                if (!q) {
                    p.x = 0;
                    p.y = this.phi0;
                    return p;
                }
                ab = 1 - q / this.qp;
                if (this.mode === this.S_POLE) {
                    ab = -ab;
                }
            }
            lam = Math.atan2(x, y);
            phi = authlat(Math.asin(ab), this.apa);
        }

        p.x = adjust_lon(this.long0 + lam);
        p.y = phi;
        return p;
    }
509350
public function forward($p)
    {
        $lon = $p->x;
        $lat = $p->y;

        $delta_lat = $lat - $this->lat0;
        $delta_lon = $lon - $this->long0;

        // 1. Calculate d_phi and d_psi    ...                          // and d_lambda
        // For this algorithm, delta_latitude is in seconds of arc x 10-5, so we need to scale to those units. Longitude is radians.
        $d_phi = $delta_lat / Common::SEC_TO_RAD * 1E-5;
        $d_lambda = $delta_lon;
        $d_phi_n = 1;  // d_phi^0

        $d_psi = 0;

        for ($n = 1; $n <= 10; $n++) {
            $d_phi_n = $d_phi_n * $d_phi;
            $d_psi = $d_psi + $this->A[$n] * $d_phi_n;
        }

        // 2. Calculate theta
        $th_re = $d_psi;
        $th_im = $d_lambda;

        // 3. Calculate z
        $th_n_re = 1;
        $th_n_im = 0;  // theta^0
        #$th_n_re1;
        #$th_n_im1;

        $z_re = 0;
        $z_im = 0;

        for ($n = 1; $n <= 6; $n++) {
            $th_n_re1 = $th_n_re * $th_re - $th_n_im * $th_im;
            $th_n_im1 = $th_n_im * $th_re + $th_n_re * $th_im;
            $th_n_re = $th_n_re1;
            $th_n_im = $th_n_im1;
            $z_re = $z_re + $this->B_re[$n] * $th_n_re - $this->B_im[$n] * $th_n_im;
            $z_im = $z_im + $this->B_im[$n] * $th_n_re + $this->B_re[$n] * $th_n_im;
        }

        // 4. Calculate easting and northing
        $p->x = ($z_im * $this->a) + $this->x0;
        $p->y = ($z_re * $this->a) + $this->y0;

        return $p;
    }
265536
public function geocentric_to_wgs84(Point $p)
    {
        Proj4php::reportDebug(sprintf(
            "geocentric_to_wgs84(%s,%s)\n",
            $p->x,
            $p->y
        ));

        if ($this->datum_type == Common::PJD_3PARAM) {
            Proj4php::reportDebug("+x=".$this->datum_params[0]."\n");
            Proj4php::reportDebug("+y=".$this->datum_params[1]."\n");
            Proj4php::reportDebug("+z=".$this->datum_params[2]."\n");
            $p->x += $this->datum_params[0];
            $p->y += $this->datum_params[1];
            $p->z += $this->datum_params[2];
        } elseif ($this->datum_type == Common::PJD_7PARAM) {
            Proj4php::reportDebug("Dx=".$this->datum_params[0]."\n");
            Proj4php::reportDebug("Dy=".$this->datum_params[1]."\n");
            Proj4php::reportDebug("Dz=".$this->datum_params[2]."\n");
            Proj4php::reportDebug("Rx=".$this->datum_params[3]."\n");
            Proj4php::reportDebug("Ry=".$this->datum_params[4]."\n");
            Proj4php::reportDebug("Rz=".$this->datum_params[5]."\n");
            Proj4php::reportDebug("M=".$this->datum_params[6]."\n"); 
            $Dx_BF = $this->datum_params[0];
            $Dy_BF = $this->datum_params[1];
            $Dz_BF = $this->datum_params[2];
            $Rx_BF = $this->datum_params[3];
            $Ry_BF = $this->datum_params[4];
            $Rz_BF = $this->datum_params[5];
            $M_BF = $this->datum_params[6];

            $p->x = $M_BF * ($p->x - $Rz_BF * $p->y + $Ry_BF * $p->z) + $Dx_BF;
            $p->y = $M_BF * ($Rz_BF * $p->x + $p->y - $Rx_BF * $p->z) + $Dy_BF;
            $p->z = $M_BF * (-$Ry_BF * $p->x + $Rx_BF * $p->y + $p->z) + $Dz_BF;
        }
    }
1315925
public function forward($p)
    {
        //alert("ll2m coords : ".coords);
        $lon = $p->x;
        $lat = $p->y;

        // convert to radians
        if ($lat * Common::R2D > 90.0 &&
                  $lat * Common::R2D < -90.0 &&
                  $lon * Common::R2D > 180.0 &&
                  $lon * Common::R2D < -180.0
        ) {
            Proj4php::reportError( "merc:forward: llInputOutOfRange: " . $lon . " : " . $lat );
            return null;
        }

        if (abs( abs( $lat ) - Common::HALF_PI ) <= Common::EPSLN) {
            Proj4php::reportError( "merc:forward: ll2mAtPoles" );
            return null;
        } else {
            if ($this->sphere) {
                $x = $this->x0 + $this->a * $this->k0 * Common::adjust_lon( $lon - $this->long0 );
                $y = $this->y0 + $this->a * $this->k0 * log( tan( Common::FORTPI + 0.5 * $lat ) );
            } else {
                $sinphi = sin( lat );
                $ts = Common::tsfnz( $this->e, $lat, $sinphi );
                $x = $this->x0 + $this->a * $this->k0 * Common::adjust_lon( $lon - $this->long0 );
                $y = $this->y0 - $this->a * $this->k0 * log( $ts );
            }

            $p->x = $x;
            $p->y = $y;

            return $p;
        }
    }
1050610
public function inverse($p)
    {
        $x = $p->x;
        $y = $p->y;

        $delta_x = $x - $this->x0;
        $delta_y = $y - $this->y0;

        // 1. Calculate z
        $z_re = $delta_y / $this->a;
        $z_im = $delta_x / $this->a;

        // 2a. Calculate theta - first approximation gives km accuracy
        $z_n_re = 1;
        $z_n_im = 0;  // z^0
        $z_n_re1;
        $z_n_im1;

        $th_re = 0;
        $th_im = 0;

        for ($n = 1; $n <= 6; $n++) {
            $z_n_re1 = $z_n_re * $z_re - $z_n_im * $z_im;
            $z_n_im1 = $z_n_im * $z_re + $z_n_re * $z_im;
            $z_n_re = $z_n_re1;
            $z_n_im = $z_n_im1;
            $th_re = $th_re + $this->C_re[$n] * $z_n_re - $this->C_im[$n] * $z_n_im;
            $th_im = $th_im + $this->C_im[$n] * $z_n_re + $this->C_re[$n] * $z_n_im;
        }

        // 2b. Iterate to refine the accuracy of the calculation
        //        0 iterations gives km accuracy
        //        1 iteration gives m accuracy -- good enough for most mapping applications
        //        2 iterations bives mm accuracy
        for ($i = 0; $i < $this->iterations; $i++) {
            $th_n_re = $th_re;
            $th_n_im = $th_im;
            $th_n_re1;
            $th_n_im1;

            $num_re = $z_re;
            $num_im = $z_im;

            for ($n = 2; $n <= 6; $n++) {
                $th_n_re1 = $th_n_re * $th_re - $th_n_im * $th_im;
                $th_n_im1 = $th_n_im * $th_re + $th_n_re * $th_im;
                $th_n_re = $th_n_re1;
                $th_n_im = $th_n_im1;
                $num_re = $num_re + ($n - 1) * ($this->B_re[$n] * $th_n_re - $this->B_im[$n] * $th_n_im);
                $num_im = $num_im + ($n - 1) * ($this->B_im[$n] * $th_n_re + $this->B_re[$n] * $th_n_im);
            }

            $th_n_re = 1;
            $th_n_im = 0;
            $den_re = $this->B_re[1];
            $den_im = $this->B_im[1];

            for ($n = 2; $n <= 6; $n++) {
                $th_n_re1 = $th_n_re * $th_re - $th_n_im * $th_im;
                $th_n_im1 = $th_n_im * $th_re + $th_n_re * $th_im;
                $th_n_re = $th_n_re1;
                $th_n_im = $th_n_im1;
                $den_re = $den_re + $n * ($this->B_re[$n] * $th_n_re - $this->B_im[$n] * $th_n_im);
                $den_im = $den_im + $n * ($this->B_im[$n] * $th_n_re + $this->B_re[$n] * $th_n_im);
            }

            // Complex division
            $den2 = $den_re * $den_re + $den_im * $den_im;
            $th_re = ($num_re * $den_re + $num_im * $den_im) / $den2;
            $th_im = ($num_im * $den_re - $num_re * $den_im) / $den2;
        }

        // 3. Calculate d_phi              ...                                    // and d_lambda
        $d_psi = $th_re;
        $d_lambda = $th_im;
        $d_psi_n = 1;  // d_psi^0

        $d_phi = 0;
        for ($n = 1; $n <= 9; $n++) {
            $d_psi_n = $d_psi_n * $d_psi;
            $d_phi = $d_phi + $this->D[$n] * $d_psi_n;
        }

        // 4. Calculate latitude and longitude
        // d_phi is calcuated in second of arc * 10^-5, so we need to scale back to radians. d_lambda is in radians.
        $lat = $this->lat0 + ($d_phi * Common::SEC_TO_RAD * 1E5);
        $lon = $this->long0 + $d_lambda;

        $p->x = $lon;
        $p->y = $lat;

        return $p;
    }
710369
function (p, p1, p2, sqDist) {
		var x = p1.x,
		    y = p1.y,
		    dx = p2.x - x,
		    dy = p2.y - y,
		    dot = dx * dx + dy * dy,
		    t;

		if (dot > 0) {
			t = ((p.x - x) * dx + (p.y - y) * dy) / dot;

			if (t > 1) {
				x = p2.x;
				y = p2.y;
			} else if (t > 0) {
				x += dx * t;
				y += dy * t;
			}
		}

		dx = p.x - x;
		dy = p.y - y;

		return sqDist ? dx * dx + dy * dy : new L.Point(x, y);
	}
851397
public function forward( $p ) {
        
        $lon = $p->x;
        $lat = $p->y;
        
        /* Forward equations
          ----------------- */
        $dlon = Common::adjust_lon( $lon - $this->long0 );
        $x = $this->x0 + $this->a * $dlon;
        $y = $this->y0 + $this->a * log( tan( (Common::PI / 4.0) + ($lat / 2.5) ) ) * 1.25;

        $p->x = $x;
        $p->y = $y;
        
        return $p;
    }
538639
public static float pointRectDistance (IRectangle r, IPoint p) {
        return FloatMath.sqrt(pointRectDistanceSq(r, p));
    }
400860
public function inverse($p)
    {
        /*
        $dlon;
        $xx;
        $yy;
        $xys;
        $c1;
        $c2;
        $c3;
        $al;
        $asq;
        $a1;
        $m1;
        $con;
        $th1;
        $d;
        */

        $p->x -= $this->x0;
        $p->y -= $this->y0;

        $con = Common::PI * $this->R;

        $xx = $p->x / $con;
        $yy = $p->y / $con;

        $xys = $xx * $xx + $yy * $yy;

        $c1 = -abs( $yy ) * (1.0 + $xys);
        $c2 = $c1 - 2.0 * $yy * $yy + $xx * $xx;
        $c3 = -2.0 * $c1 + 1.0 + 2.0 * $yy * $yy + $xys * $xys;

        $d = $yy * $yy / $c3 + (2.0 * $c2 * $c2 * $c2 / $c3 / $c3 / $c3 - 9.0 * $c1 * $c2 / $c3 / $c3) / 27.0;

        $a1 = ($c1 - $c2 * $c2 / 3.0 / $c3) / $c3;
        $m1 = 2.0 * sqrt( -$a1 / 3.0 );
        $con = ((3.0 * $d) / $a1) / $m1;

        if (abs($con) > 1.0) {
            if ($con >= 0.0) {
                $con = 1.0;
            } else {
                $con = -1.0;
            }
        }

        $th1 = acos($con) / 3.0;

        if ($p->$y >= 0) {
            $lat = (-$m1 * cos($th1 + Common::PI / 3.0) - $c2 / 3.0 / $c3) * Common::PI;
        } else {
            $lat = -(-m1 * cos($th1 + Common::PI / 3.0) - $c2 / 3.0 / $c3) * Common::PI;
        }

        if (abs($xx) < Common::EPSLN) {
            $lon = $this->$long0;
        }

        $lon = Common::adjust_lon(
            $this->long0 + Common::PI * ($xys - 1.0 + sqrt(1.0 + 2.0 * ($xx * $xx - $yy * $yy) + $xys * $xys)) / 2.0 / $xx
        );

        $p->x = $lon;
        $p->y = $lat;

        return $p;
    }
452160
public function inverse($p)
    {
        /*
        $rh;  // height above ellipsoid	
        $z;  // angle
        $sinz;
        $cosz; // sin of z and cos of z	
        $temp;
        $con;
        $lon;
        $lat;
        */

        $p->x -= $this->x0;
        $p->y -= $this->y0;
        $rh = sqrt($p->x * $p->x + $p->y * $p->y);

        if($rh > $this->a + 0.0000001) {
            Proj4php::reportError("orthoInvDataError");
        }

        $z = Common::asinz($rh / $this->a);

        $sinz = sin($z);
        $cosz = cos($z);

        $lon = $this->long0;

        if (abs($rh) <= Common::EPSLN) {
            $lat = $this->lat0;
        }

        $lat = Common::asinz($cosz * $this->sin_p14 + ($p->y * $sinz * $this->cos_p14) / $rh);
        $con = abs($this->lat0) - Common::HALF_PI;

        if (abs(con) <= Common::EPSLN) {
            if ($this->lat0 >= 0) {
                $lon = Common::adjust_lon($this->long0 + atan2($p->x, -$p->y));
            } else {
                $lon = Common::adjust_lon($this->long0 - atan2(-$p->x, $p->y));
            }
        }
        $con = $cosz - $this->sin_p14 * sin($lat);

        $p->x = $lon;
        $p->y = $lat;

        return $p;
    }
707919
public function forward($p)
    {
        $lon = $p->x;
        $lat = $p->y;

        $dlon = Common::adjust_lon( $lon - $this->long0 );
        $dlat = Common::adjust_lat( $lat - $this->lat0 );

        $p->x = $this->x0 + ($this->a * $dlon * $this->rc);
        $p->y = $this->y0 + ($this->a * $dlat );

        return $p;
    }
976826
public static GeoPosition getPosition(Point2D pixelCoordinate, int zoom, TileFactoryInfo info)
    {
        // p(" --bitmap to latlon : " + coord + " " + zoom);
        double wx = pixelCoordinate.getX();
        double wy = pixelCoordinate.getY();
        // this reverses getBitmapCoordinates
        double flon = (wx - info.getMapCenterInPixelsAtZoom(zoom).getX()) / info.getLongitudeDegreeWidthInPixels(zoom);
        double e1 = (wy - info.getMapCenterInPixelsAtZoom(zoom).getY())
                / (-1 * info.getLongitudeRadianWidthInPixels(zoom));
        double e2 = (2 * Math.atan(Math.exp(e1)) - Math.PI / 2) / (Math.PI / 180.0);
        double flat = e2;
        GeoPosition wc = new GeoPosition(flat, flon);
        return wc;
    }
1181362
private Point2d reflect(Point2d p, double x0, double y0, double x1, double y1) {

        double dx, dy, a, b;

        dx = (x1 - x0);
        dy = (y1 - y0);

        a = (dx * dx - dy * dy) / (dx * dx + dy * dy);
        b = 2 * dx * dy / (dx * dx + dy * dy);

        return new Point2d(a * (p.x - x0) + b * (p.y - y0) + x0, b * (p.x - x0) - a * (p.y - y0) + y0);
    }
580610
def for_point(cls, point, zoom):
        """"""
        latitude, longitude = point.latitude_longitude
        return cls.for_latitude_longitude(latitude=latitude, longitude=longitude, zoom=zoom)