// Euler_Verlet.cpp : Defines the entry point for the console application.
//

//#include "stdafx.h"
#include <stdexcept>
#include <iostream>
#include <fstream>
using namespace std;

//const int N = 2;
const int N = 8;

//funkcja, która oblicza przyspieszenie
double rzut(int i, double* y, double t)
{
	static const double g = 9.81; //m i s

	double wynik = 0;
	switch(i)
	{
	case 0:
		wynik = y[1];		
		break;
	case 1:
		wynik = g;
		break;

	default:
		throw runtime_error("Zły numer równania");
	}

	return wynik;
}

double oscylator(int i, double* y, double t)
{
	static const double k = 1;

	double wynik = 0;
	switch (i)
	{
	case 0:
		wynik = y[1];
		break;
	case 1:
		wynik = -k*y[0];
		break;

	default:
		throw runtime_error("Zły numer równania");
	}

	return wynik;
}

const double jednostkiMasy = 5.97219E24;
const double jednostkaOdległości = 384403000;
const double jednostkaCzasu = 24 * 60 * 60;

#define SQR(x)((x)*(x))

double grawitacja2(int i, double* y, double t)
//y[0] = x1, y[1] = vx1, y[2] = y1, y[3] = vy1
//y[4] = x2, y[5] = vx2, y[6] = y2, y[7] = vy2
{
	static const double G = 6.6742867E-11; //m^3/kg/s^2
	static const double M = jednostkiMasy;
	static const double m = 0.0123*jednostkiMasy;

	double odległość_x = y[4] - y[0];
	double odległość_y = y[6] - y[2];
	double kwadrat_odległości = SQR(odległość_x) + SQR(odległość_y);
	double odległość = sqrt(kwadrat_odległości);
	double siła = G*M*m / kwadrat_odległości;
	double _x = odległość_x / odległość;
	double _y = odległość_y / odległość;

	double wynik = 0;
	switch (i)
	{
	//Ziemia
	case 0:
		wynik = y[1];
		break;
	case 1:
		wynik = siła *_x / M;
		break;
	case 2:
		wynik = y[3];
		break;
	case 3:
		wynik = siła *_y / M;
		break;

	//Księżyc
	case 4:
		wynik = y[5];
		break;
	case 5:
		wynik = -siła *_x / m;
		break;
	case 6:
		wynik = y[7];
		break;
	case 7:
		wynik = -siła *_y / m;
		break;

	default:
		throw runtime_error("Nieprawidlowy numer rownania");
	}

	return wynik;
}

double* odeint_Euler(int N, double(*f)(int i, double* y, double t), double* y, double t, double h, double* y_nast)
{
	for (int i = 0; i < N; i++) y_nast[i] = y[i] + h*f(i, y, t);
	return y_nast;
}

double* odeint_MidPoint(int N, double(*f)(int i, double* y, double t), double* y, double t, double h, double* y_nast)
{
	double* y_tmp = new double[N];
	double k1, k2;

	for (int i = 0; i < N; ++i)
	{
		k1 = h*f(i, y, t);
		y_tmp[i] = y[i] + 0.5*k1;
	}

	for (int i = 0; i < N; ++i)
	{
		k2 = h*f(i, y_tmp, t + 0.5*h);
		y_nast[i] = y[i] + k2;
	}

	delete[] y_tmp;

	return y_nast;
}

double* odeint_RK4(int N, double(*f)(int i, double* y, double t), double* y, double t, double h, double* y_nast)
{
	double* y_tmp = new double[N];
	double* k1 = new double[N];
	double* k2 = new double[N];
	double* k3 = new double[N];
	double* k4 = new double[N];

	for (int i = 0; i < N; ++i)
	{
		k1[i] = h*f(i, y, t);
		y_tmp[i] = y[i] + 0.5*k1[i];
	}
	for (int i = 0; i < N; ++i)
	{
		k2[i] = h*f(i, y_tmp, t + 0.5*h);
		y_nast[i] = y[i] + 0.5*k2[i];
	}
	for (int i = 0; i < N; ++i)
	{
		k3[i] = h*f(i, y_nast, t + 0.5*h);
		y_tmp[i] = y[i] + k3[i];
	}
	for (int i = 0; i < N; ++i)
	{
		k4[i] = h*f(i, y_tmp, t + h);
		y_nast[i] = y[i] + (k1[i] + 2.0 * k2[i] + 2.0*k3[i] + k4[i]) / 6.0;
	}

	delete[] y_tmp;
	delete[] k1;
	delete[] k2;
	delete[] k3;
	delete[] k4;

	return y_nast;
}

/*
int _tmain(int argc, _TCHAR* argv[])
{
	return 0;
}
*/

int main(int argc, char* argv)
{
	double* y = new double[N];
	double* y_nast = new double[N];

	for (int i = 0; i < N; ++i) y[i] = 0;

	//y[0] = 1;

	//warunki początkowe
	y[6] = jednostkaOdległości;
	y[5] = 1.022E3; //vx2 m/s

	//przejście do układu środka masy
	static const double M = jednostkiMasy;
	static const double m = 0.0123*jednostkiMasy;
	double Vsm_x = (M*y[1] + m*y[5]) / (m + M);
	double Vsm_y = (M*y[3] + m*y[7]) / (m + M);
	y[1] -= Vsm_x;
	y[5] -= Vsm_x;
	y[3] -= Vsm_y;
	y[7] -= Vsm_y;

	double Rsm_x = (M*y[0] + m*y[4]) / (m + M);
	double Rsm_y = (M*y[2] + m*y[6]) / (m + M);
	y[0] -= Rsm_x;
	y[4] -= Rsm_x;
	y[2] -= Rsm_y;
	y[6] -= Rsm_y;

	const double h = 1E-3 * jednostkaCzasu;
	const double tmax = 365 * jednostkaCzasu;

	ofstream plik_wy("..\\Debug\\wyniki.dat");	

	for (double t = 0; t < tmax; t += h)
	{
		//odeint_Euler(N, rzut, y, t, h, y_nast);
		//odeint_Euler(N, oscylator, y, t, h, y_nast);
		//odeint_MidPoint(N, oscylator, y, t, h, y_nast);
		//odeint_RK4(N, oscylator, y, t, h, y_nast);
		odeint_RK4(N, grawitacja2, y, t, h, y_nast);

		//cout << t << "\t";
		plik_wy << t << "\t";
		for (int i = 0; i < N; ++i)
		{
			y[i] = y_nast[i];
			//cout << y[i] << "\t";
			plik_wy << y[i] << "\t";
		}
		//cout << "\n";
		plik_wy << "\n";
	}

	plik_wy.close();

	delete[] y;
	delete[] y_nast;

	cout << "OK.\n\n";

	return 0;
}
