//  d B(r0,lambda0,t)/dt = mean_reversion_r(r0-r)*dB/dr + 0.5*sigma_r^2*d^2B/dr^2 + 
//              0.5*sigma_lambda^2*d^2B/dlambda^2 + rho*sigma_r*sigma_lambda*d^2B/dr*dlambda 
//              -(r+lambda)*B
//              g(t,r,lambda)*lambda + f(t,r,lambda)
// where      g(t,r,lambda)   - expected value of integral 
//            f(t,r,lambda) -   continuouse source/sink
//          
//
//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=-0.1e0;

int nPointsX=200;
int nPointsY=200;

double x_min=-0.0;
double x_max=1.0;

double l_min=-0.0;
double l_max=1.0;


double intervalX=x_max-x_min;
double intervalL=l_max-l_min;
double dx=intervalX/(double)  (nPointsX-1);
double dl=intervalL/(double)  (nPointsY-1);

long counter_call_back=0;
double mean_Reverting_r=0.05e0;
double mean_Reverting_l=0.05e0;
double rho = -0.2;

double r0=0.09e0;//0.09e0;
double rate_of_meanreversion_r=2.0;
double rate_of_meanreversion_l=2.0;
double *x_arr, *l_arr;

double pi=acos(-1.0e0);
double final_payoff(double T, double r, double l)
{
      return 1.0e0;
}
double continuouse_payoff(double t, double r, double l)
{
      return 0.0e0;
}
double default_payoff(double t, double r, double l, double pi)
{
      return 0.0e0;
}
double expected_default_payoff(double t, double r, double l)
{
      return 0.0e0;
}
double sigma_r(double t, double x){
     double sigma=0.3e0;
     return sigma;
}
double sigma_l(double t, double l){
     double sigma=0.3e0;
     return sigma;
}


void rhs(long neqn, double* y, double t,double *result){
     //double y_Lft, y_Rgt, y_Top, y_Bottom, y_Center;
     double sigma_r_tmp, sigma_l_tmp;
     double drift_r=0.0e0;
     double drift_l=0.0e0;
     

     double dx_sq_inv=1.0e0/(dx*dx);
     double dl_sq_inv=1.0e0/(dl*dl);
     double dx_inv=1.0e0/dx;
     double dl_inv=1.0e0/dl;
     double dx_dl_inv=dx_inv*dl_inv;


//---------------   not finished boundary  conditions yet -----------------
     int i,j,k;
     i=0;
     for ( j=0; j<nPointsY; j++ ) {// rate=x_min
        k=j*nPointsX +i;
        y[k]=y[k+1];//change boundary conditions here dy/dr=0
     }
     i=nPointsX-1;
     for ( j=0; j<nPointsY; j++ ) {// rate=x_max
        k=j*nPointsX +i;
        y[k]=y[k-1];//change boundary conditions here dy/dr=0
     }
     j=0;
     for ( i=0; i<nPointsY; i++ ) {// default rate=l_min
        k=j*nPointsX +i;
        y[k]=exp(-x_arr[i]*t);
         /*
		   double x= x_arr[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);
	 */
     }
     j=nPointsY-1;
     for ( i=0; i<nPointsY; i++ ) {// default rate=l_max
        k=j*nPointsX +i;
        y[k]=expected_default_payoff(t,x_arr[i],l_arr[j]);
     }


     //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 ( i=1; i<nPointsX-1; i++ ) {     
       for ( j=1; j<nPointsY-1; j++ ) {
        k=j*nPointsX +i;

	drift_r = -0.5*rate_of_meanreversion_r*(mean_Reverting_r - y[k])*(y[k+1] - y[k-1])*dx_inv;
        drift_l = -0.5*rate_of_meanreversion_l*(mean_Reverting_l - y[k])*(y[k+nPointsX] - y[k-nPointsX])*dl_inv;

        sigma_r_tmp=sigma_r(t,x_arr[i]);
        sigma_l_tmp=sigma_l(t,l_arr[j]);
        result[k] = drift_r + drift_l 
           +0.5e0*sigma_r_tmp*sigma_r_tmp*(y[k+1]-2.0e0*y[k]+ y[k-1])*dx_sq_inv 
           +0.5e0*sigma_l_tmp*sigma_l_tmp*(y[k+nPointsX]-2.0e0*y[k]+ y[k-nPointsX])*dl_sq_inv
           +rho*sigma_r_tmp*sigma_l_tmp*dx_dl_inv*((y[k+nPointsX+1]-y[k+nPointsX-1])-(y[k-nPointsX+1]-y[k-nPointsX-1]))
           -(x_arr[i]+l_arr[j])*y[k]
           +l_arr[j]*expected_default_payoff(t,x_arr[i],l_arr[j])
           +continuouse_payoff(t,x_arr[i],l_arr[j]);
       }
     }

     //this part is for output only.
     //If you want efficiency-comment it out
     if ( counter_call_back%25 == 0 && t_old < t ) {
        t_old=t;
        for ( int k=0; k < neqn; k++ ) {
           int j=k/nPointsX;
           int i=k-j*nPointsX;
           fprintf( solution, "%f %f %f %f\n", t, x_arr[i], l_arr[j], y[k] );
        };
     };
     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/(nPointsX*nPointsY)); //this number is just from experiment
                                    //one can use Gershgorin's theoreme to estimete max_eigenvalue
   };

   double y[1000000],z0[1000000],z1[1000000],z2[1000000],z3[1000000], oldEigenVector[1000000];
   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);

      int isComputeEigenValue=(1==0);
      long neqn=nPointsX*nPointsY;      
      double t=0.0e0, tend=1.0e0, h0=0.00001e0, atol=0.001e0, rtol=0.001e0;

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


      //----- initial values ----//
      x_arr=new double[nPointsX];
      l_arr=new double[nPointsY];



      for ( int k=0; k < neqn; k++ ) {
         int j=k/nPointsX;
         int i=k-j*nPointsX;
         x_arr[i]=x_min+dx*((double)i);
         l_arr[j]=l_min+dl*((double)j);
         y[k]=final_payoff(1.0e0, x_arr[i], l_arr[j]);
         fprintf( payoffFunction, "%f %f %f\n", x_arr[i], l_arr[j], y[k] );
      };    

      // start computation
      dumka3(neqn, &t, tend, h0,
         atol, rtol, rhs, cour,
         y, z0, z1, z2, z3, oldEigenVector, isComputeEigenValue);
      
      fclose( solution );
      fclose( solution_exact );
      fclose( payoffFunction );
      return 0;
   };//end of main
   // remember - you need dumka3.c to run this