openboundary.c

开发边界条件下用元胞自动机模拟交通流的基本情况

 /////*改进的（谭惠丽：《2002年物理学报》）NS-MODEL 在开放性边界条件下程序*////////
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <stdio.h>
#define Vmax  2
#define N  1024
#define T  10000
int Max(int x, int y)
{
    int z;
    z=x>y?x:y;
    return(z);
} 
int Min(int a, int b)
{
	int z;
    z=a<b?a:b;
    return (z);
}
int rule(int a[],int x[],int v[])    //////*0-1024个格点内的车子速度，位置演化函数*////////
{  
    int i,k,n,gap,m=0;
	int local[N+2],vel[N+2],kk[N+2];
    double p;
	for(i=0;i<N+2;i++) { local[i]=0;   vel[i]=0;}
	for (i=0;i<N+1;i++)///速度更新
	{   
		if(a[i]==1)
		{
                         n=0;
			 k=i+1;
			 while(k<N+2)
			 {
				 if(a[k]==1)
				 { n=k; break;}
				 k++;
			 }
			 if(n==0)     { gap=Vmax+1;   v[k-1]=Vmax;}//没车挡
			 else         { gap=n-i-1;    v[N+1]=0;}//有车挡
		        p=(double)(rand())/(double)(RAND_MAX);
	
		        if(v[i]>=gap+v[k-1])                                     ////条件一   
			{   
			    if(p<=0.5)            vel[i]=Max(gap+v[k-1]-1,0);
		            else                  vel[i]=gap+v[k-1];
			}
		        else                                                    ////条件二
			{ 
			    if(p<=0.5)           vel[i]=v[i];
			    else                 vel[i]=Min(v[i]+1,Vmax);
			}
		}
	}///

	for (i=0;i<N+1;i++)///位置更新
	{   
		if(a[i]==1)
		{
			local[i]=x[i]+vel[i];  
		}
	}////
   
	for(i=0;i<N+2;i++)          { kk[i]=0;  x[i]=0;  v[i]=0;  }
	for(i=0;i<N+2;i++) 
	{
		if(a[i]==1)
		{
	         if(local[i]<=N)     
			 { kk[local[i]]=1; x[local[i]]=local[i]; v[local[i]]=vel[i];}
		}
	}
    for(i=0;i<N+2;i++)    a[i]=kk[i];
    for(i=1;i<=N;i++)   
	if(a[i]==1)  m=m+1;
    return(m);  
}
			        
///////*主函数*////////////////////////
void main ()
{ 
	FILE*fp1;
    long t;
    int i,sum,sample;   
    double Pin,Pout=1.00,m,n;
    int a[N+2],x[N+2],v[N+2];
    int rand(void);
    srand(time(NULL)); 
	fp1=fopen("e:\\覃惠丽-density(2v-0.5).txt","a+");
    for(Pin=0.05;Pin<1.0005;Pin=Pin+0.05) /////*左端注人汽车的概率*////////
    {   
		double  G1=0.00;
		for(sample=0;sample<20;sample++)           ///*20个样品*///////
		{ 	
			double av1=0.00;				    	  
            for(i=0;i<N+2;i++)    ////初值
			{a[i]=0;x[i]=0;v[i]=0;}
	  	
            for(t=0;t<T;t++) /////*汽车位置，速度按时步演化*////////
			{  
				sum=0;
			    n=(double)(rand())/(double)(RAND_MAX); /////////*在0位，按概率Pin产生车子*////////
	            m=(double)(rand())/(double)(RAND_MAX); /////////*在N+1位，按概率1-Pout产生阻挡车子*////////
			
		    if(m<=Pout )    a[N+1]=0;
	            else            a[N+1]=1;
		        
                if(n<=Pin)       //////*第一种情况以Pin的概率产生车子*////////
				{  	
				 
				    a[0]=1;  x[0]=0;  v[0]=Vmax;
				
                    sum=rule(a,x,v);
				
				    if(a[0]==1)   /////////*如果注入的车子在t时步速度为零，消除它*/////////
					a[0]=0;
				}
                else     sum=rule(a,x,v); 
				
				if(t>=T-1000)      /////////////*对最后的10000步的密度求平均*///////////        
				av1=av1+sum/(double)(1000*N);
				
			}					
		   G1=G1+av1/(double)(20);    /////////////*对20个样品求平均*//////////
		     
		}
        fprintf(fp1,"%8f,%8f\n",Pin,G1);
	}
	fclose(fp1);
}