#test of gravity corrections

sub free_air_normal ($){
  my $h = shift; #height of the free air correction

#h is station elevation - datum elevation
# note - for a station located above the datum
#        the free air corrected gravity should be
#        greater than the observed gravity

# This subroutine uses the normal linear model for
# free air corrections described in most textbooks

  my $free_air_correction = 0.3086*$h; #mgal

  return $free_air_correction;
};

sub free_air_precise ($$) {

	 my $h = shift; #height of the free air correction
   my $phi = shift; #geodetic latitude (degrees)
   my $pi = 3.14159265358979;

#h is station elevation - datum elevation
# note - for a station located above the datum
#        the free air corrected gravity should be
#        greater than the observed gravity

# This subroutine uses the normal linear model for
# free air corrections described in most textbooks
# and relies on the GRS80 ellipsoid as recommended by NAGDB
# working group 2003

  $phi = $phi*$pi/180;
  my $free_air_correction = (0.3087691-0.0014398*sin($phi)**2)*$h + 7.2125e-8*$h*$h; #mgal

  return $free_air_correction;

};

sub bullard_a ($$){

# This subroutine calculates the simple
# Bouguer gravity correction (flat slab)
# for a slab thickness of h and density
	my $h = shift; #meters
	my $density = shift; #kg/m3

	my $pi = 3.14159265358979;
	my $big_G = 6.673e-11;
	my $to_mgal = 1e5;

#h is station elevation - datum elevation
# note - for a station located above the datum
#        the Bouguer corrected gravity should be
#        less than the free air corrected gravity

	my $gba = 2*$pi*$big_G * $to_mgal * $density *$h;
	return $gba;


	};

sub bullard_b ($){

# This subroutine calculates the correction to the simple
# Bouguer gravity correction (flat slab)
# for a slab thickness of h and density to account
# for the curvature of the earth (spherical cap
# approximation to a radius of 166.7 km)

	my $h = shift; #meters

	#constants from Cogbill (1979)
	#constants are claculated for an
	#earth radius of 6371 km and density of 2670 kg /m3
  my $A = 1.464139e-3;
  my $B = 3.533047e-7;
  my $C = 1.002709e-13;
  my $D = 3.002407e-18;

#h is station elevation - datum elevation
# note - the Bullard B correction is subtracted from observations

	my $gbb = $A*$h - $B*$h**2 + $C*$h**3 + $D*$h**4;
	return $gbb;


	};

#######Main Program
$observed_grav = -800;
$phi = 45;
$ele = 4000;
$density = 2670;

	$free1 = free_air_normal ($ele);
	$free2 = free_air_precise($ele, $phi);
	$gba = bullard_a ($ele, $density);
  $gbb = bullard_b ($ele);
  $corrected_grav = $observed_grav + $free2 - $gba - $gbb;
print "elevation: $ele observed gravity: $observed_grav mgal\n";
print "normal free air correction: $free1 mgal\n";
print "precise free air correction: $free2 mgal\n";
print "bullard a correction: $gba mgal\n";
print "bullard b correction: $gbb  mgal\n";

print "corrected gravity: $corrected_grav  mgal\n";