//the mixed traffic flow (fast and slow) with open boundary conditions:
#include "stdafx.h"
#include<iostream>
#include<stdlib.h>
#include <time.h>
#using <mscorlib.dll>
using namespace std;
using namespace System;
int x[2000],i,v[2000],gap[2000],V_slow_max,V_fast_max,x1[2000],v1[2000],sum=0;
float p=0;
int rule(int l,int m,int V_max,float p_noise);
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],c_temporary[2000],q,h,N;
	double f,p1,p_injection,p_extinction,p_injection_NS,p_each_step,p_steps,p_each_step_average,p_total=0,p_average;
	float R,V_EachStep_Sum,V_EachStep_Ave,V_Steps_Sum=0,V_Steps_Ave,V_Samples_Sum=0,V_Samples_Ave,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 of slow vehicle V_slow_max:          "<<endl;
	cin>>V_slow_max;
	cout<<"Please enter the maximum of velocity of fast vehicle V_fast_max:          "<<endl;
	cin>>V_fast_max;
	//	cout<<"Please enter the stochastic braking probability p:          "<<endl;
	//	cin>>p;
	//	cout<<"Please enter the injection rate p_injection:          "<<endl;
	//	cin>>p_injection;
	cout<<"Please enter the injection rate of NS p_injection_NS:          "<<endl;
	cin>>p_injection_NS;
	cout<<"Please enter the extinction rate p_extinction:          "<<endl;
	cin>>p_extinction;
	cout<<"Please enter f:          "<<endl;
	cin>>f;
	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(p<=0.88)
	{
		cout<<"p="<<p<<endl;
		p_injection=0.05;
		while(p_injection<=1.01)
		{
			cout<<"p_injection="<<p_injection<<"  ";			
			p_total=0;
			V_Samples_Sum=0;

			for(h=1;h<=N;h++)
			{

				for(i=0;i<=L;i++)//產(chǎn)生初始構(gòu)型
				{
					c[i]=0;		
				}
				for(i=0;i<Math::Round(L*p1*f);i++)
				{
					q=int((float)rand()/(float)RAND_MAX*(L-1)+1);
					if(c[q]==0)
						c[q]=3;
					else
						i=i-1;
				}
				for(i=0;i<Math::Round(L*p1*(1-f));i++)
				{
					q=int((float)rand()/(float)RAND_MAX*(L-1)+1);
					if(c[q]==0)
						c[q]=5;
					else
						i=i-1;
				}
				sum=0;
				R=(float)rand()/(float)RAND_MAX;
				if(R<=p_injection)//在格點i=0處以p_injection概率隨機(jī)注入一輛v=V的車；
				{
					if(R<=p_injection_NS*p_injection)
					{
						c[0]=3;
						sum=sum+1;
						v[sum]=V_slow_max;
						x[sum]=0;
					}
					else
					{
						c[0]=5;
						sum=sum+1;
						v[sum]=V_fast_max;
						x[sum]=0;           
					}  
				}

				for(i=1;i<=L;i++)
				{
					if(c[i]!=0)
					{
						sum=sum+1;
						x[sum]=i;
						v[sum]=0;
						//	cout<<v[sum];
						//cout<<c[i];
					}
					//	else
					//	cout<<".";
				}	
				/*	if(c[0]!=0)
				cout<<"sum="<<sum-1;
				else
				cout<<"sum="<<sum;
				cout<<endl;*/
				//輸出t=0時步的構(gòu)形；
				/*		for(i=1;i<=sum;i++)
				{
				x1[i]=x[i];
				v1[i]=v[i];
				}*/
				p_steps=0;
				V_Steps_Sum=0;
				for(int k=1;k<time_steps;k++)
				{	
					for(i=1;i<=sum;i++)
					{
						x1[i]=x[i];
						v1[i]=v[i];
					}

					for(i=1;i<=sum;i++)
					{
						if(i==sum)
						{
							if((float)rand()/(float)RAND_MAX<=(1-p_extinction))//在格點i=L+1處以1-p_extinction的概率產(chǎn)生一輛堵塞車;
								gap[i]=L+1-x1[i]-1;
							else
								gap[i]=V_fast_max;
						}
						else
							gap[i]=x1[i+1]-x1[i]-1;
					}

					for(i=0;i<=L;i++)
					{
						c_temporary[i]=c[i];
						c[i]=0;
					}
					for(i=sum;i>=1;i--)
					{
						if(c_temporary[x[i]]==3)
						{
							rule(v1[i],gap[i],V_slow_max,p);
							x[i]=x[i]+v[i];
							if(x[i]>L)
								sum=sum-1; 
							c[x[i]]=3;    
						}
						else
						{
							rule(v1[i],gap[i],V_fast_max,p);
							x[i]=x[i]+v[i];
							if(x[i]>L)
								sum=sum-1; 
							c[x[i]]=5; 
						}
					}
					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<=sum;i++)
					{
						x1[i]=x[i];
						v1[i]=v[i];
					}
					sum=0;
					c[0]=0;  
					R=(float)rand()/(float)RAND_MAX;
					if(R<=p_injection)//在格點i=0處以p_injection概率隨機(jī)注入一輛v=V的車；
					{
						if(R<=p_injection_NS*p_injection)
						{
							c[0]=3;
							sum=sum+1;
							v[sum]=V_slow_max;				
							x[sum]=0; 
						}
						else
						{
							c[0]=5;
							sum=sum+1;
							v[sum]=V_fast_max;				
							x[sum]=0; 
						}//在格點i=0處以p_injection概率隨機(jī)注入一輛v=V的車；  
						for(i=1;i<=L;i++)
						{
							if(c[i]!=0)
							{
								sum=sum+1;
								v[sum]=v1[sum-1];
								x[sum]=x1[sum-1];
								//	cout<<v[sum];
								//cout<<c[i];
							}
							//	else
							//		cout<<".";
						}
						/*	if(c[0]!=0)
						cout<<"sum="<<sum-1;
						else
						cout<<"sum="<<sum;
						cout<<endl;*/

					}
					else
					{
						for(i=1;i<=L;i++)
						{
							if(c[i]!=0)
							{
								sum=sum+1;
								v[sum]=v1[sum];	
								x[sum]=x1[sum];
								//	cout<<v[sum];
								//cout<<c[i];
							}	
							//	else
							//		cout<<".";
						}
						/*	if(c[0]!=0)
						cout<<"sum="<<sum-1;
						else
						cout<<"sum="<<sum;
						cout<<endl;*/
					}
					if(k>=(t-1))
					{
						if(c[0]!=0)
							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;
						V_EachStep_Sum=0;
						for(i=1;i<=sum;i++)
							V_EachStep_Sum=V_EachStep_Sum+v[i];
						if(c[0]!=0)
						{
							V_EachStep_Sum=V_EachStep_Sum-v[1];
							V_EachStep_Ave=(float)(V_EachStep_Sum)/(float)(sum-1);
						}
						else
							V_EachStep_Ave=(float)(V_EachStep_Sum)/(float)(sum);
						//cout<<"V_EachStep_Ave="<<V_EachStep_Ave<<endl;
						V_Steps_Sum=V_Steps_Sum+V_EachStep_Ave;
					}

					/*			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;
				V_Steps_Ave=V_Steps_Sum/(float)(time_steps-t+1);
				V_Samples_Sum=V_Samples_Sum+V_Steps_Ave;
				//cout<<"V_Samples_Sum="<<V_Samples_Sum<<endl;
			}
			p_average=p_total/(float)(N);
			cout<<"p_average="<<p_average<<"  ";
			V_Samples_Ave=V_Samples_Sum/(float)(N);
			cout<<"V_Samples_Ave="<<V_Samples_Ave<<"  ";
			J=V_Samples_Ave*p_average;
			cout<<"J="<<J<<endl;
			FILE *fp;
			if((fp=fopen("E:\\the mixed traffic flow (fast and slow)with open boundary conditions.dat","a+"))==NULL)

			{
				cout<<"Can not open this file."<<endl;
				exit(0);
			}
			fprintf(fp,"%f  %f  %f  %f  %f\n",p,p_injection,p_average,V_Samples_Ave,J);
			fclose(fp);

			p_injection=p_injection+0.05;				
		}
		p=p+0.125;
	}
	return 0;
}
int rule(int l,int m,int V_max,float p_noise)//p_noise剎車概率
{	
	l=Min((l+1),V_max);	
	l=Min(l,m);
	if((float)(rand())/(float)(RAND_MAX)<p_noise)
		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;
}