// the fundamental diagram of NS.cpp 
#include "stdafx.h"
#include<iostream>
#include<stdlib.h>
#include <time.h>
using namespace std;
int	x[2000],i,v[2000],gap[2000],V,x1[2000],v1[2000],a[2000],b[2000],sum;
float p2;
int	rule(int l,int m);
int	Max(int	a,int b);
int	Min(int	d,int e);
int	_tmain(int argc, _TCHAR* argv[])
{
	int	L,j,flg,time_steps,t,V1,N,h,c[2000],q;
	double p1=0;
	float V2,V3=0,V4,V_total,V_average,J;
	cout<<"Please enter the number of sites L:			"<<endl;
	cin>>L;
	cout<<"Please enter the maximum of velocity V:			"<<endl;
	cin>>V;
	cout<<"Please enter the stochastic braking probability p2:			"<<endl;
	cin>>p2;
	cout<<"Please enter the time_steps:			 "<<endl;
	cin>>time_steps;
	cout<<"Please enter t:			"<<endl;
	cin>>t;
	cout<<"Please enter the number of samples N:		  "<<endl;
	cin>>N;
	srand( (unsigned)time( NULL	) );
	while(p1<=1.01)
	{
		cout<<"p1="<<p1<<"  ";
		V_total=0;
		for(h=1;h<=N;h++)
		{
			for(i=0;i<L;i++)
			{
				c[i]=0;		
			}
			for(i=0;i<int(L*p1);i++)
			{
				q=int(rand()*(L-1)/RAND_MAX);
				if(c[q]==0)
					c[q]=1;
				else
					i=i-1;
			}
			sum=0;
			for(i=0;i<L;i++)
			{
				if(c[i]==1)
				{
					sum=sum+1;
					x[sum]=i;
				    v[sum]=0;
				}
				
			}	
			V3=0;
			for(int	k=1;k<=time_steps;k++)
			{


				for(i=1;i<=sum;i++)
				{
					a[i]=x[i];
					b[i]=v[i];
				}
				for(i=1;i<=sum;i++)
				{
					if(i==sum)
						gap[i]=a[1]-a[i]-1;
					else
						gap[i]=a[i+1]-a[i]-1;
					if(gap[i]<0)
						gap[i]=gap[i]+L;
				}
				for(i=1;i<=sum;i++)
				{
					rule(b[i],gap[i]);
				}
				for(i=1;i<=sum;i++)
				{
					if((x[i]+v[i])>L)
						x[i]=x[i]+v[i]-L;
					else 
						x[i]=x[i]+v[i];
				}

				if(k>=t)
				{
					V1=0;
					for(i=1;i<=sum;i++)
						V1=V1+v[i];
					V2=(float)(V1)/(float)(sum);
					V3=V3+V2;
				}
			}

			V4=V3/(float)(time_steps-t+1);
			V_total=V_total+V4;
		}
		V_average=V_total/(float)(N);
		cout<<"V_average="<<V_average<<"  ";
		J=V_average*p1;
		cout<<"J="<<J<<endl;
		FILE *fp;
		if((fp=fopen("E:\\file.dat","a+"))==NULL)
		{
			cout<<"Can not open this file."<<endl;
			exit(0);
		}
		fprintf(fp,"%f,%f,%f\n",p1,V_average,J);
		fclose(fp);

		p1=p1+0.01;
	}
	return 0;		
}
int	rule(int l,int m)
{
	//if(l<V&&m>(l+1))
	//l=l+1;
	l=Min((l+1),V);
	//if((m+1)<l)
	//l=gap[i];	
	l=Min(l,m);
	if((float)(rand())/(float)(RAND_MAX)<p2)
		l=Max((l-1),0);
	v[i]=l;
	return v[i];
}
int	Max(int	a,int b)
{
	int	c;
	c=a>b?a:b;
	return c;
}
int	Min(int	d,int e)
{
	int	f;
	f=d<e?d:e;
	return f;
}