<?php
    
// +----------------------------------------------------------------------+
    // | PHP version 4                                                        |
    // +----------------------------------------------------------------------+
    // | Copyright (c) 2003-2004 Michal Migurski                              |
    // +----------------------------------------------------------------------+
    // | This source file is subject to version 3.0 of the PHP license,       |
    // | that is bundled with this package in the file LICENSE, and is        |
    // | available through the world-wide-web at the following url:           |
    // | http://www.php.net/license/3_0.txt.                                  |
    // | If you did not receive a copy of the PHP license and are unable to   |
    // | obtain it through the world-wide-web, please send a note to          |
    // | license@php.net so we can mail you a copy immediately.               |
    // +----------------------------------------------------------------------+
    // | Author: Michal Migurski <mike@teczno.com>                            |
    // +----------------------------------------------------------------------+
    //
    // $Id: Lambert_Azimuthal_EqualArea.php,v 1.3 2004/11/24 19:18:44 migurski Exp $
    /* vim: set expandtab tabstop=4 shiftwidth=4 softtabstop=4: */

    // theta: latitude (N/S)
    // lambda: longitude (E/W)
    // DON'T FORGET: deg2rad()
    
    
require_once('Linear.php');
    
   
/** Lambert_Azimuthal_EqualArea_Projection
    * calculates Lambert Azimuthal Equal-Area Projection
    * http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html
    */
    
class Lambert_Azimuthal_EqualArea_Projection extends Linear_Projection
    
{
        var 
$transform;
        var 
$inverse;
    
        var 
$lat0;
        var 
$lon0;
        
        function 
Lambert_Azimuthal_EqualArea_Projection($lat0=0$lon0=0)
        {
            
$this->initializeTransform();

            
$this->lat0 $lat0;
            
$this->lon0 $lon0;
        }
        
       
/** _GPSToMap
        * calculates map x,y coordinate from gps lat/long coordinates
        *
        * @param    lat     float       latitude, in radians
        * @param    lon     float       longitude, in radians
        *
        * @return   array   'x', 'y' elements, map coordinates
        */
        
function _GPSToMap($lat$lon)
        {
            
$k sqrt(/ (+ (sin($this->lat0) * sin($lat)) + (cos($this->lat0) * cos($lat) * cos($lon $this->lon0))));
        
            
$x $k cos($lat) * sin($lon $this->lon0);
            
$y $k * ((cos($this->lat0) * sin($lat)) - (sin($this->lat0) * cos($lat) * cos($lon $this->lon0)));
    
            
$vector = new Math_Matrix(array(array($x$y1)));
            
$vector->multiply($this->transform);

            return array(
'x' => $vector->getElement(00),
                         
'y' => $vector->getElement(01));
        }
    
       
/** _mapToGPS
        * calculates gps lat/long coordinate from map coordinates.
        *
        * @param    x       float       horizontal map position
        * @param    y       float       vertical map position
        *
        * @return   array   'lat', 'lon' elements, gps coordinates in radians
        */
        
function _mapToGPS($x$y)
        {
            
$vector = new Math_Matrix(array(array($x$y1)));
            
$vector->multiply($this->inverse);
            
            
$x $vector->getElement(00);
            
$y $vector->getElement(01);

            
$p sqrt(pow($x2) + pow($y2));
            
$c asin(0.5 $p);
        
            if (
$p == 0) {
                
$lat $this->lat0;
                
$lon $this->lon0;
            } else {
                
$lat asin((cos($c) * sin($this->lat0)) + (($y sin($c) * cos($this->lat0)) / $p));
                
$lon $this->lon0 atan(($x sin($c)) / (($p cos($this->lat0) * cos($c)) - ($y sin($this->lat0) * sin($c))));
            }
    
            return array(
'lat' => $lat'lon' => $lon);
        }

    }

?>