//the fundamental diagram of the mixed traffic flow ( NS and FI):
#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_NS_max,V_FI_max,x1[2000],v1[2000],a[2000],b[2000],sum=0;
float p_NS,p_FI;
int rule_NS(int l,int m,int V_max,float p_noise);
int rule_FI(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 L,j,flg,time_steps,c[2000],c_temporary[2000],q,t,N,h,test;
	double p1,f=0,V_EachStep_Sum,V_EachStep_Ave,V_Steps_Sum=0,V_Steps_Ave,V_Samples_Sum,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 NS vehicle V_NS_max:          "<<endl;
	cin>>V_NS_max;
	cout<<"Please enter the maximum of velocity of FI vehicle V_FI_max:          "<<endl;
	cin>>V_FI_max;
	cout<<"Please enter the stochastic braking probability of slow p_NS:          "<<endl;
	cin>>p_NS;
	cout<<"Please enter the stochastic braking probability of fast p_FI:          "<<endl;
	cin>>p_FI;
//	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_FI<=1.01)
	//	while(p_NS<=1.01)
	while(f<=1.01)
	{
		//		cout<<"p_FI="<<p_FI<<endl;
		//		cout<<"p_NS="<<p_NS<<endl;
		cout<<"f="<<f<<endl;
		p1=0.01;
		while(p1<=1.01)
		{
			cout<<"p1="<<p1<<"  ";
			V_Samples_Sum=0;
			for(h=1;h<=N;h++)
			{
				for(i=1;i<=L;i++)//產生初始構型
				{
					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;
				for(i=1;i<=L;i++)//輸出初始構型
				{
					if(c[i]!=0)
					{
						sum=sum+1;
						x[sum]=i;
						v[sum]=0;
						///		cout<<v[sum];
					}
					//	else
					//		cout<<".";
				}	
				//		cout<<endl;	
				V_Steps_Sum=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<=L;i++)
					{
						c_temporary[i]=c[i];
						c[i]=0;
					}
					for(i=1;i<=sum;i++)
					{
						if(c_temporary[x[i]]==3)
						{
							rule_NS(b[i],gap[i],V_NS_max,p_NS);
							if((x[i]+v[i])>L)			
								x[i]=x[i]+v[i]-L; 
							else 
								x[i]=x[i]+v[i];
							c[x[i]]=3;    
						}
						else
						{
							rule_FI(b[i],gap[i],V_FI_max,p_FI);
							if((x[i]+v[i])>L)			
								x[i]=x[i]+v[i]-L;               
							else 
								x[i]=x[i]+v[i];
							c[x[i]]=5; 
						}
					}		
					/*			for(j=1;j<=L;j++)
					{
					flg=0;
					for(i=1;i<=sum;i++)
					{
					if(j==x[i])
					{					
					cout<<v[i];
					flg=1;
					}
					}
					if(flg==0)
					cout<<".";

					}
					cout<<endl;*/
					if(k>=(t-1))
					{
						V_EachStep_Sum=0;
						for(i=1;i<=sum;i++)
							V_EachStep_Sum=V_EachStep_Sum+v[i];
						V_EachStep_Ave=V_EachStep_Sum/(L*p1);			
						V_Steps_Sum=V_Steps_Sum+V_EachStep_Ave;
					}


				}
				V_Steps_Ave=V_Steps_Sum/(float)(time_steps-t+1);		
				V_Samples_Sum=V_Samples_Sum+V_Steps_Ave;
			}
			V_Samples_Ave=V_Samples_Sum/N;
			cout<<"V_Samples_Ave="<<V_Samples_Ave<<"  ";
			J=V_Samples_Ave*p1;
			cout<<"J="<<J<<endl;
			FILE *fp;
			if((fp=fopen("E:\\f對基本圖的影響（P_NS=0.5,P_FI=0.3）.dat","a+"))==NULL)			
			{
				cout<<"Can not open this file."<<endl;
				exit(0);
			}
			fprintf(fp,"%f,%f,%f,%f\n",f,p1,V_Samples_Ave,J);
			fclose(fp);
			p1=p1+0.01;
		}		
		//	p_FI=p_FI+0.1;
		//	p_NS=p_NS+0.1;
		f=f+0.2;

	}
	return 0;
}
int rule_NS(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 rule_FI(int l,int m,int V_max,float p_noise)
{
	if(m>=V_max)
	{
		l=V_max;
		if((float)(rand())/(float)(RAND_MAX)<p_noise)
			l=Max((l-1),0);
	}
	else
		l=m;
	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;
}