Gravity simulator under Earth's surface

<?php
/*
    Method:
1. Construct spherical lattice that will be contained within
another larger sphere that we will be testing the gravity at -
at different points.

2. Move the point from the centre of the large, containing sphere
to an outside location of any choosing.

3. Display the results at the test points on a graph.

4. Win.

    Method for constructing spherical lattice:
1. inputs, radius of the larger sphere (a: large whole number, preferably a multiple of 10), radius of the interior
spheres (b: small whole number that should evenly and naturally divide the radius of the larger sphere in whole number ratios)

2. A valid internal sphere is one whose surface (radius) does not exceed the boundaries of the larger containing sphere, and does not over lap with any other spheres. For the sake of simplicity (not accuracy, this is an estimate), I am going to treat the interior spheres as cubes (mainly because I am going to be constructing a cubic lattice), this is because using pythagoras' theorem to test whether the corners of the cubes exceed the distance from the center of the sphere to the surface of the larger containing sphere. (pythagoras' theorem in 3 dimensions r^2 = (x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2

3. Build the sphere using a top down view (treating the XY plane as the top view), start at position (b, b, sqrt(a^2-(2*b)^2))
then copy this location into the other 3 quadrants. These points will be loaded into an array and given unique identifiers.

4. Move to the next level down, (b, b, [start y]-b*2*1), the times by 1 indicates where the level itteration number goes, this will cause the level to decrement down (I know it's a tautology) by a diameter each time. Copy this result to the other 3 quadrants once more.

5. Move to the other test points on that level. Do this by using the simple pixel mapping that I used in the image printer. Move to the next X value (curPx + 2b), and then test the 8 points around the cube, (x+b, y+b, z+b), (x-b, y+b, z+b), (x-b, y-b, z+b)... etc. If the distance to any of these points from (0, 0, 0) is greater than the radius of the larger containing circle, the sphere will not be created and the testing will move to the next Y row, also testing this cube's points. When the Y test fails, move down to the next level. Copy these results to the other 3 quadrants.

6. Note, when initial testing of the next level fails, end constructing the cube lattice.

    How to move the point (easy):This is basic modulus and common factors. The crests of each of the frequencies line up
1. Start at (0, 0, 0). Work out the force vectors (F = GM/r^2), and split this into the X, Y, Z values (so each X, Y, Z value will be true for F^2 = X^2 + Y^2 + Z^2), load these values into array. Force[testpoint][x] = F.

2. Once the array has been filled, work out the sum of the X, Y, Z values for force, the result will be the next force. Return F = X^2 + Y^2 + Z^2.

    GRAPH RESULTS:
1. X-axis is the distance from the centre of the sphere, Y-axis is the resultant net force.

Now I gotta write this son of a bitch...
*/

function cubeLatticeArray($a, $b){
  //$a is the radius of the larger sphere, $b is the radius of the smallers spheres that are fitting into the larger sphere.
  //This function will only solve for one eighth of the cube. I will later have this function copy these results over
  //to the other 7 parts of the sphere.
  $radius = pow($a, 2);
    $points[0] = array('x' => $b, 'y' => $b, 'z' => sqrt(pow($a, 2)-pow(2*$b, 2))-$b); //I solved this line for myself, stfu
    $pass = true;
    for($i = 1; $pass == true; $i++){
        $current = end($points);
        if($radius > (pow($current['x']+$b, 2)+pow($current['y']+3*$b, 2)+pow($current['z']+$b, 2))){
            //Move to the next Y cube relative to the last exact position
            $points[$i] = array('x' => $current['x'], 'y' => $current['y']+2*$b, 'z' => $current['z']);
        }elseif($radius > (pow($current['x']+3*$b, 2)+pow(2*$b, 2)+pow($current['z']+$b, 2))){
            //Move to the next X cube relative to the last Z and X position (not Y)
            $points[$i] = array('x' => $current['x']+2*$b, 'y' => $b, 'z' => $current['z']);
        }elseif($radius > (pow(2*$b, 2)+pow(2*$b, 2)+pow($current['z']-3*$b, 2)) && $current['z']-3*$b >= 0){
            //Move to the next Z cube relative to the last Z value (not x or y)
            $points[$i] = array('x' => $b, 'y' => $b, 'z' => $current['z']-2*$b);
        }else{
            $pass = false;
            $total = $i;
        }
    }
    foreach($points as $key => $value){
        $total++;
        $points[$total]['x'] = -$value['x'];
        $points[$total]['y'] = -$value['y'];
        $points[$total]['z'] = -$value['z'];
        $total++;
        $points[$total]['x'] = -$value['x'];
        $points[$total]['y'] = $value['y'];
        $points[$total]['z'] = -$value['z'];
        $total++;
        $points[$total]['x'] = -$value['x'];
        $points[$total]['y'] = -$value['y'];
        $points[$total]['z'] = $value['z'];
        $total++;
        $points[$total]['x'] = $value['x'];
        $points[$total]['y'] = -$value['y'];
        $points[$total]['z'] = -$value['z'];
        $total++;
        $points[$total]['x'] = -$value['x'];
        $points[$total]['y'] = $value['y'];
        $points[$total]['z'] = $value['z'];
        $total++;
        $points[$total]['x'] = $value['x'];
        $points[$total]['y'] = -$value['y'];
        $points[$total]['z'] = $value['z'];
        $total++;
        $points[$total]['x'] = $value['x'];
        $points[$total]['y'] = $value['y'];
        $points[$total]['z'] = -$value['z'];
    }
    return $points;
}

