#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define G 6.672592e-11         // gravitaional constant
#define ALT 100000.0           // altitude in m (100km)
#define SATMASS 250.0          // satellite mass (250kg)
#define MASS 5.23533e8         // mass of our attractant


double pythag(double a, double b);
double gravforce(double m1, double m2, double r);
double fz(double f, double a, double b);
double getmarez(double x, double y, double z, double m1, double m2);

int main (void)
{

  double satx,saty,fzhat,d,x,y;
  int  loop;
  FILE *out;

  loop=0;
  satx=saty=fzhat=d=x=y=0.0;

  /******************************************************
   * First we will do an iron sphere 1km across located *
   * in the centre of our 200x200km field at (100,100,0)*
   * which is in metres (1d5,1d5,0.0). Such a sphere of *
   * pure iron has a density of 7874 kgm^-3 and so has a*
   * mass of approximately 5.23533d8kg.                 *
   ******************************************************/
 
  out=(FILE *) fopen("iron.dat","wt");
  satx=0.0;
  while (satx<200100.0) // one pass over the field
    {
      saty=0;          // start at y origin

      loop=0;          // set up counter

      while(loop<10)  //go 10 times for 10 paths
	{

	  x=1e5-satx;  // get x displacement
	  y=1e5-saty;  // get y displacement
	  d=pythag(x,y); // get dispacement in xy plane
	  d=pythag(ALT,d); // get total displacement in 3-space

	  fzhat=gravforce(SATMASS,MASS,d); // get the total force
	  fzhat=fz(fzhat,ALT,pythag(x,y)); // get the z componet of the force

	  fprintf(out,"%G",fzhat);  // print result
	  if(loop<10)
	    fprintf(out," ");    //space them out
	  loop++;        // increment counter
	  saty+=(2.0e5/10);// next orbit path
	}
      fprintf(out,"\n"); // new line in file
      satx+=100.0;     // next hundred metres
    }

  fclose(out);

  /**********************************************************
   * Now we will do a 2e5x100m mass of mare represented by  *
   * the mineral ilmenite which has a density of 4790kgm^-3 *
   * which, at a mass of 5.23533d8kg makes it a depth of    *
   * about 5mm! We will compute the force as if the ilmenite*
   * is spread as 10 equally massed and distributed bodies  *
   * in a line across the spacecraft's path. This makes the *
   * math easier.                                           *
   **********************************************************/
 
  out=(FILE *) fopen ("mare.dat","wt");

  satx=0.0;
  while(satx<200100.0)  // one pass over the field
    {
      saty=0.0;         // start at the y origin.
      loop=0;

      while(loop<10)   // go 10 times
	{
	  fzhat=getmarez(satx,saty,ALT,SATMASS,MASS);  // do this in a subroutine since it's complicated
	  fprintf(out,"%G",fzhat); // print result to file
	  if (loop<10)
	    fprintf(out," ");  // space them out
	  saty+=(2.0e5/10);   // next orbit path
	  loop++;             // increment counter

	}
      fprintf(out,"\n");         //new line in file next line
      satx+=100.0;    //next hundred metres
    }
  fclose(out);

}


double pythag (double a, double b)
{

  double c;

  c=(a*a)+(b*b);
  c=sqrt(c);

  return c;

}


double gravforce(double m1, double m2, double r)
{

  double f;

  f=m1m2;
  f=f/(r*r);
  f*=G;

  return f;

}

double fz(double f, double a, double b)
{
  double fzhat;

  fzhat=f/(sin(atan(a/b)));

  return fzhat;

}

double getmarez(double x, double y, double z, double m1, double m2)
{

  /* we're going to do this brute forceish, not particularly elegant math
     basically, compute each of the 10 masses seperately */

  double f, fzhat, xdisp,ydisp,d,mass,massy;
  int loop;

  mass=m2/10.0;  // set our mass at 1/10 total mass
  fzhat=0; // clear our force store
  massy=0; // sart masses at origin
  loop=0; // clear counter
  while(loop<10)
    {

      xdisp=1e5-x;  // get x displacement
      ydisp=massy-y;  // get y displacement
      d=pythag(xdisp,ydisp); // get dispacement in xy plane
      d=pythag(z,d); // get total displacement in 3-space
  
      f=gravforce(m1,mass,d); // get the total force
      fzhat+=fz(f,ALT,pythag(x,y)); // get the z componet of the force
      
      massy+=2e4;     // next mass
      loop++;   //increment counter
    }

  return fzhat; // send result

}