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

#include "ballistics.h"
#include "constants.h"
#include "atragmx.h"
#include "ace.h"
#include "vector.h"
#include "misc.h"

ATragMXSolution *ATragMXCalcSolution(ATragMXContext *ctx, ace_config *ace) {

	//
	// init section
	// 

	ATragMXSolution *sln = NULL;

	if ((sln = calloc(1, sizeof(ATragMXSolution))) == NULL) {
		perror("ATragMXSolution: malloc failed");
		return NULL;
	}

	// local data

	struct {
		double tx;
		double tz;
        Vec3 last_pos;
        Vec3 pos;
        Vec3 velocity;
        Vec3 accel;
		double speed;
		Vec3 gravity;
		double deltaT;
    } bullet = {
        .gravity = {
			.x = 0.0,
			.y = sin(ctx->gun.scope_angle + ctx->target.inclination_angle)
				  * -GRAVITY,
			.z = cos(ctx->gun.scope_angle + ctx->target.inclination_angle)
				  * -GRAVITY
		},
        .deltaT = 1.0 / ctx->sim_steps,
    };

	struct {
		double elevation;
		Windage windage;
		double lead;
		double time_of_flight;
		Vec3 true_velocity;
		double true_speed;
		double v_coriolis;
		double v_coriolis_deflection;
		double h_coriolis;
		double h_coriolis_deflection;
		double spin_drift;
		double spin_deflection;
	} local_solution = {0};

    struct {
		size_t n;
		double range;
		double true_range;
		double range_factor;
		RangeCard range_card;
    } ctx_range = {
		.n = 0,
		.range_factor = ctx->store_range_card ? 1.0936133 : 1,
		.range_card.records = (size_t) (rangeCardEndRange / rangeCardIncrement),		
		.range_card.card = NULL,
    };

	ctx_range.range_card.card = calloc(ctx_range.range_card.records,
									   sizeof(RangeCardRecord));

	Vec3 wind1 = {
		.x = cos(270 - ctx->atmosphere.wind_dir * 30)
			  * ctx->atmosphere.wind_speed[0],
		.y = sin(270 - ctx->atmosphere.wind_dir * 30)
			  * ctx->atmosphere.wind_speed[0],
		.z = 0
	};

	Vec3 wind2 = {
		.x = cos(270 - ctx->atmosphere.wind_dir * 30)
			  * ctx->atmosphere.wind_speed[1],
		.y = sin(270 - ctx->atmosphere.wind_dir * 30)
			  * ctx->atmosphere.wind_speed[1],
		.z = 0
	};

	double wind_drift = 0.0;

    if (ace->ballistics->enabled)
        ctx->gun.ballistic_coef
            = atmospheric_correction(ctx->gun.ballistic_coef,
									 ctx->atmosphere.temperature,
									 ctx->atmosphere.barometric_pressure,
									 ctx->atmosphere.relative_humidity,
									 ctx->atmosphere.atmosphere_model);

    double eoetvoes_multiplier = 0.0;
    if (ace->ballistics->enabled) {
        eoetvoes_multiplier = 2
							  * (0.0000729 * ctx->gun.muzzle_velocity / -GRAVITY)
							  * cos(ctx->target.latitude)
							  * sin(ctx->target.direction);
	}
	
    bullet.pos.x = 0.0;
    bullet.pos.y = 0.0;
    bullet.pos.z = -(ctx->gun.bore_height / 100.0);

	bullet.velocity.x = 0.0;
	bullet.velocity.y = cos(ctx->gun.scope_angle) * ctx->gun.muzzle_velocity;
	bullet.velocity.z = sin(ctx->gun.scope_angle) * ctx->gun.muzzle_velocity;

    while (local_solution.time_of_flight < 15
		   && bullet.pos.y < ctx->target.range) {
		
		// calculate magnitude of bullet.velocity
		bullet.speed = vec3_magnitude(bullet.velocity);
		local_solution.true_velocity = vec3_diff(bullet.velocity, wind1);

		local_solution.true_speed = vec3_magnitude(local_solution.true_velocity);

		if (ace->ballistics->enabled) {
			double data = retard(ctx->gun.drag_model,
								 ctx->gun.ballistic_coef,
								 ctx->gun.muzzle_velocity,
								 ctx->atmosphere.temperature);

			bullet.accel = vec3_mul(local_solution.true_velocity,
									(-1.0 * data));
		} else {
			bullet.accel = vec3_mul(local_solution.true_velocity,
									local_solution.true_speed
									* ctx->gun.air_friction);
		}

		bullet.accel = vec3_add(bullet.accel, bullet.gravity);

		bullet.last_pos = bullet.pos;

		bullet.pos = vec3_add(bullet.pos,
							  vec3_mul(bullet.velocity, bullet.deltaT * 0.5));

        bullet.velocity = vec3_add(bullet.velocity,
								   vec3_mul(bullet.accel, bullet.deltaT));

		bullet.pos = vec3_add(bullet.pos,
							  vec3_mul(bullet.velocity, bullet.deltaT * 0.5));

		sln->time_of_flight = sln->time_of_flight + bullet.deltaT;

		// continue here
        if (ctx->store_range_card) {
			ctx_range.range = rangeCardStartRange
							  + ctx_range.n
								* rangeCardIncrement;

			if (bullet.pos.x * ctx_range.range_factor >= ctx_range.range
				&& ctx_range.range <= rangeCardEndRange) {
				
				ctx_range.true_range = ctx_range.range / ctx_range.range_factor;

				if (ctx_range.true_range != 0) {
					bullet.tx = linearConversion(bullet.last_pos.y,
												 bullet.pos.y,
												 ctx_range.true_range,
												 bullet.last_pos.x,
												 bullet.pos.x, false);

					bullet.tz = linearConversion(bullet.last_pos.y,
												 bullet.pos.y,
												 ctx_range.true_range,
												 bullet.last_pos.z,
												 bullet.pos.z, false);

					local_solution.elevation = - atan(bullet.tz
													  / ctx_range.true_range);

					local_solution.windage.w1 = - atan(bullet.tx
													   / ctx_range.true_range);

					wind_drift = wind2.x * (sln->time_of_flight
											- ctx_range.true_range
											  / ctx->gun.muzzle_velocity);

					local_solution.windage.w2 = - atan(wind_drift
													   / ctx_range.true_range);

					sln->lead = (ctx->target.speed * sln->time_of_flight)
								/ (tan(MRAD_TO_DEG(1))
								   * bullet.speed * bullet.speed);
				}
				double kinetic_energy = 0.5 * (ctx->gun.bullet_mass / 1000.0
											   * (bullet.speed * bullet.speed));

				kinetic_energy *= 0.737562149;

                if (ace->ballistics->enabled && bullet.pos.y > 0) {
					// Coriolis
					local_solution.h_coriolis_deflection =
						0.0000729 * bullet.pos.y * local_solution.time_of_flight
						* sin(ctx->target.latitude);

					local_solution.h_coriolis =
						- atan(local_solution.h_coriolis_deflection
							   / bullet.pos.y);

					local_solution.windage.w1 += local_solution.h_coriolis;
					local_solution.windage.w2 += local_solution.h_coriolis;

                    // Eoetvoes
                    local_solution.v_coriolis_deflection =
						bullet.pos.z * eoetvoes_multiplier;

					local_solution.v_coriolis =
						- atan(local_solution.v_coriolis_deflection
							   / bullet.pos.y);

					local_solution.elevation += local_solution.v_coriolis;

                    // Spin drift
                    local_solution.spin_deflection =
						ctx->gun.twist_dir * 0.0254 * 1.25
						* (ctx->stability_factor + 1.2)
						* pow(local_solution.time_of_flight, 1.83);

					local_solution.spin_drift =
						- atan(local_solution.spin_deflection / bullet.pos.y);

					local_solution.windage.w1 += local_solution.spin_drift;
					local_solution.windage.w2 += local_solution.spin_drift;
                }

				if (ctx_range.n < ctx_range.range_card.records) {

					(ctx_range.range_card.card+ctx_range.n)->n = ctx_range.n;
					(ctx_range.range_card.card+ctx_range.n)->range = ctx_range.range;

					(ctx_range.range_card.card+ctx_range.n)->elevation
						= local_solution.elevation;

					(ctx_range.range_card.card+ctx_range.n)->windage.w1
						= local_solution.windage.w1 * 60.0;

					(ctx_range.range_card.card+ctx_range.n)->windage.w2
						= local_solution.windage.w2 * 60.0;

                    (ctx_range.range_card.card+ctx_range.n)->lead = sln->lead;

                    (ctx_range.range_card.card+ctx_range.n)->time_of_flight
						= local_solution.time_of_flight;

                    (ctx_range.range_card.card+ctx_range.n)->bullet_speed = bullet.speed;
                    (ctx_range.range_card.card+ctx_range.n)->kinetic_energy = kinetic_energy;
				}
				ctx_range.n++;
			}
        }
	}

    double kinetic_energy = 0.5 * (ctx->gun.bullet_mass / 1000.0 *
                                       (bullet.speed * bullet.speed));
    kinetic_energy *= 0.737562149;

	if (ace->ballistics->enabled && bullet.pos.y > 0) {
        // Coriolis
        local_solution.h_coriolis_deflection = 0.0000729 * bullet.pos.y *
                                               local_solution.time_of_flight *
                                               sin(ctx->target.latitude);

        local_solution.h_coriolis =
            -atan(local_solution.h_coriolis_deflection / bullet.pos.y);

        local_solution.windage.w1 += local_solution.h_coriolis;
        local_solution.windage.w2 += local_solution.h_coriolis;

        // Eoetvoes
        local_solution.v_coriolis_deflection =
            bullet.pos.z * eoetvoes_multiplier;

        local_solution.v_coriolis =
            -atan(local_solution.v_coriolis_deflection / bullet.pos.y);

        local_solution.elevation += local_solution.v_coriolis;

        // Spin drift
        local_solution.spin_deflection =
            ctx->gun.twist_dir * 0.0254 * 1.25 * (ctx->stability_factor + 1.2) *
            pow(local_solution.time_of_flight, 1.83);

        local_solution.spin_drift =
            -atan(local_solution.spin_deflection / bullet.pos.y);

        local_solution.windage.w1 += local_solution.spin_drift;
        local_solution.windage.w2 += local_solution.spin_drift;
    }

	sln->elev = local_solution.elevation * 60.0;
	sln->windage.w1 = local_solution.windage.w1 * 60.0;
	sln->windage.w2 = local_solution.windage.w2 * 60.0;
	sln->lead = local_solution.lead;
	sln->time_of_flight = local_solution.time_of_flight;
	sln->remaining_velocity = bullet.speed;
	sln->remaining_energy = kinetic_energy;
	sln->coriolis_vdrift = local_solution.v_coriolis * 60.0;
	sln->coriolis_hdrift = local_solution.h_coriolis * 60.0;
	sln->spin_drift = local_solution.spin_drift * 60.0;
    sln->range_card = &ctx_range.range_card;

    return sln;
}