/**
 * @brief motion of a single particle in a homogenous magnetic field simulator 
 * @author Michael Roth <xrothm00@stud.fit.vutbr.cz>
 * @version 1.0.5
 * @since 0.5
*/
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include "physics.h"
#include "cont_sim.h"
#include "simulator.h"

int main(/* int argc, char ** argv */) {

	// initialize the random generator
	srand(RANDOM_MODE ? time(NULL) : 42);

    // comment out experiments that we don't want to execute
    sim1();
    sim2();
    sim3();
	sim4();

    // end the program
    return 0;

} // end main()

int sim1() {

    // create the particle structures
    Particle p1, p2;

    // setup the particles
    p1.id = 1;
    p2.id = 2;

    p1.time = 0.0;
    p2.time = 0.0;
    p1.end_time = 100.0;
    p2.end_time = 100.0;
    
    p1.m_0 = 1.0;
    p2.m_0 = 1.0;
    p1.charge = 1.0;
    p2.charge = 1.0;
    p1.st[I_VX] = 1.0;
    p2.st[I_VX] = 1.0;

    // calculate the relativistic constants for particles
    calculate_gqm(&p1);
    calculate_gqm(&p2);

    // setup the output .csv file 
    FILE *f1 = fopen("sim1_rk4.csv", "w");
    if(f1 == NULL) return -1;

    fprintf(f1, "time,x,y,z,v_x,v_y,v_z,id\n");

    // main loop of simulation, particle 1
    simulation_loop(p1, STEP_SIZE, f1, NULL, rk4);
    
    // close the output file 
    fclose(f1);

    // setup the output file
    FILE *f2 = fopen("sim1_euler.csv", "w");
    if(f2 == NULL) return -1;
    fprintf(f2, "time,x,y,z,v_x,v_y,v_z,id\n");

    // main loop of simulation
    simulation_loop(p2, STEP_SIZE, f2, NULL, euler);

    // close the output file
    fclose(f2);

    // end the program
    return 0;

} // end sim1

int sim2() {

    Particle p1, p2, p3;

    p1.id = 1; // proton-ish particle
    p2.id = 2; // neutron-ish particle
    p3.id = 3; // electron-ish particle

    p1.time = 0.0;
    p2.time = 0.0;
    p3.time = 0.0;

    p1.end_time = 2;
    p2.end_time = 2;
    p3.end_time = 2;

    p1.m_0 = 1;
    p2.m_0 = 1; 
    p3.m_0 = 0.1;

    p1.charge = 1;
    p2.charge = 0;
    p3.charge = -1;

    p1.st[I_VY] = 1;
    p2.st[I_VY] = 1;
    p3.st[I_VY] = 1;
    p1.st[I_VX] = 0;
    p2.st[I_VX] = 0;
    p3.st[I_VX] = 0;
    p1.st[I_VZ] = 0;
    p2.st[I_VZ] = 0;
    p3.st[I_VZ] = 0;

    calculate_gqm(&p1);
    calculate_gqm(&p2);
    calculate_gqm(&p3);

    // setup the output files
    FILE * traj1 = fopen("sim2_trajectory_proton.csv", "w");
    FILE * traj2 = fopen("sim2_trajectory_neutron.csv", "w");
    FILE * traj3 = fopen("sim2_trajectory_electron.csv", "w");
    
    if(traj1 == NULL || traj2 == NULL || traj3 == NULL) return -1;

    fprintf(traj1, "time,x,y,z,v_x,v_y,v_z,id\n");
    fprintf(traj2, "time,x,y,z,v_x,v_y,v_z,id\n");
    fprintf(traj3, "time,x,y,z,v_x,v_y,v_z,id\n");

    // main loop of simulation
    simulation_loop(p1, STEP_SIZE / 5, traj1, NULL, rk4);
    simulation_loop(p2, STEP_SIZE / 5, traj2, NULL, rk4);
    simulation_loop(p3, STEP_SIZE / 5, traj3, NULL, rk4);

    // close the output files
    fclose(traj1);
    fclose(traj2);
    fclose(traj3);

    // end the program
    return 0;

} // end sim2.

