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

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

/* Define the Datapoint data structure
   to hold the easting, northing, and mag values
   for magnetic observations and calculated mag values
   */
typedef struct {
    double northing;
    double easting;
    double mag_obs;
    double mag_calc;
} DataPoint;

/* 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.northing;
    y = y - sphere.easting;

    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;
}

// Function prototype for calculating RMSE
double calculate_rmse(double *obs, double *calc, int size) {
    double sum = 0.0;
    for (int i = 0; i < size; i++) {
        double diff = obs[i] - calc[i];
        sum += diff * diff;
    }
    return sqrt(sum / size);
}

// Function to perform the brute force search
// prints depth and estimated RMSE
double brute_force(DataPoint *mag_anomaly, Sphere sphere, int size, double high, double low) {
	double calc[size];
	double obs[size];
	//printf("%lf %lf\n",mag_anomaly[500].northing, mag_anomaly[500].easting);
	//printf("%lf, %lf\n", sphere.a, sphere.minc);
    for (int j=0; j<100; j++){
    	
    	  // calculate the depth over interval low to high
    	  sphere.z = low + (float) j * (high - low)/(100.0); 
        for (int i = 0; i < size; i++) {
            //mag_anomaly[i].mag_calc = magnetized_sphere_anomaly(mag_anomaly[i].northing, mag_anomaly[i].easting, sphere);
            obs[i] = mag_anomaly[i].mag_obs;
            calc[i] = magnetized_sphere_anomaly(mag_anomaly[i].northing, mag_anomaly[i].easting, sphere);
        }
       // printf("%lf %lf\n", obs[1000], calc[1000]);
        double rmse = calculate_rmse(obs, calc, size);
        printf("%lf %lf\n", sphere.z, rmse);
       }
    return 0;
}


int main () {

    //set parameters (depth will change)
    Sphere sphere1 = {0.0, 0.0, 300.0, 100.0, 45.0, 0.0, 1.0, 45.0, 0.0};
  
    /*read the observed magnetic anomaly data file  */
     // File pointer
    FILE *file = fopen("mag_sphere.dat", "r");
    if (file == NULL) {
        perror("Holy smokes...no magnetic data file found");
        return 1;
    }

    // Array to hold data
    DataPoint mag_anomaly[2000]; // Adjust size as needed
    int count = 0;
    /*set bounds on search*/
    double low = 0.1;
    double high = 1000.0;

    // Read data from file
    while (fscanf(file, "%lf %lf %lf %lf", 
                  &mag_anomaly[count].northing, 
                  &mag_anomaly[count].easting,
                  &mag_anomaly[count].mag_obs,
                  &mag_anomaly[count].mag_calc)
                  == 4) {
        count++;
        if (count > 2000) {printf("max size reached\n");}
    }

    fclose(file);
    
    int result = brute_force(mag_anomaly, sphere1, count, high,low);
    return 0;
}