Just for reference, here is the quick and dirty implementation I threw together after consulting several pages at Wikipedia:

    sub acos { atan2( sqrt(1-$_[0]*$_[0]), $_[0] ) }
    sub asin { atan2( $_[0], sqrt(1-$_[0]*$_[0]) ) }
    my $pi= atan2(0,-1);
    my @c= ( 33.943603, -118.408189, 39.17965, -76.668824 );
    my $lat1= $c[0]/180*$pi;
    my $lat2= $c[2]/180*$pi;
    my $dlong= ($c[1]-$c[3])/180*$pi;
    my $ang= acos( sin($lat1)*sin($lat2)
        + cos($lat1)*cos($lat2)*cos($dlong) );
    my $max= $ang * 3_963.19059;
    my $min= $ang * 3_949.90257;
    printf "%.1f .. %.1f miles (%.2f range)\n", $min, $max, $max-$min;
    my $h1= sin(($lat1-$lat2)/2);
    my $h2= sin($dlong/2);
    $ang= 2*asin( sqrt( $h1*$h1 + cos($lat1)*cos($lat2)*$h2*$h2 ) );
    $max= $ang * 3_963.19059;
    $min= $ang * 3_949.90257;
    printf "%.1f .. %.1f miles (%.2f range)\n", $min, $max, $max-$min;
__END__
2318.6 .. 2326.4 miles (7.80 range)
2318.6 .. 2326.4 miles (7.80 range)
[download]

And if somebody tells you that a >10% error is due to Earth being an oblate spheroid, tell them that such a discrepancy would be less than 0.34%.