int sim3() {
    Particle p;
    
    p.id = 1;
    p.time = 0.0;
    p.end_time = 200.0;

    p.m_0 = 1.0;
    p.charge = 1.0;
    p.st[I_VX] = 1.0;
    p.st[I_VY] = 0.0;
    p.st[I_VZ] = 0.02;

    calculate_gqm(&p);

    FILE *traj = fopen("sim3_trajectory.csv", "w");
    FILE *detector = fopen("sim3_detector_log.csv", "w");

    if(traj == NULL || detector == NULL) return -1;

    // setup the output files
    fprintf(traj, "time,x,y,z,v_x,v_y,v_z,id\n");
    fprintf(detector, "id,time,x,y,z,v_x,v_y,v_z,charge\n");

    // main loop of simulation
    simulation_loop(p, STEP_SIZE, traj, detector, rk4);

    // close the output files
    fclose(traj);
    fclose(detector);

    // end the program
    return 0;

} // end sim3

int sim4() {
	// step set to one picosecond
	double step = pow(10.0, -12);
	
	// demo particle
	//Particle p;

	Particle particles[N_PARTICLES];

	for(int i = 0; i < N_PARTICLES; i++) {
		// set the id of each particle
		particles[i].id = i + 1;

		// set time and end time of each particle
		particles[i].time = 0.0;
		particles[i].end_time = 15 * pow(10.0, -9);

		// randomly pick, which of 4 select particles is created
		// and assign invariant mass and charge accordingly
		switch(rand() % 4) {
			case 0: // case proton
				particles[i].m_0 = PROTON_M0;
				particles[i].charge = PROTON_Q * ELEMENTARY_CHARGE;
				break;
			case 1: // case charmed lambda
				particles[i].m_0 = LAMBDA_CHARM_M0;
				particles[i].charge = LAMBDA_CHARM_Q * ELEMENTARY_CHARGE;
				break;
			case 2: // case sigma minus
				particles[i].m_0 = SIGMA_MINUS_M0;
				particles[i].charge = SIGMA_MINUS_Q * ELEMENTARY_CHARGE;
				break;
			case 3: // case bottom sigma
				particles[i].m_0 = SIGMA_BOTTOM_MINUS_M0;
				particles[i].charge = SIGMA_BOTTOM_MINUS_Q * ELEMENTARY_CHARGE;
				break;
			default:
				break;
		} // end switch

		// create random velocity vector
		double radius = LIGHT_SPEED * rand_float(0.8, 0.95); // speed from 0,8c to 0,9c
		double azimuth = rand_float(5, 7) * PI / 12; // azimuth from PI/3 to 2PI/3
		double inclination = rand_float(5,7) * PI / 12; // inclination from PI/3 to 2PI/3

		// calculate the remaining parameters
		set_velocity_from_spherical(&particles[i], radius, azimuth, inclination);
		calculate_gqm(&particles[i]);
	} // end for

	FILE* traj = fopen("sim4_trajectory.csv", "w");
	FILE* detector = fopen("sim4_detector_log.csv", "w");

	if (traj == NULL || detector == NULL) return -1;

	// setup the output files
	fprintf(traj, "time, x, y, z, v_x, v_y, v_z, id\n");
	fprintf(detector, "id, time, x, y, z, v_x, v_y, v_z, charge\n");

	//simulation_loop(p, step, traj, detector, rk4);

	// main loop of simulation
	for (int i = 0; i < N_PARTICLES; i++) {
		printf("p%d: m_o = %10e, q = %10e, v = %10e\n", 
			particles[i].id, particles[i].m_0, 
			particles[i].charge, get_speed(&particles[i]));
		simulation_loop(particles[i], step, traj, detector, rk4);
	} // end for

	// close the output files
	fclose(traj);
	fclose(detector);

	// end the program
	return 0;
} // end sim4

void set_velocity_from_spherical(Particle *p, double radius, double azimuth, double inclination) {
	// TODO make a graph of the reference axes
	p->st[I_VX] = radius * sin(azimuth) * sin(inclination);
	p->st[I_VY] = radius * cos(azimuth) * sin(inclination);
	p->st[I_VZ] = radius * cos(inclination);
} // end velocity_from_spherical

double rand_float(double min, double max) {
	return (double) rand() / (double)(RAND_MAX / (max - min)) + min;
} // end rand_float

double get_speed(Particle* p) {
	return sqrt(pow(p->st[I_VX], 2) + pow(p->st[I_VY], 2) + pow(p->st[I_VZ], 2));
} // end get_speed