//This code assumes relatively large diffusion coefficient (I used 0.3)
//if you make it small, than well-known oscillatory effect arise
//due  to  central difference approximartion of the first 
//order derivative, so keep this in mind and if you want small vall, than
//change approximation of the first order derivative
//or use "particle methods"
//also you need "dumka.c" to run this
//depending on your compiler you can attach it in different manner
//the code suppose tobeself-explanatory,rather than very efficient, but if 
//you can sqeese more than 20% efficiency out of it-let me know please.

#include "stdafx.h"//commen this line out if you are not unver Miceosoft Visual Studio
#include <stdio.h>
#include <math.h>

#ifndef max
#define max(a,b)
#endif

FILE *payoffFunction;
FILE *solution;
FILE *solution_exact;

double t_old;
int nPointsX=300;
double x_max=1.5;
double x_min=0.0;
double intervalX=x_max-x_min;
long counter_call_back=0;
double mean_Reverting=0.05e0;
//double r0=0.09e0;//0.09e0;
double strike=0.025e0;
double rate_of_meanreversion=2.0;

double pi=acos(-1.0e0);

//put instanteneouse square of vol here, 
//   for example use 0.3*0.3*x
//   for SDE:
//   dx = a(x,t)dt + 0.3 *sqrt(x)*dW
//
double sigma_sq(double t, double x){
     double sigma=0.3e0;
     return sigma*sigma;
}

//
//right hand side goes here dV/dt = d(x0-x)*v/dx  + d^2 sigma^2*v/d x^2
//
void rhs(long neqn, double* y, double t,double *result){
     double y_Lft, y_Rgt;
     double x;
     double dx=intervalX/(double)  (nPointsX+1);

     double dx_sq_inv=1.0e0/(dx*dx);
     double dx_inv=1.0e0/dx;

     double x_Rgt, x_Left;
     //calculation of the right hand side by method of lines
     //also read comments above about sigma_sq-it got to be large
     //if you sigma is small change approximation of the first derivative 
     //to avoid well-known oscillatory effect
     for( int i=0; i < neqn; i++ ) {
        x=x_min+dx*(double)i;
        x_Rgt=x+dx;
        x_Left=x-dx;
        if( i==0 ) {
          y_Lft=2.0e0*y[i]-y[i+1];//null gamma boundary
        } else  {
          y_Lft=y[i-1];
        }
        if( i==neqn-1 ){
          y_Rgt=2.0e0*y[i]-y[i-1];//zero second derivative boundary
        } else {
          y_Rgt = y[i+1];
        }
        result[i]=rate_of_meanreversion*(mean_Reverting-x)*(y_Rgt- y_Lft)*(0.5e0*dx_inv) + 
                  0.5e0*sigma_sq(t,x)*(y_Rgt - 2.0e0*y[i]+y_Lft)*dx_sq_inv - x*y[i];
     }

     //this part is for output only.
     //If you want efficiency-comment it out
     if ( counter_call_back%1 == 0 && t_old < t ) {
        t_old=t;
        double tau=1.0e0-t;
        fprintf(solution, "%f",tau );
        for ( int m=0; m < neqn; m++ ) {
           x=x_min+dx*(double)m;
           fprintf( solution, " %f",y[m]);
        };
        fprintf(solution,"\n");

	//below is output of exact solution
        fprintf(solution_exact, "%f", tau );
        for ( int i=0; i <neqn; i++ ) {
           double x=x_min+dx*(double)i;
           double B=(1.0e0-exp(-rate_of_meanreversion*t))/rate_of_meanreversion;
		   double rate_of_meanreversion_sq=rate_of_meanreversion*rate_of_meanreversion;
           double A=exp(
                      (B-t)
                      *(rate_of_meanreversion_sq*mean_Reverting- 0.5e0*sigma_sq(t,x) )/
                      rate_of_meanreversion_sq
                      - sigma_sq(t,x)*B/(4.0e0*rate_of_meanreversion_sq)
                      );

           double bond_price= A*exp(-B*x);
           fprintf(solution_exact, " %f", bond_price);
         };
         fprintf(solution_exact,"\n");
     };
     counter_call_back++;
   };

   /*
   *   The subroutine courgives an estimation of the spectral
   *   radius of the Jacobian matrix of the problem rhoand returns
   *   2/rho, which is the maximum step in explicit Euler method.
   *   Another way to estimate -try Euler method with this step
   *   and if it is stable - use this step as return value here
   */
   double cour(double* y, double t,
          void(f(long neqn, double* y, double t, double* result))) {
             return (2.0/250.0e0); //this number is just from experiment
                                    //one can use Gershgorin's theoreme to estimete max_eigenvalue
   };
/*
   double cour_integral(double* y, double t,
           void(f(long neqn, double* y, double t, double* result))) {
              return (2.0e0/(2.0e0*250.0e0));
   }
*/

   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==1);
      long neqn=nPointsX;
      int i;
      double t=0.0e0, tend=1.0e0, h0=0.00001e0, atol=0.001e0, rtol=0.001e0;
      //----- initial values ----//
      double x;
      double dx=intervalX/(double) (nPointsX+1);

      errno_t err = fopen_s( &payoffFunction, "payoffFunation", "w+" );
      errno_t err_final = fopen_s( &solution, "solution", "w+" );
      errno_t err_exact = fopen_s( &solution_exact, "solution_exact", "w+" );

      for ( i=0; i < neqn; i++ ) {
         x=x_min+dx*((double)(i));
         //y[i]=0.0e0;
         y[i]=1.0e0;
         fprintf( payoffFunction, "%f %f \n", x, y[i] );
      };
      fclose( payoffFunction );


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

      fclose( solution );
      fclose( solution_exact );
   };//end of main
   // remember - you need dumka3.c to run this