gistfile1.php

<?php

class GEO
{
   const
    PI                = M_PI,
    TWOPI             = 6.28318530718,
    DE2RA             = 0.01745329252,
    RA2DE             = 57.2957795129,
    ERAD              = 6378.135,
    ERADM             = 6378135.0,
    AVG_ERAD          = 6371.0,
    FLATTENING        = 0.00335281066,
                                  // Earth flattening (WGS '84)
    EPS               = 0.000000000005,
    KM2MI             = 0.621371,
    GEOSTATIONARY_ALT = 35786.0;    // km
}

function intersectPoint($line1start, $line1end, $line2start, $line2end)
{
    $lat1A = $line1start['lat'];
    $lon1A = $line1start['lng'];
    $lat1B = $line1end['lat'];
    $lon1B = $line1end['lng'];

    $lat2A = $line2start['lat'];
    $lon2A = $line2start['lng'];
    $lat2B = $line2end['lat'];
    $lon2B = $line2end['lng'];

   $isInside = false;

   $v2 = $v1 = array();

   

   $d1 = distance($lat1A, $lon1A, $lat1B, $lon1B);
   $d2 = distance($lat2A, $lon2A, $lat2B, $lon2B);

   //
   // for path 1, setting up my 2 vectors, $v1 is vector
   // from center of the Earth to point A, $v2 is vector
   // from center of the Earth to point B.
   //
   $v1[0] = cos(deg2rad($lat1A)) * cos(deg2rad($lon1A));
   $v1[1] = sin(deg2rad($lat1A));
   $v1[2] = cos(deg2rad($lat1A)) * sin(deg2rad($lon1A));

   $v2[0] = cos(deg2rad($lat1B)) * cos(deg2rad($lon1B));
   $v2[1] = sin(deg2rad($lat1B));
   $v2[2] = cos(deg2rad($lat1B)) * sin(deg2rad($lon1B));

   //
   // $n1 is the normal to the plane formed by the three points:
   // center of the Earth, point 1A, and point 1B
   //
   $n1[0] = ($v1[1]*$v2[2]) - ($v1[2]*$v2[1]);
   $n1[1] = ($v1[2]*$v2[0]) - ($v1[0]*$v2[2]);
   $n1[2] = ($v1[0]*$v2[1]) - ($v1[1]*$v2[0]);

   //
   // for path 2, setting up my 2 vectors, $v1 is vector
   // from center of the Earth to point A, $v2 is vector
   // from center of the Earth to point B.
   //
   $v1[0] = cos(deg2rad($lat2A)) * cos(deg2rad($lon2A));
   $v1[1] = sin(deg2rad($lat2A));
   $v1[2] = cos(deg2rad($lat2A)) * sin(deg2rad($lon2A));

   $v2[0] = cos(deg2rad($lat2B)) * cos(deg2rad($lon2B));
   $v2[1] = sin(deg2rad($lat2B));
   $v2[2] = cos(deg2rad($lat2B)) * sin(deg2rad($lon2B));

   //
   // $n1 is the normal to the plane formed by the three points:
   // center of the Earth, point 2A, and point 2B
   //
   $n2[0] = ($v1[1]*$v2[2]) - ($v1[2]*$v2[1]);
   $n2[1] = ($v1[2]*$v2[0]) - ($v1[0]*$v2[2]);
   $n2[2] = ($v1[0]*$v2[1]) - ($v1[1]*$v2[0]);

   //
   // $v3 is perpendicular to both normal 1 and normal 2, so
   // it lies in both planes, so it must be the line of
   // intersection of the planes. The question is: does it
   // go towards the correct intersection point or away from
   // it.
   //
   $v3a[0] = ($n1[1]*$n2[2]) - ($n1[2]*$n2[1]);
   $v3a[1] = ($n1[2]*$n2[0]) - ($n1[0]*$n2[2]);
   $v3a[2] = ($n1[0]*$n2[1]) - ($n1[1]*$n2[0]);

   //
   // want to make $v3a a unit vector, so first have to get
   // magnitude, then each component by magnitude
   //
   $m = sqrt($v3a[0]*$v3a[0] + $v3a[1]*$v3a[1] + $v3a[2]*$v3a[2]);
   $v3a[0] /= $m;
   $v3a[1] /= $m;
   $v3a[2] /= $m;

   //
   // calcu$lating intersection points 3A & 3B. A & B are
   // arbitrary designations right now, $later we make A
   // the one close to, or within, the paths.
   //
   $lat3A = asin($v3a[1]);
   if (($lat3A > GEO::EPS) || (-$lat3A > GEO::EPS))
      $lon3A = asin($v3a[2]/cos($lat3A));
   else
      $lon3A = 0.0;

   $v3b[0] = ($n2[1]*$n1[2]) - ($n2[2]*$n1[1]);
   $v3b[1] = ($n2[2]*$n1[0]) - ($n2[0]*$n1[2]);
   $v3b[2] = ($n2[0]*$n1[1]) - ($n2[1]*$n1[0]);

   $m = sqrt($v3b[0]*$v3b[0] + $v3b[1]*$v3b[1] + $v3b[2]*$v3b[2]);
   $v3b[0] /= $m;
   $v3b[1] /= $m;
   $v3b[2] /= $m;

   $lat3B = asin($v3b[1]);
   if (($lat3B > GEO::EPS) || (-$lat3B > GEO::EPS))
      $lon3B = asin($v3b[2]/cos($lat3B));
   else
      $lon3B = 0.0;

   //
   // get the distances from the path endpoints to the two
   // intersection points. these values will be used to determine
   // which intersection point lies on the paths, or which one
   // is closer.
   //
   $d1a3a = distance($lat1A, $lon1A, $lat3A, $lon3A);
   $d1b3a = distance($lat1B, $lon1B, $lat3A, $lon3A);
   $d2a3a = distance($lat2A, $lon2A, $lat3A, $lon3A);
   $d2b3a = distance($lat2B, $lon2B, $lat3A, $lon3A);

   $d1a3b = distance($lat1A, $lon1A, $lat3B, $lon3B);
   $d1b3b = distance($lat1B, $lon1B, $lat3B, $lon3B);
   $d2a3b = distance($lat2A, $lon2A, $lat3B, $lon3B);
   $d2b3b = distance($lat2B, $lon2B, $lat3B, $lon3B);

   if (($d1a3a < $d1) && ($d1b3a < $d1) && ($d2a3a < $d2)
                    && ($d2b3a < $d2))
   {
      $isInside = true;
   }
   else if (($d1a3b < $d1) && ($d1b3b < $d1) && ($d2a3b < $d2)
                         && ($d2b3b < $d2))
   {
      //
      // 3b is inside the two paths, so swap 3a & 3b
      //
      $isInside = true;

      $m = $lat3A;
      $lat3A = $lat3B;
      $lat3B = $m;
      $m = $lon3A;
      $lon3A = $lon3B;
      $lon3B = $m;
   }
   else
   {
      //
      // figure out which one is closer to the path
      //
      $d1 = $d1a3a + $d1b3a + $d2a3a + $d2b3a;
      $d2 = $d1a3b + $d1b3b + $d2a3b + $d2b3b;

      if ($d1 > $d2)
      {
         //
         // Okay, we are here because 3b {$lat3B,$lon3B} is closer to
         // the paths, so we need to swap 3a & 3b. the other case is
         // already the way 3a & 3b are organized, so no need to swap
         //
         $m = $lat3A;
         $lat3A = $lat3B;
         $lat3B = $m;
         $m = $lon3A;
         $lon3A = $lon3B;
         $lon3B = $m;
      }
   }

   //
   // convert the intersection points to degrees
   //
   $lat3A *= GEO::RA2DE;
   $lon3A *= GEO::RA2DE;
   $lat3B *= GEO::RA2DE;
   $lon3B *= GEO::RA2DE;

   if($lat3A && $lon3A)
   {
       return array('lat' => $lat3A, 'lng' => $lon3A, 'e' => array($lat3B, $lon3B), 'i' => ($isInside ? 'true' : 'false'));
   }else{
       return null;
   }
   
}