sub gbox ($$$$$$$$$$) {

my $x0 = shift; #north location of observation pt
my $y0 = shift; #east location of observation point
my $z0 = shift; #elevation of observation point (positive down)
my $x1 = shift; #northing of south side of prism
my $y1 = shift; #easting of west side of prism
my $z1 = shift; #depth to top of prism
my $x2 = shift; #northing of north side of prism
my $y2 = shift; #easting of east side of prism
my $z2 = shift; #depth to base of prism
my $density = shift; #density contrast represented by prism

#distances are in meters
#density is in kg/m3

my @isign = (-1,1);
my $big_G = 6.67e-11;
my $twopi = 6.2831853;
my $to_mgal = 1e5;

$x[0] = $x0-$x1;
$y[0] = $y0-$y1;
$z[0] = $z0-$z1;
$x[1] = $x0-$x2;
$y[1] = $y0-$y2;
$z[1] = $z0-$z2;

$sum = 0;

for ($i=0; $i<2; $i++) {
  for ($j=0; $j<2; $j++) {
    for ($k=0; $k<2; $k++) {

      $rijk=sqrt($x[$i]*$x[$i]+$y[$j]*$y[$j]+$z[$k]*$z[$k]);
      $ijk= $isign[$i]*$isign[$j]*$isign[$k];

      $arg1 = atan2($x[$i]*$y[$j],$z[$k]*$rijk);
      if ($arg1 < 0){$arg1+=$twopi};

      $arg2=$rijk+$y[$j];
      $arg3=$rijk+$x[$i];

      if ($arg2<=0){$arg2=1e-6};
      if ($arg3<=0){$arg3=1e-6};

      $arg2 = log($arg2);
      $arg3 = log($arg3);

      $sum += $ijk*($z[$k]*$arg1-$x[$i]*$arg2 - $y[$j]*$arg3);
    };
  };
};

$gravity = $density*$big_G*$to_mgal*$sum;
return $gravity;

};


####### Main Program #############
#

open OUTPUT, ">  gbox.dat " or die "Can't open output file 1: $!";

# define the extent of the buried prism
$n1 = -200;
$n2 = 200;
$e1 = -200;
$e2 = 200;
$z1 = 100;
$z2 = 200;

$obsz = 0;
$den_contrast = 1000;

for ($obsn = -1000; $obsn < 1010; $obsn += 25){
	for ($obse = -1000; $obse < 1010; $obse += 25){
     my $grav = gbox($obsn,$obse,$obsz,$n1,$e1,$z1,$n2,$e2,$z2, $den_contrast);
     print OUTPUT "$obse $obsn $grav\n";
  };
};

system "xyz2grd gbox.dat -Ggbox.grd -I25 -R-1000/1000/-1000/1000";
system "makecpt -Chaxby -T-1/4/0.01 > gbox.cpt";
system "grdimage gbox.grd -X5i -JX4i -R -Cgbox.cpt -E200 -K > gbox.ps";
system "grdcontour gbox.grd -JX4i -R -C0.1 -O -B100a200:\"Easting (m)\":/100a200:\"Northing (m)\":WSne -O -K >> gbox.ps";
system "psxy prism.dat -: -JX4i -R -W2p  -O >> gbox.ps";