/*
*			BlackScholesExp.c:
*
*		Sample code for AMS-691 homework.
*
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define		sqr(x) 		((x)*(x))
#define		min(x,y)	(((x) < (y)) ? (x) : (y))

static void expiry_func(double,double*,double*,int);

main(int argc, char **argv)
{ 
	double *x,*u_old,*u_new;
	double time,dx,dt,dt_cfl;
	int i,n;
	FILE *xg_file;
	double E = 1.0;
	double gamma = 0.08;
	double sigma = 0.2;
	double CFL = 0.7;
	double print_time;
	int is_print_time;
	char fname[100];

	/* computational domain and mesh parameters */
	double L = 0.0;
	double U = 10.0;
	int  mesh_size = 401;

	x = (double*)malloc(mesh_size*sizeof(double));
	u_old = (double*)malloc(mesh_size*sizeof(double));
	u_new = (double*)malloc(mesh_size*sizeof(double));

	/* Initialization of states */

	dx = (U - L)/(mesh_size - 1);

	dt_cfl = dx/gamma/10.0; // for hyperbolic part
	dt_cfl = min(dt_cfl,sqr(dx)/0.5/sqr(sigma*10.0))/2.0; 
				// for parabolic part
	dt_cfl *= CFL;

	for (i = 0; i < mesh_size; i++)
	    x[i] = L + i*dx;

	/* Set the initial condition */

	expiry_func(E,x,u_old,mesh_size);
	expiry_func(E,x,u_new,mesh_size);


	/* Time loop */

	n = 0;   time = 0.0;

        sprintf(fname,"bs_data-%3.1f",time);
        xg_file = fopen(fname,"w");
        fprintf(xg_file,"\"Option price vs. stock value\"\n");
        for (i = 0; i < mesh_size/4; i++)
            fprintf(xg_file,"%f  %f\n",x[i],u_old[i]);
        fclose(xg_file);

	print_time = 1.0;
	while (time <= 4.0)
	{
	    dt = dt_cfl;
	    is_print_time = 0;
	    if (time < print_time && time+dt >= print_time)
	    {
		dt = print_time - time;
		print_time += 1.0;
		is_print_time = 1;
	    }
	    for (i = 1; i < mesh_size-1; ++i)
	    {
		double a = -gamma*x[i];
		double b = 0.5*sqr(sigma*x[i]);
		u_new[i] = u_old[i] - a*(u_old[i+1] - u_old[i])*dt/dx
			+ b*(u_old[i+1] - 2.0*u_old[i] + u_old[i-1])*dt/sqr(dx)
			- dt*gamma*u_old[i];
	    }
	    n++;
	    time += dt;
	    for (i = 1; i < mesh_size-1; i++)
	    {
		u_old[i] = u_new[i];
	    }
	    u_old[0] = 0.0;
	    u_old[mesh_size-1] = u_old[mesh_size-2] + dx;
	    if (is_print_time == 1)
	    {
		sprintf(fname,"bs_data-%3.1f",time);
		xg_file = fopen(fname,"w");
		fprintf(xg_file,"\"Option price vs. stock value\"\n");
	    	for (i = 0; i < mesh_size/4; i++)
		    fprintf(xg_file,"%f  %f\n",x[i],u_old[i]);
		fclose(xg_file);
	    }
	}
	printf("n = %d\n",n);
}

static void expiry_func(
	double E,
	double *x,
	double *u,
	int mesh_size)
{
	int i;
	double arg;

	for (i = 0; i < mesh_size; ++i)
	{
	    arg = x[i];
	    if (arg > E) u[i] = arg - E;
	    else u[i] = 0.0;
	}
}