nsf.cpp

用CA模拟交通流基本模型NaSch模型,能再现交通流的基本自由相和阻塞相

#include<stdlib.h>
#include<iostream.h>
#include<time.h>
#include<math.h>
#include<fstream.h>
#include<iomanip.h>
#define L 300
struct car
{
	int max_v;
	int x;
	int v;
};
struct car now[2000],old[2000],temp[2000];
void evolve(int, struct car *,struct car *);
void print(int, struct car*);
void sort(struct car*,struct car*); 
int max(int x, int y)
{
	int z;
	z=x>y?x:y;
	return z;
}
int min(int x, int y)
{
	int z;
	z=x<y?x:y;
	return z;
} 
void main()
{
	 srand((unsigned)time(NULL));
	 int i,j,d;
	 double density,aver,density1,aver1,flux;
	 ofstream file1("d:\\Nsflow'.txt");
   for(d=1;d<=100;d++)
   {
	  for(i=1;i<=L;i++)
	   {
		   old[i].x=now[i].x=temp[i].x=0;
		   old[i].max_v=now[i].max_v=temp[i].max_v=0;
		   old[i].v=now[i].v=temp[i].v=0;
	   }
	  int t,m,tmp[2000],n,count=0;
	  for(i=1;i<L;i++)
		 tmp[i]=L+2;
    	for(i=1;i<=L;i++)
		{
	    	m=i;	
		   t=L*rand()/RAND_MAX+1;
		   for(j=1;j<=L;j++)
		   if(t==tmp[j])
		   {i=i-1;break;}
		   if(i!=m)
		     continue;
	    	tmp[i]=t;
		    count++;		   
			old[t].max_v=5;
			old[t].v=5;
			old[t].x=t;	
			if(count>=L*d/100)
		    	break;
		}
	    sort(old,temp);
		 aver=0,density=0,flux=0;
	   for(int time_step=0;time_step<15000;time_step++)
	   {
		   n=0;aver1=0;
	     for(i=1,j=1;i<=L;i++)
		   if(old[i].x!=0)
		   {		
			now[j].max_v=old[i].max_v;
			now[j].v=old[i].v;
			now[j].x=old[i].x;
			old[i].max_v=0;
			old[i].v=0;
			old[i].x=0;
            j++;
		    n++;
		   }
	     now[n+1].max_v=now[1].max_v;
		 now[n+1].v=now[1].v;
	     now[n+1].x=now[1].x+L;
	     cout<<setw(6)<<n<<setw(6)<<count<<endl;
         evolve( n,now,old);
		 if(time_step>=10000)
		 {
			 for(int i=1;i<=L;i++)
				 aver1+=(double)old[i].v/n;
			  // print(L,old);
			 density1=(double)n/L;
	         flux+=aver1*density1/5000;
			 aver+=aver1/5000;
			 density+=density1/5000;
		 }
	   }
	   //density=(double)n/L;
	   //flux=(double)aver*density;
	  file1<<density<<" "<<aver<<" "<<flux<<endl;
	  }
     file1.close();	    
}
 void evolve(int m,struct car *pt,struct car *pt1)     /////(n,now,old)
{
     int i,T;
	 float p;
	 //p=(float)rand()/RAND_MAX;
	 for( i=1;i<=m;i++)
	 {
		  p=(float)rand()/RAND_MAX;
	   pt[i].v=min(pt[i].v+1,pt[i].max_v);
	   pt[i].v=min(pt[i].v,pt[i+1].x-pt[i].x-1);
	   if(p<0.5)
        pt[i].v=max(pt[i].v-1,0);
	   pt[i].x=pt[i].x+pt[i].v;
	 }
	for( i=1;i<=m;i++)
	{
		if((pt+i)->x>L)
			(pt+i)->x=(pt+i)->x-L;
			T=(pt+i)->x;
	 (pt1+T)->x=(pt+i)->x; 
	 (pt1+T)->v=(pt+i)->v;  
	 (pt1+T)->max_v=(pt+i)->max_v; 
	 (pt+i)->x=0;
	 (pt+i)->v=0;
	 (pt+i)->max_v=0;
	}
	 (pt+m+1)->x=0;
	 (pt+m+1)->v=0;
	 (pt+m+1)->max_v=0;
}
 void sort(struct car*ptr,struct car*ptr1)
 {
	for(int i=1;i<=L;i++)
	{
	if((ptr+i)->x>L)
		(ptr+i)->x=(ptr+i)->x-L;
	 (ptr1+(ptr+i)->x)->x=(ptr+i)->x; 
	 (ptr1+(ptr+i)->x)->v=(ptr+i)->v;  
	 (ptr1+(ptr+i)->x)->max_v=(ptr+i)->max_v;
	}
	for(i=1;i<=L;i++)
	{
	 (ptr+i)->x=(ptr1+i)->x;
	 (ptr+i)->v=(ptr1+i)->v;
	 (ptr+i)->max_v=(ptr1+i)->max_v;
	 (ptr1+i)->x=0;
	 (ptr1+i)->v=0;
	 (ptr1+i)->max_v=0;
	}
 }
 void print(int n,struct car*pr)
 {
	 for(int i=1;i<=n;i++)
	 {    if(pr[i].x!=0)
	    	cout<<"*";
	      else
		    cout<<" ";
	 }
	 cout<<endl;
 }