intersect <- function(a, b, c, d) {
    #   test two line segments for intersection.
    #   modified from segseg in "Computational Geometry in C" 
    #   by Joseph O'Rourke

    p = c(0,0)

   denom = a[1] * ( d[2] - c[2] ) +
           b[1] * ( c[2] - d[2] ) +
           d[1] * ( b[2] - a[2] ) +
           c[1] * ( a[2] - b[2] );

   # If denom is zero, then segments are parallel: handle separately. 
   if (abs(denom) <= 1.0e-10) {
      #return  ParallelInt(a, b, c, d, p);
      #c(p,0)
      return (c(p,0)) 
    }

   num =    a[1] * ( d[2] - c[2] ) +
            c[1] * ( a[2] - d[2] ) +
            d[1] * ( c[2] - a[2] );
   if ( (num == 0.0) || (num == denom) ) code = 0;
   s = num / denom;

   num = -( a[1] * ( c[2] - b[2] ) +
            b[1] * ( a[2] - c[2] ) +
            c[1] * ( b[2] - a[2] ) );
   if ( (num == 0.0) || (num == denom) ) code = 0;
   t = num / denom;

   if      ( (0.0 < s) && (s < 1.0) &&
             (0.0 < t) && (t < 1.0) ) {
     code = 1;
     }
   else if ( (0.0 > s) || (s > 1.0) ||
             (0.0 > t) || (t > 1.0) ) {
     code = 0;
     }

   p[1] = a[1] + s * ( b[1] - a[1] );
   p[2] = a[2] + s * ( b[2] - a[2] );

   c(p,code);
}

#   Intersect two lines defined by a series of segments. Assume the
#   lines have common x-values and differ only in y
#   intersect is array of (x,y) values

intersect.lines <- function (y1,y2,x) {
    intersect = array(data = NA, dim = c(2,0), dimnames = NULL)
    for (i in 2:length(x)) {
       a = c(x[i-1], y1[i-1])
       b = c(x[i], y1[i])
       c = c(x[i-1], y2[i-1])
       d = c(x[i], y2[i])
       foo = intersect(a,b,c,d)
       if (foo[3]) {
            intersect = cbind(intersect, foo[1:2])
       }
    }
    intersect
}