Gravity anomaly due to a 2D polygon

// Calculate the gravity anomaly due to a 2d polygon in cross section
// Chuck Connor
// 9/2018
// using the method of Manik Talwani and colleagues:
// Talwani et al (1959) Rapid gravity calculatons of two dimensional bodies...JGR 64:49-59. 

// USAGE:
// Add these two script lines to the html header section:
// <script src="http://cdnjs.cloudflare.com/ajax/libs/p5.js/0.5.6/p5.js"><script>
// Here, the script source is stored with the html page
// <script src="talwani.js"><script>
//
// Add this div line to the html body section:
// <div id="sketch-holder" style="position: relative;"><div>

// The following is the script source code
var density_contrast = 500.0;
var x_zero = 50;
var  x_scale = 100; //1 km = 50 screen units
var depth_zero = 200; //where the depth axis starts in screen units
var depth_scale = 100; //keep vert exaggeration 1:1
var  mgal_zero = 170;  //less than depth_zero 
var mgal_scale = 10;

//create an array of corner objects
var corner = [];

//location of corner object
//these are global so mouse state can be checked
var x = [];
var y = [];

//number of corner objects
var num_corners;

//size of corner box
var boxSize = 5;

//keep track of mouse state for each corner object
//these are global so mouse state can be checked
var locked = [];
var overBox = [];

// i indexes corner objects
// this are global so mouse state can be checked
var i;
var densitySlider;

let myFont;
function preload() {
  myFont = loadFont('talwani/p5/DejaVuSans.ttf');
}

function setup() {
  var canvas = createCanvas(600, 480);
  canvas.parent('sketch-holder');
  
  stroke(0);
  fill(0)
   .strokeWeight(.30);
  textFont(myFont, 16);
  
  densitySlider = createSlider(1,255,0);
  densitySlider.parent('sketch-holder');
  densitySlider.position(450,20);
  densitySlider.style('width', '80px');
  densitySlider.parent('sketch-holder'); 

  //noStroke();
  rectMode(RADIUS);

  //specify number and intial location of polygon corners
  num_corners = 4;
  x[0] = 300;
  y[0] = 200;
  x[1] = 500;
  y[1] = 250;
  x[2] = 300;
  y[2] = 400;
  x[3] = 100;
  y[3] = 250;

 //create corner objects and initialize mouse state for each corner
  for (i=0; i x[i] - boxSize && mouseX < x[i] + boxSize &&
      mouseY > y[i] - boxSize && mouseY < y[i] + boxSize) {
      overBox[i] = true;        
      if (!locked[i]) {
        fill(255);           
      }
    } else {
      overBox[i] = false;
    }
  }

  this.display = function () {
    // Draw the corner box and connecting line
  if (locked[i] == true) { 
     boxSize = 10;
  }else{
     boxSize = 8;
  };
  fill('red');
    rect(x[i], y[i], boxSize, boxSize);
    if (i == num_corners-1) {
       line(x[num_corners-1], y[num_corners-1], x[0], y[0]);
    }else{
       line(x[i], y[i], x[i+1], y[i+1]);
    }
  }
}

function mousePressed() {
  //check the state f each corner, highlight slected corner
  for (i=0; i width) {x[i]=width};
      if (y[i] > height) {y[i]=height};
    }
  }
}

function mouseReleased() {
  // no corner selected on mouse release
  for (i = 0; i < num_corners; i++) {
    locked[i] = false;
    fill(255); 
  }
}

function talwani (x1,x2,z1,z2,density){
  // assume that integration is clockwise
  var G=6.67e-11;
  var pi = 3.14159265358979;
  var denom, alpha, beta, factor;
  var r1sq, r2sq, term1, term2, zz, grav ;

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

  alpha = (x2-x1)/denom;
  beta = (x1*z2-x2*z1)/denom;
  factor = beta/(1+alpha*alpha);
  r1sq = (x1*x1 + z1 *z1);
  r2sq = (x2*x2 + z2 *z2);
  term1 = 0.5*(log(r2sq) - log(r1sq));
  term2 = atan2(z2,x2)-atan2(z1,x1);
  zz = factor*(term1-alpha*term2);
  grav = 2*G*density*zz*1e5;
  return grav; // return gravity in mgal
};

# USAGE:
# type: perl talwani.pl
# output file created by the code is named,
# output.dat

#######################################
# Function which
# calculates the gravity
# anomaly associated with one side of a 2D
# polygon (viewed in cross section and infinite
# in and out of plane of screen).
# INPUTS:
#  x1, x2, z1, z2, density;
#  in this order!
# OUTPUT:
#  gravity in mGals
########################################
sub talwani ($$$$$){

  my $x1= shift; #meters
  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 an output file
open OUTPUT, ">  outfile.dat " or die "Can't open output file 1: $!";

# create the polygon
# vertices must be 
# in clockwise order!
$x[0] = 0; 
$z[0] = 50;
$x[1] = 275; 
$z[1] = 50;
$x[2] = 90; 
$z[2] = 200;
$density_contrast = 1000;
$n = 3; # 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 OUTPUT "$offset $gravity\n";
  $gravity = 0;
};