#calculate magnetic anomaly for
#polygonal prism
#C. B. Connor

use warnings;
use Math::Trig;
$pi = 3.14159;

#experimental 4-sided prism
#x must be the northing (north coordinate), 
#y must be the easting (east coordinate),
#coordinates must be entered in clockwise order

#map units can be arbitary by must be consistent
#for example, all units in meters
$x[0] = 15; #north
$y[0] = 10; #east
$x[1] = 5;
$y[1] = 15;
$x[2] = -10;
$y[2] = -10;
$x[3] = -5;
$y[3] = -20;

#depth to top
$z1 = 2;
#depth to bottom
$z2 = 4;

$sides = 4;

#set magnetization (amp/m)
$minc = 45*$pi/180; #inclination down in rad
$mdec = 0*$pi/180;  #declination east in rad
$mi = 1;

#calculate direction cosines for magnetiation
$ml = cos($minc) * cos($mdec);
$mm = cos($minc) * sin($mdec);
$mn = sin($minc);

#components of magnetization in x,y,z directions
$mx = $mi * $ml;
$my = $mi * $mm;
$mz = $mi * $mn;

#set earths field
$einc = 45*$pi/180;  #inclination down in rad
$edec = 0*$pi/180;   #declination east in rad

$el = cos($einc) * cos($edec);
$em = cos($einc) * sin($edec);
$en = sin($einc); 

#proportionality constant is
# magnetic permeability * 1e9 
#to convert of nT
$prop = 400*$pi;

for ($j =0; $j<10; $j++) {
for ($k =0; $k<10; $k++) {
#px is the northing of the observation point
#py is the easting of the observation point
$px = -20 + $j*5;
$py = -20 + $k*5;

#set the volume integrals to zero
#for each observation point
$v1 = 0;
$v2 = 0;
$v3 = 0;
$v4 = 0;
$v5 = 0;
$v6 = 0;

#calculate volume integrals
for ($i=0; $i<$sides; $i++) {

    #if last side of polygon
    if($i == $sides-1) {
       $x1 = $x[$i]-$px;
       $y1 = $y[$i]-$py;
       $x2 = $x[0]-$px;
       $y2 = $y[0]-$py;
    } else {
       #all other polygon sides
       $x1 = $x[$i]-$px;
       $y1 = $y[$i]-$py;
       $x2 = $x[$i+1]-$px;
       $y2 = $y[$i+1]-$py;
    }

    #calculate geometry
    $delta_x = $x2-$x1;
    $delta_y = $y2-$y1;
    $delta_s = sqrt($delta_x**2 + $delta_y**2);

    #avoid division by zero if delta_s = 0
   # if ($delta_s == 0) {$delta_s = 0.1};
    $c = $delta_y/$delta_s;
    $s = $delta_x/$delta_s;
    $p = ($x1*$y2 - $x2*$y1)/$delta_s;

    $d1 = $x1*$s + $y1*$c;
    $d2 = $x2*$s + $y2*$c;

    $r1sq = $x1**2 + $y1**2;
    $r2sq = $x2**2 + $y2**2;
    $r11 = sqrt($r1sq + $z1**2);
    $r12 = sqrt($r1sq + $z2**2);
    $r21 = sqrt($r2sq + $z1**2);
    $r22 = sqrt($r2sq + $z2**2);

    $f = log(($r22+$z2)/($r12+$z2) * ($r11+$z1)/($r21+$z1));
    $q = log(($r22+$d2)/($r12+$d1) * ($r11+$d1)/($r21+$d2));
    $w = atan2(($z2*$d2),($p*$r22)) - atan2(($z2*$d1),($p*$r12)) - atan2(($z1*$d2),($p*$r21)) + atan2(($z1*$d1),($p*$r11));
    
    $v1 += $s*$c*$f - $c*$c*$w;
    $v2 += $s*$c*$w + $c*$c*$f;
    $v3 += $c*$q;
    $v4 += -($s*$c*$f + $s*$s*$w);
    $v5 += -($s*$q);
    $v6 += $w;
}

#calculate the components of the magnetic field
$bx = $prop*($mx*$v1 + $my*$v2 + $mz* $v3);
$by = $prop*($mx*$v2 + $my*$v4 + $mz* $v5);
$bz = $prop*($mx*$v3 + $my*$v5 + $mz* $v6);

#calculate total anomaly
#this calculation works for anomaly << total field strength
$b_total = $el*$bx + $em*$by + $en*$bz;
$bint =  int($b_total);
#uncomment these lines to check volume integral
# v4 = -($v1+$v6)
#$test = -($v1+$v6);
#print"$v4 $test\n";

#print results 
#note the format of the output is easting, northing, mag
print " $py $px $bint\n";
}

}