/* 
 * This program is free software; you can redistribute it and/or modify 
 * it under the terms of the GNU General Public License as published by 
 * the Free Software Foundation; either version 2 of the License, or 
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but 
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 
 * for more details.
 * 
 * You should have received a copy of the GNU General Public License along 
 * with this program; if not, write to the Free Software Foundation, Inc., 
 * 59
 */

# include <iostream>
# include "LuxHeader.h"

using namespace std;

void RungeKutta(double y[], double dy[]){

            double k1[var], k2[var], k3[var], k4[var], y2[var], y3[var], y4[var], *pode;
            double h = 0.05;
            pode = ode(y,dy);

            	// Calculates k1 (for all variables) -> The slope at the start of the interval (h)
            	for (int i = 0; i < var; ++i) {
                        k1[i] = *(pode + i) * h;
            		}

           	 //  Calculates k2 (for all variables) -> The slope at the midpoint of the interval (h),
            	for (int i = 0; i < var; ++i) {
                        y2[i] =y[i] + k1[i] * 0.5 ;
            		}
            	
		pode = ode(y2,dy);
            	for (int i = 0; i < var; ++i) {
                        k2[i] = *(pode + i) * h;
            		}

            	// Calculates k3 (for all variables) -> The slope at the midpoint of... using the y value (y2) determined from k2.
            	for (int i = 0; i < var; i++) {
                        y3[i]=y[i] + k2[i] * 0.5;
            		}
            
		pode = ode(y3,dy);
            	for (int i = 0; i < var; i++) {
                        k3[i] = *(pode + i) * h;
            		}

            	// Calculates k4 (for all variables) -> The slope at the end of the interval (h).
            	for (int i = 0; i < var; i++) {
                        y4[i]= y[i] + k3[i];
            		}
            
		pode=ode(y4,dy);
            	for (int i = 0; i < var; i++) {
                        k4[i] = *(pode + i) * h;
            		}
	
            	// Calculates the new y values (for all variables).
            		for (int i = 0; i < var; i++){
                        y[i] +=((k1[i] + 2*k2[i] + 2*k3[i] + k4[i])/6);
            		}
}