function gravityAtPoints($pos, $array, $radius, $mass, $total){
  //$pos is the distance away from the origin, $array is the array of spheres, $radius is the radius of the spheres, $mass is the total mass of the planet
  //$total is the total amount of spheres.
  $g = 6.6742*pow(10, -11);
  print $g."\n";
  $sphereMass = $mass/$total;
  foreach($array as $key => $value){
    if($value['x'] < $pos){
      $times = -1;
    }else{
      $times = 1;
    }
    $difx = abs($pos-$value['x']);
    $dify = abs($value['y']);
    $difz = abs($value['z']);
    $dist = pow($difx, 2)+pow($dify, 2)+pow($difz, 2);
    $force = ($g*$sphereMass)/$dist;
    $tanh = atan($dify/$difx);
    $xforce = $force*cos($tanh);
    //print "difx: $difx dify: $dify difz: $difz dist: $dist xforce: $xforce \n";
    $sumForce += $xforce*$times;
  }
  return $sumForce;
}


$radius = 10;
$sphereSize = 200;
$mass = 5.9736*pow(10, 24);
$array = cubeLatticeArray($sphereSize, $radius);

end($array);
$total = key($array);

var_dump($array);
$pos = 0;


//THE FILE PRINTER STARTS HERE

echo "Name of Your Function\n"; // :D

//Your parameters
$domain = array("neg" => 1, "pos" => 300); //Please don't use decimal values. :\
$lineThickness = 3; //Decimals are cool here.
$accuracy = 0.1; //lower the number = higher the accuracy, longer the processing time, larger the image. Please only use values of 10^x

//---Begin main program
$maxheight = 0;
$minheight = 0;
for($n = -$domain["neg"]/$accuracy; $n <= $domain["pos"]/$accuracy; $n++){
  $j++;
  $i = gravityAtPoints(abs($n*$accuracy), $array, $radius, $mass, $total)/-10000000; //testing how many steps that it will take for the number $n to reach 1.
  if($i > $maxheight || $j == 0){ //if the result is larger than all the  rest, the new result will become the largest result.
    $maxheight = $i;
  }
  if($i < $minheight || $j == 0){
    $minheight = $i;
  }
  echo $n*$accuracy." = $i\n";
  $x[$j] = $i; //$n is the position along the x-axis, $i is the position along the y-axis.
}
$x[0] = $x[1]; //<same as the line below me>
$x[$j+1] = $x[$j]; //This is because of line 72, where $dif is set. It will try to reference an array with no entry, so dif would otherwise become $x[$CurPx]
//var_dump($x);

//---begin writing file
unlink("image.ppm"); //file that we will be editing, we are are deleting the file so we have a clean file to work on
$fh = fopen("image.ppm", 'a+'); //creates the image file that we will be editting and opens the file for appending

$domain["dif"] = $domain["neg"]+$domain["pos"];
$width = floor($domain["dif"]/$accuracy); //width of the image is the largest number that was tested
$height = floor(abs($minheight - $maxheight)+1); //height is the largest amount of steps that it took to reach zero

$startWrite = "P3\n# image.ppm\n".$width." ".$height."\n255"; //fomatting the ppm file
fwrite($fh, $startWrite);

$j = 0;
for($k = -$domain["neg"]/$accuracy; $k < $width*$height; $k++){
  $CurPx = $j%$width; //basic modulus to decide when a new line must be made
  if($j%5 == 0){ //basic formatting decides when a new line must be placed
    $Nl = " \n";
  }else{
    $Nl = " ";
  }
  if($CurPx == 0){
    $CurPy++;
  }
  $dist = $x[$CurPx] - ($maxheight-$CurPy+1);
  $dif = abs($x[$CurPx]-$x[$CurPx+1])+$lineThickness;
  if($dist > -$dif && $dist < $dif){ //basic algebra, fuck I hate zero being counted as an integer
    $red = "255"; //if the if statement above reflects the array which was made on line 23, the pixel will be white
    $green = "255";
    $blue = "255";
  }elseif($CurPy >= $maxheight+1){ //if Y >= 0
    $clr = str_pad(abs(255-abs(floor(255*(($CurPy-$maxheight-2)/$minheight)))), 3, "0", STR_PAD_LEFT);
    if(0 == $CurPx%($width/$domain["dif"])){
      $red = "128";
      $green = "128";
      $blue = "128";
    }elseif($CurPx <= $domain["neg"]/$accuracy){
      $blue = $clr;
      $red = $clr;
      $green = "000";
    }else{
      $green = $clr;
      $blue = "000";
      $red = "000";
    }
  }elseif($CurPy < $maxheight+1){ //if Y < 0
    $clr = str_pad(abs(floor(255*($CurPy/$maxheight))), 3, "0", STR_PAD_LEFT);
    if(0 == $CurPx%($width/$domain["dif"])){
      $red = "128";
      $green = "128";
      $blue = "128";
    }elseif($CurPx <= $domain["neg"]/$accuracy){
      $blue = $clr;
      $red = "000";
      $green = "000";
    }else{
      $red = $clr;
      $blue = "000";
      $green = "000";
    }
  }
  $write = "$Nl$red $green $blue"; //the line that is to be outputted
    fwrite($fh, $write); //write the pixel to the file.
    $j++;
}
fclose($fh); //close the file, done!
?>