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

// Define the Sphere data structure
typedef struct {
    double easting;
    double northing;    
    double z;
    double a;
    double minc;
    double mdec;
    double mi;
    double einc;
    double edec;
} Sphere;

/* Function definition for magnetized_sphere_anomaly
   -calculate the magnetic anomaly due to a buried sphere. 
   Parameters  x = an observation point northing coordinate (m)
               y = an observation point easting coordinate (m)
               (x and y are relative to the sphere center, which is assumed to be (0,0))
   Sphere includes
               easting = easting coordinate of the sphere center (m)
               northing = northing coordinate of the sphere center (m)
               z = depth the center of the sphere (m)
               a = radius of the sphere
               minc = inclination of the vector of magnetization (deg)
               mdec = declination of the vector of magnetization (deg)
               mi = intensity of magnetization (amp/m)
               einc = inclination of the earths magnetic field (deg)
               edec = declination of the earths magnetic field (deg)
        return the total anomaly of the magnetic field
*/
double magnetized_sphere_anomaly(double x, double y, Sphere sphere) {
    double minc = sphere.minc * M_PI / 180.0; // Inclination down in rad
    double mdec = sphere.mdec * M_PI / 180.0; // Declination east in rad

    double ml = cos(minc) * cos(mdec);
    double mm = cos(minc) * sin(mdec);
    double mn = sin(minc);

    double dpm = (4.0 / 3.0) * M_PI * pow(sphere.a, 3) * sphere.mi;

    double mx = sphere.mi * ml;
    double my = sphere.mi * mm;
    double mz = sphere.mi * mn;

    double einc = sphere.einc * M_PI / 180.0; // Inclination down in rad
    double edec = sphere.edec * M_PI / 180.0; // Declination east in rad

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

    double prop = 400 * M_PI;

    double z = - sphere.z;
    
    x = x - sphere.easting;
    y = y - sphere.northing;

    double r2 = x * x + y * y + z * z;
    double r = sqrt(r2);
    double r5 = pow(r, 5);

    double dotprod = x * mx + y * my + z * mz;

    double bx = prop * dpm * (3 * dotprod * x - r2 * mx) / r5;
    double by = prop * dpm * (3 * dotprod * y - r2 * my) / r5;
    double bz = prop * dpm * (3 * dotprod * z - r2 * mz) / r5;

    return el * bx + em * by + en * bz;
}

/* This function randomly draws one sample from
   a normal distribution with mean and standard deviation
   using the Box-Muller method
*/
double generate_normal_sample(double mean, double stddev) {
    // Generate two uniform random numbers between 0 and 1
    double u1 = ((double)rand() + 1) / ((double)RAND_MAX + 1);
    double u2 = ((double)rand() + 1) / ((double)RAND_MAX + 1);

    // Apply the Box-Muller transform
    double z0 = sqrt(-2.0 * log(u1)) * cos(2.0 * M_PI * u2);

    // Scale and shift to get the desired mean and standard deviation
    return mean + z0 * stddev;
}



int main () {

    //set parameters
    Sphere sphere1 = {0.0, 0.0, 300.0, 100.0, 45.0, 0.0, 1.0, 45.0, 0.0};
  
    double y = 0.0;
    double sum = 0.0; 
    for (int i=0;i<2000; i++){
    	  double x = (float) i - 1000.0; // x is the horizontal distance across the sphere
        // nt is the magnetic anomaly due to the sphere - with normal noise added
        double nt = magnetized_sphere_anomaly(x,y, sphere1) + generate_normal_sample(0.0,25.0);
        printf("%lf %lf %lf %lf\n", x,y,nt,0.0);
    }
    
   
    return 0;
}