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

//#include "stdafx.h"
//cc dumka3.c dumkaC.c -lm
#include <stdio.h>
#include "math.h"
FILE *solution;
#ifndef max
#define max(a,b)            (((a) > (b)) ? (a) : (b))
#endif
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*		Linear compementary problem for American oprions 
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*
*   This program solves a system of ODEs resulting from the 
*    space discretization of 
*   the Black-Sholes equations for american Call:
*    
*   (u_t + 0.5*sigma^2*x^2*u_xx + (r-D_0)*x*u_x - ru) >= 0 ("=" - For European) 
*   (u(x,t) - CallPO(x)) >= 0
*   (u_t + 0.5*sigma*x^2*u_xx + (r-D_0)*x*u_x - ru)*(u(x,t) - CallPO(x)) = 0
*   Payoff function CallPO(x) = max(x-E,0) (initial condition)
*  
*   1. Conservative form of the BS equation:
*   u_t + (0.5*sigma^2*x^2*u_x +(r-D_0-sigma^2)*x*u)_x -
*              -(r-D_0-sigma^2)*u - ru >= 0 ("=" - For European)
*
* The algorithm dumka3 and dumka4 is described in the article:
* Medovikov A.A. High order explicit methods for parabolic 
*                equations. BIT, V38,No2,pp.372-390,1998 
*
* Many thanks Prof. Lebedev V.I., Prof. Wanner G. and Prof. Hairer E.
* and Prof. A. Abdulle
*/

//parameters of the problem
double E = 10.e0;
double D0 = 0.2e0;
double r = 0.25e0;
double sigma = 0.8e0;
double interval = 30.e0;
long nPoints = 99;

void rhs(long neqn, double* y, double t, double* result ) {
	int i;
	// ---------- discretisation  -------
	// this is not good - non conservative scheme - change it
	// 0.5*sigma^2*x^2*u_xx + (r-D_0)*x*u_x - ru
	double coef1, coef2, coef_Lft, coef_Ctr, coef_Rgt;
	double y_Lft,y_Rgt;
	double x;
	double dx = interval/(double) (nPoints+1);
	// if american option uncomment lines below y=max(y,Payoff)
	/*
	for ( i=0; i < neqn; ++i) {
		x=dx*((double)(i));
		double a = max(x-E,0.e0);
		double b = y[i];
    	y[i]=max(b,a);
	}
	*/

    // conservative form
	// (0.5*sigma^2*x^2*u_x +(r-D_0-sigma^2)*x*u)_x -
    //          -(r-D_0-sigma^2)*u - ru	
	double sigma2=sigma*sigma;
	coef1 = 0.5e0*sigma2/(dx*dx);
	coef2 = 0.5e0*(r-D0-sigma2)/(dx);
	double x_05,x05;
	for ( i=0; i < neqn; i++ ) {
		x_05=dx*((double)i+0.5e0);
		x05=dx*((double)i+1.5e0);
		x   =0.5e0*(x_05+x05);
		coef_Lft = coef1*x_05*x_05-coef2*x_05;
	    coef_Ctr = -coef1*(x_05*x_05 + x05*x05) + coef2*(x05 - x_05) - (r-D0-sigma2)-r;
        coef_Rgt = coef1*x05*x05+coef2*x05;
		
		if ( i==0 ) {
			y_Lft = 0.e0;
		} else {
			y_Lft = y[i-1];
		}
		if ( i==neqn-1 ) {
			//y_Rgt=x+dx-E*exp(-r*t);
			y_Rgt=max(x+dx-E,0.e0)*exp(-r*t);// // this is test boundary condition fixx it
			//y_Rgt= y[i];
			//y_Rgt= max(x+dx-E,0.e0);
		} else {
			y_Rgt = y[i+1];
		}
		result[i]=(coef_Lft*y_Lft+coef_Ctr*y[i]+coef_Rgt*y_Rgt);
	}
	
};
/*
*    The subroutine cour gives an estimation of the spectral 
*    radius of the Jacobian matrix of the problem rho and returns
*    2/rho, which is tha maximum possible step in explicit Euler method.
-
*/			
double cour(double* y, double t, 
			void (f(long neqn, double* y, double t, double* result))) {
				return (2.e0/1.134250e+004);
			};
//int _tmain(int argc, _TCHAR* argv[])
int main(int argc, char** argv)
{
	void dumka3(long neqn, double* t, const double tend, const double h0,
		const double atol, const double rtol, 
		void (f(long neqn, double* y, double t, double* result)),
		double (cour(double* y, double t, void (f(long neqn, double* y, double t, double* result)))),
		double* y,double* z0,double* z1,double* z2,
		double* z3,double* oldEigenVector,
		int isComputeEigenValue);
	double y[400],z0[400], z1 [400],z2[400],z3[400], oldEigenVector[400];
	int isComputeEigenValue =(1==0);
	long neqn = nPoints;
	int i;
	double  t = 0.e0, tend = 1.e0, h0 = 0.00001e0, atol = 0.001e0, rtol = 0.001e0;

	// ----- initial values -----
	double x;
	double dx = interval/(double) (nPoints+1);
	solution = fopen( "payoffFunction", "w+" );
	for ( i=0; i < neqn; i++ ) {
		x=dx*((double)(i+1));
		y[i]=max(x-E,0.e0);
		fprintf( solution, "%e %e \n", x, y[i] ); 
	}; 
    fclose( solution );

	// start computation
	dumka3(neqn, &t,  tend, h0,
		atol, rtol, rhs, cour, 
		y,z0,z1,z2, z3, oldEigenVector, isComputeEigenValue);

	// output result
	solution = fopen( "solution", "w+" );
	for ( i=0; i < neqn; i++ ) {
		x=dx*((double)(i+1));
		fprintf( solution, "%e %e \n", x, y[i] ); 
	}; 
	fclose( solution );
	return 0;
}