The Buried Sphere

#!/usr/bin/perl
# USAGE: perl buried_sphere.pl > buried_sphere.dat

use warnings;
use strict;
use Math::SigFigs;

# Define some example constants
my $sphere_radius = 4000; #m
my $sphere_depth = 5000; #m (sphere_depth > sphere_radius)
my $sphere_density_contrast = 500; #kg/m3

for (my $i=0; $i<=10; $i++) {
  my $offset = $i * 1000; #horizontal distance (m)
  my $gravity = gsphere($sphere_radius, $sphere_depth, $offset, $sphere_density_contrast);

  print "$offset ";
  print FormatSigFigs($gravity,3), "\n";
}

###################################################################
#subroutine to calculate the gravity anomaly due to a buried sphere

# inputs:
#  sphere radius (m), 
#  sphere depth (m), 
#  calculation point: 
#      horizontal distance from center of sphere to measurement(m)
#  density constrast (kg/m3)
# output:
#   gravity anomaly (mGal) 
################################################################### 
sub gsphere ($$$$) {

  my $pi = 3.14159265358979;
  my $a = shift;           # radius of the sphere (m)
  my $z = shift;           # depth of the center of the sphere (m)
  my $x = shift;           # horizontal distance of the 
                           # location of the gravity observation (m)
                         
  my $del_density = shift; # density cotrast (kg/m3)
  my $big_G = 6.67e-11;    #gravitational const.
  my $to_mgal = 1e5;       #convert SI to mGal

  my $V = (4/3)* $pi * $a ** 3;
  my $grav = $to_mgal * $big_G * $del_density * $V * $z * (($z*$z + $x*$x) ** (-3/2));

  return $grav;
};

// 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 with the html page
// <script src="canvas_gravity_sphere.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 x = 100;
var y = 100;
var radius =  50; //in screen units, 1 km = 100 screen units
var x_zero = 50;
var x_scale = 100; //1 km = 100 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 = 20;
var density_contrast = 500.0;

let myFont;
function preload() {
	// Font file is DejaVuSans.ttf,
	// copied to web page directory
  myFont = loadFont('DejaVuSans.ttf');
}
     
function setup() {
  var canvas = createCanvas(600,550);
 
  // Move the canvas so it's inside our div tag:
  // div id="sketch-holder"
  canvas.parent('sketch-holder');

  fill(0)
   .strokeWeight(.25);
  textFont(myFont, 16);
}

function draw() {
  background('#006747');
 
 var gravity;
 var x1km, y1km, r;
 var del_grav;
 var offset;
 var x1,  z1, m;
  // lerp() calculates a number between two numbers at a specific increment. 
  // The amt parameter is the amount to interpolate between the two values 
  // where 0.0 equal to the first point, 0.1 is very near the first point, 0.5 
  // is half-way in between, etc.  
  
  // Here we are moving 5% of the way to the mouse location each frame
  x = lerp(x,mouseX,0.05);
  y = lerp(y,mouseY,0.05);
  
  if (x < x_zero+radius) {x = x_zero + radius;}
  if (x > width - radius) {x = width - radius;} 
  if (y < depth_zero+radius) {y = depth_zero + radius;}
  
  ellipse(x, y, 2*radius, 2*radius);
  
  // calculate and draw gravity
  // count in screen units
  for (offset = x_zero; offset < width; offset += 10) { 
    
     x1 = (x) - offset;
     z1 = (y); 

     x1km = x1 / x_scale * 1000.0;
     y1km = (z1 - depth_zero) / depth_scale * 1000.0;
     r = radius / x_scale * 1000.0;
  
     gravity =gsphere(r, x1km, y1km, density_contrast);

    // convert gravity to integer
    del_grav = gravity * mgal_scale;
     m = mgal_zero - del_grav;
    // plot point
    ellipse(offset, m, 6); 
  }

  //draw axes and labels
  drawAxes(x_zero, x_scale, depth_zero, depth_scale, mgal_zero, mgal_scale);
}

function drawAxes(x_zero, x_scale, depth_zero, depth_scale, mgal_zero, mgal_scale) {
 fill('#CFC493');
 textAlign(RIGHT);
 for (var i = 0; i < 8; i++) {
  var d = i * depth_scale + depth_zero;
  line(x_zero, d, x_zero - 10, d);
  text(i, x_zero - 12, d + 5);
 }
 
 //draw and label mgal axis
 textAlign(RIGHT);
 for (var i = 0; i < 8 ; i++) {
  var d= mgal_zero - 2 * i * mgal_scale;
  line(x_zero, d, x_zero - 10, d);
  text(i * 2, x_zero - 12, d + 5);
 }
 
  //draw label the horizontal axis
  textAlign(CENTER);
  for (var i = 0; i < 12; i++) {
    var x = i * x_scale + x_zero;
    line(x, depth_zero - 35, x, depth_zero - 25);
    line(x, depth_zero-5, x, depth_zero + 5);
    var x_value =  i;
    text(x_value, x, depth_zero - 10);
  }
 
  // rotate(3.14159/2.0);
   textAlign(CENTER);
   text("z (km)", 100, 350);
   text("x (km)", 500, depth_zero+25);  
   text("g (mgal)", 100, 100);   
   stroke(255);
   fill(255);
   //depth arrow
   line(75, 360, 75, 400);
   triangle(70, 400, 80, 400, 75, 425);
    //horiontal arrow
   line(520, depth_zero+20, 555, depth_zero+20);
   triangle(555, depth_zero + 15, 555, depth_zero + 25, 580, depth_zero + 20);
    //mgal arrow
   line(75,80, 75, 40);
   triangle(70,40,80,40,75,15);
}

/************************* 
calculate the gravity
  anomaly due to a sphere

  g-sphere = 4/3 rho G pi a^3 (1/(x^2+^2) (z)/(sqrt(x^2+z^2))
  distance units in meters
  returns mgal

  a = radius of sphere;
  x1 is horizontal offset
  z1 is depth; 
***************************/
function gsphere(a, x2, z2, del_density){

  var g;
  var pi;
  var grav;
  
  //define consts.
  g = 6.67e-11;
  pi = 3.14159265358979;
  grav = 4.0/3.0 * pi * g * del_density * 1.0e5 * a * a * a * 1/(x2 * x2 + z2 * z2) * z2/(sqrt(x2 * x2 + z2 * z2));

  return grav; //return gravity in mgal
}