// the fundamental diagram of NS model with open conditions.cpp 
#include "stdafx.h"
#include"stdafx.h"
#include<iostream>
#include<stdlib.h>
#include <time.h>
using namespace std;
int x[2000],i,v[2000],gap[2000],V=4,x1[2000],v1[2000],sum=0;
double p2=0,g;
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,c[2000],q,h,N,V1;
	double p1,p_injection,p_extinction,p_each_step,p_steps,p_each_step_average,p_total=0,p_average;
	float V2,V3=0,V4,V_total=0,V_average,J;
	cout<<"Please enter the number of sites L:          "<<endl;
	cin>>L;
	cout<<"Please enter the density of cars per site p1:          "<<endl;
	cin>>p1;
	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 injection rate p_injection:          "<<endl;
	//	cin>>p_injection;
	cout<<"Please enter the extinction rate p_extinction:          "<<endl;
	cin>>p_extinction;
	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(V<=10)
	while(p2<=0.88)
	{
		cout<<"P="<<p2<<endl;
		p_injection=0.05;
		while(p_injection<=1.01)
			//		p_extinction=0.05;
			//	while(p_extinction<=1.01)

		{
			cout<<"p_injection="<<p_injection<<"  ";
			//	cout<<"p_extinction="<<p_extinction<<"  ";
			p_total=0;
			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((float)rand()/(float)RAND_MAX*(L-1)+1);     //隨機產生一個從1到L的數；
					if(c[q]==0)
						c[q]=1;
					else
						i=i-1;
				}
				sum=0;
				if((float)rand()/(float)RAND_MAX<=p_injection)
				{
					c[0]=1;
					sum=sum+1;
					v[sum]=V;
					x[sum]=0;           
				}                      //在格點i=0處以p_injection概率隨機注入一輛v=V的車；
				for(i=1;i<=L;i++)
				{
					if(c[i]==1)
					{
						sum=sum+1;
						x[sum]=i;
						v[sum]=0;
						//cout<<v[sum];
					}
					//else
					//cout<<".";
				}	
				/*	if(c[0]==1)
				cout<<"sum="<<sum-1;
				else
				cout<<"sum="<<sum;
				cout<<endl;*/

				//輸出t=0時步的構形；
				for(i=1;i<=sum;i++)
				{
					x1[i]=x[i];
					v1[i]=v[i];
				}
				p_steps=0;
				V3=0;
				for(int k=1;k<time_steps;k++)
				{		

					for(i=1;i<=sum;i++)
					{
						if(i==sum)
						{
							if((float)rand()/(float)RAND_MAX<=(1-p_extinction))//在格點i=L+1處以1-p_extinction的概率產生一輛堵塞車;
								gap[i]=L+1-x1[i]-1;
							else
								gap[i]=V;
						}
						else
							gap[i]=x1[i+1]-x1[i]-1;
					}
					for(i=1;i<=sum;i++)
					{
						rule(v1[i],gap[i]);
					}
					for(i=sum;i>=1;i--)
					{
						x[i]=x[i]+v[i];
						if(x[i]>L)
							sum=sum-1;           // 以p_extinction的概率消失一輛車；
					}
					if(x[1]==0)
					{
						for(i=1;i<sum;i++)
						{
							x[i]=x[i+1];
							v[i]=v[i+1];
						}
						sum=sum-1;             //如果入射車在下一時步的速度為零，那么這輛車被刪除；
					}
					for(i=1;i<=L;i++)
					{
						flg=0;
						for(j=1;j<=sum;j++)
						{
							if(i==x[j])
							{
								c[i]=1;
								flg=1;
							}				
						}
						if(flg==0)
						{
							c[i]=0;
						}
					}

					for(i=1;i<=sum;i++)
					{
						x1[i]=x[i];
						v1[i]=v[i];
					}
					sum=0;
					c[0]=0;
					if((float)rand()/(float)RAND_MAX<=p_injection)
					{
						c[0]=1;
						sum=sum+1;
						v[sum]=V;
						x[sum]=0;         //在格點i=0處以p_injection概率隨機注入一輛v=V的車；  
						for(i=1;i<=L;i++)
						{
							if(c[i]==1)
							{
								sum=sum+1;
								v[sum]=v1[sum-1];
								x[sum]=x1[sum-1];
								//cout<<v[sum];
							}
							//else
							//cout<<".";
						}
						/*	if(c[0]==1)
						cout<<"sum="<<sum-1;
						else
						cout<<"sum="<<sum;
						cout<<endl;*/


					}
					else
					{
						for(i=1;i<=L;i++)
						{
							if(c[i]==1)
							{
								sum=sum+1;
								v[sum]=v1[sum];	
								x[sum]=x1[sum];
								//cout<<v[sum];
							}	
							//else
							//cout<<".";
						}
						/*if(c[0]==1)
						cout<<"sum="<<sum-1;
						else
						cout<<"sum="<<sum;
						cout<<endl;*/

					}
					if(k>=(t-1))
					{
						if(c[0]==1)
							p_each_step=(float)(sum-1)/(float)L;
						else
							p_each_step=(float)sum/(float)L;
						//cout<<"p_each_step="<<p_each_step;
						p_steps=p_steps+p_each_step;
						V1=0;
						for(i=1;i<=sum;i++)
							V1=V1+v[i];
						if(c[0]==1)
						{
							V1=V1-v[1];
							V2=(float)(V1)/(float)(sum-1);
						}
						else
							V2=(float)(V1)/(float)(sum);
						//cout<<"V2="<<V2<<endl;
						V3=V3+V2;
					}

					for(i=1;i<=sum;i++)
					{
						x1[i]=x[i];
						v1[i]=v[i];
					}
				}
				p_each_step_average=p_steps/(float)(time_steps-t+1);
				p_total=p_total+p_each_step_average;
				//cout<<"p_total="<<p_total<<"  "<<endl;
				V4=V3/(float)(time_steps-t+1);
				V_total=V_total+V4;
				//cout<<"V_total="<<V_total<<endl;
			}
			p_average=p_total/(float)(N);
			cout<<"p_average="<<p_average<<"  ";
			V_average=V_total/(float)(N);
			cout<<"V_average="<<V_average<<"  ";
			J=V_average*p_average;
			cout<<"J="<<J<<endl;
			FILE *fp;
			if((fp=fopen("E:\\開放邊條下的NS模型[文獻]V_Max=5,（β=1，P不同）.dat","a+"))==NULL)			
			{
				cout<<"Can not open this file."<<endl;
				exit(0);
			}
			fprintf(fp,"%f  %f  %f  %f  %f\n",p2,p_injection,p_average,V_average,J);
			fclose(fp);

			p_injection=p_injection+0.05;
			//		p_extinction=p_extinction+0.05;
		}
		//	V=V+1;
		p2=p2+0.125;
	}
	return 0;
}
int rule(int l,int m)
{
	l=Min((l+1),V);
	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;
}