#include <stdio.h>
#include <math.h>
#include "numcip.h"

/*constants for ran11(long *)*/
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

long seed;//seed for random generator

//Linear SDE parameters
double xinit;
double lambda,mu;
double combp;

double ran11(long *idum);
double gasdev(long *idum);

double winc(double *,int,int);
double exsoln(double tt,double ww);

main()
{
	int i,N;
	int k,M;
	double dt,sdt;
	double t,w;
	double *dw;

	double x;
	double Dt,Dw;
	int L;
	int R;

	double *sum;
	char filename[25];

	FILE *outfile;

	// user io
	printf("Please input random seed [-xxx]:\n");
	scanf("%ld",&seed);
	N=1024;
	xinit=1.0;
	lambda=2.0;mu=1.0;
	combp=lambda-0.5*mu*mu;
	R=4;

	dt=(double)(1.0/N);
	sdt=sqrt(dt);//this is the most important part in modeling stochastic process

	// vector storage for Brownian increments
	dw=dvector(1,N);

	L=N/R;
	Dt=dt*R;

	// for calculate the mean value
	sum=dvector(1,L);
	M=10000;


	for(k=1;k<=3;k++)//we will plot the first 3 sample of the solution
	{
	// Generate Brownian Path and the exact solution to the linear SDE
		sprintf(filename,"sde_theo%d.dat",k);
		outfile=fopen(filename,"w");
		fprintf(outfile,"tt\ttheo\n");
		fprintf(outfile,"%lg\t%lg\n",0.0,xinit);

		t=0.0;
		w=0.0;
		for(i=1;i<=N;i++)
			{
			dw[i]=sdt*(double)gasdev(&seed);
			t+=dt;
			w+=dw[i];
			fprintf(outfile,"%lg\t%lg\n",t,exsoln(t,w));
			}
		fclose(outfile);

		//Do the SDE with Euler-Maruyama method
		sprintf(filename,"sde_em%d.dat",k);
		outfile=fopen(filename,"w");
		fprintf(outfile,"tt\txx\n");
		fprintf(outfile,"%lg\t%lg\n",0.0,xinit);
		

		x=xinit;
		t=0.0;
		for(i=1;i<=L;i++)
			{
			Dw=winc(dw,R*(i-1)+1,R*i);
			t+=Dt;
			x+=(lambda*x*Dt+mu*x*Dw);

			sum[i]+=x;//summing to calculate the mean
			fprintf(outfile,"%lg\t%lg\n",t,x);
			}
		fclose(outfile);
	}

	//calculating other realizations of the random process for mean value 
	for(k=4;k<=M;k++)
	{
	// Generate Brownian Path 
		for(i=1;i<=N;i++)
			{
			dw[i]=sdt*(double)gasdev(&seed);
			}

		//Do the SDE with Euler-Maruyama method
		x=xinit;
		t=0.0;
		for(i=1;i<=L;i++)
			{
			Dw=winc(dw,R*(i-1)+1,R*i);
			t+=Dt;
			x+=(lambda*x*Dt+mu*x*Dw);
			
			sum[i]+=x;//summing to calculate the mean
			}
	}

	//calculate the mean of the stochastic process
	outfile=fopen("sde_mean.dat","w");
	fprintf(outfile,"tt\txx\n");
	fprintf(outfile,"%lg\t%lg\n",0.0,xinit);


	for(i=1;i<=L;i++)
		{
		sum[i]/=M;
		fprintf(outfile,"%lg\t%lg\n",i*Dt,sum[i]);
		}
	fclose(outfile);
	
	return 0;
}

double winc(double *ddw,int start,int stop)
//this calculated the increment in the Brownian path
{
	int i;
	double temp;

	temp=0.0;
	for(i=start;i<=stop;i++)
		temp+=ddw[i];

	return temp;
}

double exsoln(double tt,double ww)
//this calculates the exact solution to the linear SDE
{
	return xinit*exp(combp*tt+mu*ww);
}

//
// random number generators
//

double ran11(long *idum)
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	double temp;
	
	if(*idum <= 0 || !iy)
		{
		if(-(*idum) < 1)
			*idum = 1;
		else
			*idum = -(*idum);
		for(j=NTAB+7;j>=0;j--)
			{
			k=(*idum)/IQ;
			*idum=IA*(*idum-k*IQ)-IR*k;
			if(*idum<0)
				*idum += IM;
			if(j<NTAB)
				iv[j]= *idum;
			}
		iy=iv[0];
		}
	k=(*idum)/IQ;
	*idum=IA*(*idum-k*IQ)-IR*k;
	if(*idum < 0)
		*idum += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j]= *idum;
	if((temp=AM*iy) >RNMX)
		return RNMX;
	else
		return temp;
}

double gasdev(long *idum)
{
	static int iset=0;
	static double gset;
	double fac,r,v1,v2;
    
	if  (iset == 0) {
		do {
			v1=2.0*ran11(idum)-1.0;
			v2=2.0*ran11(idum)-1.0;
			r=v1*v1+v2*v2;
		} while (r >= 1.0);
		fac=sqrt(-2.0*log(r)/r);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} else {
		iset=0;
		return gset;
	}
}