use Math::MatrixReal;
#linear inversion for density of a buried sphere

sub gsphere ($$$) {

#subroutine to calculate gravity anomaly due to buried sphere
my $pi = 3.14159265358979;
my $r = shift; #input sphere radius in km
my $z = shift; #input sphere depth in km
my $x = shift; #input offset (location of gravity obs) in km

my $big_G = 6.67e-11;
my $to_mgal = 1e5;
my $R_m = 1000 * $r;
my $z_m = 1000 * $z;
my $x_m = 1000 * $x;

my $V = (4/3)* $pi * $R_m ** 3;

$grav = $to_mgal * $big_G * $V * $z_m* ((($z_m*$z_m + $x_m * $x_m))**(-3/2));

return $grav;
};

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

open OUTPUT, ">  outfile.dat " or die "Can't open output file 1: $!";
open OUTPUT2, ">  mod_file.dat " or die "Can't open output file 2: $!";

my $r = 1.0;
my $z = 2.0;
my $i=0;

for ($offset=-15.0; $offset<15.0; $offset+=0.1)
{

$x[$i] = $offset;
$gravity[$i] = 0.2 - 0.1*$offset - 0.003*$offset*$offset + 200*gsphere($r, $z, $offset-2)+ 0.5*rand();

print OUTPUT "$x[$i] $gravity[$i] \n";
$i++
};


system "psxy outfile.dat -JX6i -R-15/15/-3/3 -SC0.01i -G0 -W1p -B2a4:\"offset (km)\":/0.1a0.5:\"gravity anomaly (mgal)\":WSne -K > outfile.ps";

$n=$i;
$i=0;

#build the A mtx
$A = new Math::MatrixReal($n,4);
$grav_data = new Math::MatrixReal($n,1);

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

	$A->assign($i,1,1);
	$A->assign($i,2,$x[$i]);
	$A->assign($i,3,$x[$i]*$x[$i]);
	$A->assign($i,4,gsphere($r, $z, $x[$i]-2));

	$grav_data->assign($i,1,$gravity[$i]);
	}
$AT = $A->new(4,$n);

$AT->transpose($A);
$AAT = $AT * $A;
$inv_AAT = $AAT->inverse;
$AAAT = $inv_AAT * $AT;
$p = $AAAT*$grav_data;
print "p=\n$p";

$c1 = $p->element(1,1);
$c2 = $p->element(2,1);
$c3 = $p->element(3,1);
$c4 = $p->element(4,1);

for ($offset=-15.0; $offset<15.0; $offset+=0.1)
{

$x[$i] = $offset;
$mod_gravity[$i] = $c1 + $c2*$offset + $c3*$offset*$offset + $c4*gsphere($r, $z, $offset-2);

print OUTPUT2 "$x[$i] $mod_gravity[$i] \n";
$i++
};

system "psxy mod_file.dat -JX6i -R-15/15/-3/3 -O -W2p >> outfile.ps";
system "rm mod_file.dat";
system "rm outfile.dat";