LOCAL
     dir2, dir3, dir4 : direction :=
                 dummy_gri || direction([1.0, 0.0, 0.0]);
     val, mag         : REAL;
   END_LOCAL;

   IF (p1.dim <> 3) THEN
     RETURN(?);
   END_IF;
   REPEAT i := 1 TO 3;
     dir2.direction_ratios[i] := p2.coordinates[i] - p1.coordinates[i];
     dir3.direction_ratios[i] := p3.coordinates[i] - p1.coordinates[i];
     dir4.direction_ratios[i] := p4.coordinates[i] - p1.coordinates[i];
     mag := dir4.direction_ratios[i]*dir4.direction_ratios[i];
  END_REPEAT;
  mag := sqrt(mag);
  val := mag*dot_product(dir4, cross_product(dir2, dir3).orientation);
  RETURN(val);