sub talwani ($$$$$){

#calculate the gravity
#anomaly associated with one side of a 2D
#polygon (viewed in cross section and infinite
#in and out of plane fof screen)

#pass subroutine: x1, x2, z1, z2, density;
# in that order
my $x1= shift;
my $x2= shift;
my $z1= shift ;
my $z2= shift;
my $del_density = shift;

# assume that integration is clockwise

my $G=6.67e-11;
my $pi = 3.14159265358979;

#avoid singularites
if ($x1 == 0) {$x1=0.0001};
if ($x2 == 0) {$x2=0.0001};
if (($x2-$x1)==0) {$x2 = $x1 - 0.0001};
my $denom = $z2-$z1;
if ($denom == 0){$denom = 1e-6};


my $alpha = ($x2-$x1)/$denom;
my $beta = ($x1*$z2-$x2*$z1)/$denom;
my $factor = $beta/(1+$alpha*$alpha);

my $r1sq = ($x1*$x1 + $z1 *$z1);
my $r2sq = ($x2*$x2 + $z2 *$z2);

my $term1 = 0.5*(log($r2sq) - log($r1sq));
my $term2 = atan2($z2,$x2)-atan2($z1,$x1);

my $zz = $factor*($term1-$alpha*$term2);

$grav = 2*$G*$del_density*$zz*1e5;

return $grav; #return gravity in mgal
};


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

open OUTPUT, ">  outfile.dat " or die "Can't open output file 1: $!";
#create the polygon
#must be clockwise
$x[0]= 100.1; $z[0] = 50.001;
$x[1]= 200.001; $z[1] = 100.001;
$x[2]= 90.001; $z[2] = 200.001;
$x[3]= 50.001; $z[3] = 110.001;
$density_contrast = 1000;

$n=4; #number of polygon sides

for ($offset=-300; $offset<500; $offset+=25){

for ($i=0; $i<$n; $i++) {

if($i==$n-1) {
	$i2=0
  } else {
  $i2=$i+1
  };

	$x1 = $x[$i]-$offset;
  $x2 = $x[$i2]-$offset;
  $z1 = $z[$i];
  $z2 = $z[$i2];


  $gravity += talwani($x1, $x2, $z1, $z2, $density_contrast);
#print "i=$i, i2 = $i2, gravity = $gravity\n";
};

print OUTPUT "$offset $gravity\n";
$gravity =0;

	};

system "psxy outfile.dat -JX6i/3i -Y3i -R-300/500/-5/5 -W1p -B25a50/0.1a0.5:\"gravity anomaly (mgal)\":Wsne   -K > gpoly.ps";
system "psxy gpoly_model.dat -JX6i/3i -Y-3i -R-300/500/-300/0 -W1p -B25a50:\"offset (m)\":/25a50:\"depth (m)\":WSne   -O >> gpoly.ps";
system "rm outfile.dat";