netlim1.cpp

分步傅利叶方法解非线性方程

/* This program is to calculate the BER(Bit Error Rate) of shuffleNet and
   Manhattan Street Network with single-buffer deflection routing.
    p(e)=p(e/1)p(1)+...+p(e/n)p(n)+...;
    where p(n) is the hop distribution probability;    */


# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <conio.h>
# include <iostream.h>

struct snode
{
     char          flag;   //if flag=1, care node, else don't care node
     unsigned int  nxt1;   //if flag=1, nxt1 is the preferred path
     unsigned int  nxt2;
} ;

struct element      // It is a structure of transition matrix element
{
    char    flag;   // if flag=1, care node
    int     i1;     // if flag=1, i1 is the preferred node, q1 and q2 are
    int     i2;     // the probability to i1 and i2 respectively
    double   q1;     // if flag=1, q1=1-p, q2=p; else q1=q2=0.5;
    double   q2;
};

struct element   T[1024];  //T is the transition probability matrix
int   dc[1024];            //It stores the number of don't care node
int   NN;  	           //The total number of nodes
int   dcn;                 //The number of care nodes

double const  PI=3.1415926;

void ShuffleNet();
void MSN();
void TP(double *, double *);  //TP is the function T*P, that is matrix time vector
int mod(int x, int y);
double  calcupc(int, double);
void  Gs2(int, double **, double *, double);
double  **TwoArrayAlloc(int, int);
void  TwoArrayFree(double **);
double  conber(double, double, int, double, int); //conditional BER p(e/n)
double  conint(double, double, int);      //conditional intersymbol interfere
double  qfunc(double);  //Q function;
double  simp2(double (*funcp)(double),double, double, double);//simpleson integration
double  fun(double);   //used by qfunc()
double  erfc(double);
double  fun1(double);  //used by erfc();
double range(double, double); //to determine the soliton pulse width T/tau;
double funx(double, double, double);//used by range();




void main()
{
   char     *fnam, *fnam1;
   int      ch, flag, flag1, i, j, h;
   double   pdc2, pdc1, D, g, u, a, pc, pc0, pa0, pa1, pu, pcm, judge;
   double   P1[1024], P0[1024];   //P1 and P0 store the P(k) and P(k-1) respectivly
   double   z, m, l, R, pber, pcon, pn;
   //l is the distance between nodes in km;
   //R is the bit rate in Gb/s;
   //pber is the total BER, and pcon is the conditional BER p(e/n);
   FILE     *fp, *fp1;

   clrscr();
   gotoxy(1,5);
   printf("     -----------------------------------------------------------\n");
   printf("     |   This program is to calculate the BER(Bit Error Rate)  |\n");
   printf("     |   of ShuffleNet and Manhattan Street Network with       |\n");
   printf("     |   single-buffer deflection routing. Where l=1km         |\n");
   printf("     |   Designed by Xie Chongjin in Sep. 1997.                |\n");
   printf("     -----------------------------------------------------------\n");

   do
   {
      printf("\nPlease choise the ShuffleNet or Manhattan Street Network");
      printf("\n1.ShuffleNet;  2.Manhattan Street Network:   ");
      scanf("%d",&ch);
   }while(ch!=1&&ch!=2);

   if(ch==1)
      ShuffleNet();
   if(ch==2)
      MSN();

   do
   {
      printf("\n1.Soliton system. 2.Linear system.   Please chose:");
      scanf("%d",&flag1);
   }while(flag1!=1&&flag1!=2);

   if(flag1==1)
   {
      flag=1;
      do
      {
	 printf("\nPlease choise considered noise sources:");
	 printf("\n1.All noise sources.\n2.Gordon-Haus jitter only:");
	 printf("\n3.Roman self frequency shift induced by ASE noices only:");
	 printf("\n4.Energy BER produced by ASE noices only:\n");
	 scanf("%d",&flag);
      }while(flag!=1&&flag!=2&&flag!=3&&flag!=4);
   }

   fnam=(char *)malloc(80);
   cout<<"\nPlease give the output data file name:\n";
   scanf("%s",fnam);
   if((fp=fopen(fnam,"wt"))==NULL)
   {
      printf("file %s create error!\n",fnam);
      exit(1);
   }
   free(fnam);
   if(flag1==1)
   {
      fnam1=(char *)malloc(80);
      cout<<"\nPlease give the output pulse width file name:\n";
      scanf("%s",fnam1);
      if((fp1=fopen(fnam1,"wt"))==NULL)
      {
	 printf("file %s create error!\n",fnam1);
	 exit(1);
      }
      free(fnam1);
   }

   printf("\nPlease wait!!!");

//calculate network performance with the g
      g=1.0;
      pdc1=0.0;
      pdc2=0.0;
      pc=1.0;
      pc0=1.0;
      pcm=1.0;
      //calculate pdc;
      do
      {
	 pdc1=pdc2;
	 P1[0]=0.0;
	 for(i=1;i<NN;i++)
	    P1[i]=1.0/(double)(NN-1);
	 for(i=1;i<NN;i++)
	 {
	    if(T[i].flag==1)
	    {
	       T[i].q2=0.25*pc0*pcm*0.5*(1-pdc1);
	       T[i].q1=1.0-T[i].q2;
	    }
	 }

	 h=0;
	 D=0.0;
	 pdc2=0.0;
	 for(i=0;i<NN;i++)
	    P0[i]=P1[i];
	 TP(P1,P0);
	 h++;
	 D=D+(double)h*(P1[0]-P0[0]);
	 for(i=0;i<dcn;i++)
	 {
	    j=dc[i];
	    pdc2=pdc2+P0[j];
	 }

	 for(i=1;i<NN;i++)
	 {
	    if(T[i].flag==1)
	    {
	       T[i].q2=0.25*pc*pcm*0.5*(1-pdc1);
	       T[i].q1=1.0-T[i].q2;
	    }
	 }

	 do
	 {
	    for(i=0;i<NN;i++)
	       P0[i]=P1[i];
	    TP(P1,P0);
	    h++;
	    D=D+(double)h*(P1[0]-P0[0]);
	    for(i=0;i<dcn;i++)
	    {
	       j=dc[i];
	       pdc2=pdc2+P0[j];
	    }
	    judge=fabs(P1[0]-P0[0]);
	    for(i=1;i<NN;i++)
	    {
	       if(judge<fabs(P1[i]))
		  judge=fabs(P1[i]);
	    }
	 }while(judge>1.0e-14);

	 pdc2=pdc2/D;
	 a=1.0/D;
	 u=(sqrt(a*a+g*g*(1-a)*(1-a))-a)/(g*(1-a)*(1-a));
	 pc=u*(1-a)+u*a*g+(1-u)*g;
	 pu=u*(1-pdc1);
	 pcm=2.0*calcupc(3,pu);
	 pa0=(1-u)*(1-u)+2*u*(1-u)*a+u*u*a*a;
	 pa1=2*u*(1-u)*(1-a)+2*u*u*a*(1-a);
	 pc0=pa1/(pa0+pa1);

      }while(fabs(pdc2-pdc1)>1.0e-4);

//calculate the hop number that p(n)=1.0e-14
      pdc1=pdc2;
      P1[0]=0.0;
      for(i=1;i<NN;i++)
	 P1[i]=1.0/(double)(NN-1);
      for(i=1;i<NN;i++)
      {
	 if(T[i].flag==1)
	 {
	    T[i].q2=0.25*pc0*pcm*0.5*(1-pdc1);
	    T[i].q1=1.0-T[i].q2;
	 }
      }

      h=0;
      for(i=0;i<NN;i++)
	 P0[i]=P1[i];
      TP(P1,P0);
      h++;

      for(i=1;i<NN;i++)
      {
	 if(T[i].flag==1)
	 {
	    T[i].q2=0.25*pc*pcm*0.5*(1-pdc1);
	    T[i].q1=1.0-T[i].q2;
	 }
      }

      do
      {
	 for(i=0;i<NN;i++)
	    P0[i]=P1[i];
	 TP(P1,P0);
	 h++;
	 judge=fabs(P1[0]-P0[0]);
      }while(judge>1.0e-13);
      l=1.0;
      z=l*(double)h;

      //calculate the hop distribution and BER
      if(flag1==1)
      {
	 for(R=20;R<=200;R=R+10)
	 {
	    m=range(R,z);
	    pber=0.0;
	    pdc1=pdc2;
	    P1[0]=0.0;
	    for(i=1;i<NN;i++)
	       P1[i]=1.0/(double)(NN-1);
	    for(i=1;i<NN;i++)
	    {
	       if(T[i].flag==1)
	       {
		  T[i].q2=0.25*pc0*pcm*0.5*(1-pdc1);
		  T[i].q1=1.0-T[i].q2;
	       }
	    }

	    h=0;
	    for(i=0;i<NN;i++)
	       P0[i]=P1[i];
	    TP(P1,P0);
	    h++;
	    pn=P1[0]-P0[0];
	    pcon=conber(R,l,h,m,flag);
	    pber=pber+pn*pcon;

	    for(i=1;i<NN;i++)
	    {
	       if(T[i].flag==1)
	       {
		  T[i].q2=0.25*pc*pcm*0.5*(1-pdc1);
		  T[i].q1=1.0-T[i].q2;
	       }
	    }

	    do
	    {
	       for(i=0;i<NN;i++)
		  P0[i]=P1[i];
	       TP(P1,P0);
	       h++;
	       pn=P1[0]-P0[0];
	       pcon=conber(R,l,h,m,flag);
	       pber=pber+pn*pcon;
	       judge=fabs(pn);
	       for(i=1;i<NN;i++)
	       {
		  if(judge<fabs(P1[i]))
		     judge=fabs(P1[i]);
	       }

	    }while(judge>1.0e-14);
	    printf("R=%lf has been calculated\n",R);
	    fprintf(fp,"%lf\t%e\n",R,pber);
	    fprintf(fp1,"%lf\t%lf\n",R,m);
	 }
	 fclose(fp1);
      }
      //calculate the intersymbol crosstalk
      if(flag1==2)
      {
	 for(R=5;R<=200;R=R+5)
	 {
	    l=1.0;
	    pber=0.0;
	    pdc1=pdc2;
	    P1[0]=0.0;
	    for(i=1;i<NN;i++)
	       P1[i]=1.0/(double)(NN-1);
	    for(i=1;i<NN;i++)
	    {
	       if(T[i].flag==1)
	       {
		  T[i].q2=0.25*pc0*pcm*0.5*(1-pdc1);
		  T[i].q1=1.0-T[i].q2;
	       }
	    }

	    h=0;
	    for(i=0;i<NN;i++)
	       P0[i]=P1[i];
	    TP(P1,P0);
	    h++;
	    pn=P1[0]-P0[0];
	    pcon=conint(R,l,h);
	    pber=pber+pn*pcon;

	    for(i=1;i<NN;i++)
	    {
	       if(T[i].flag==1)
	       {
		  T[i].q2=0.25*pc*pcm*0.5*(1-pdc1);
		  T[i].q1=1.0-T[i].q2;
	       }
	    }

	    do
	    {
	       for(i=0;i<NN;i++)
		  P0[i]=P1[i];
	       TP(P1,P0);
	       h++;
	       pn=P1[0]-P0[0];
	       pcon=conint(R,l,h);
	       pber=pber+pn*pcon;
	       judge=fabs(pn);
	       for(i=1;i<NN;i++)
	       {
		  if(judge<fabs(P1[i]))
		     judge=fabs(P1[i]);
	       }

	    }while(judge>1.0e-14);
	    printf("R=%lf has been calculated\n",R);
	    fprintf(fp,"%lf\t%e\n",R,pber);
	 }
      }

      fclose(fp);
      printf("\nOK!!!");
}


/***********************************************************************
  This program is to calculate the transition probability matrix of
  ShuffleNet
***********************************************************************/
void ShuffleNet()
{
   int            i, j, k, r, c, r1, c1, rmax, d, base;
   struct snode   nodes[1024];

   printf("Please input paremeter k(1-7):");
   scanf("%d",&k);
   while(k<1||k>7)
   {
      printf("The k inputed is beyond the scope, please input again k(1-7):");
      scanf("%d",&k);
   }

      dcn=0;
      rmax=(int)pow(2.0,(double)k);
      NN=k*rmax;
      base=rmax-1;
      for(c=0;c<k;c++)
      {
	 for(r=0;r<rmax;r++)
	 {
	    i=c*rmax+r;
	    if(c!=0||r!=0)
	    {
	       if(c==0)
		  d=k;
	       else
		  d=mod(k-c,k);
	       r1=r<<d;
	       r1=r1&base;
	       if(r1==0)
		  nodes[i].flag=1;
	       else
	       {
		  nodes[i].flag=0;
		  dcn++;
	       }
	       c1=mod(c+1,k);
	       r1=r<<1;
	       r1=r1&base;
	       nodes[i].nxt1=c1*rmax+r1;
	       nodes[i].nxt2=nodes[i].nxt1+1;
	    }
	 }
      }
      T[0].flag=0;
      T[0].i1=0;
      T[0].i2=1;
      T[0].q1=1.0;
      T[0].q2=0.0;
      j=0;
      for(i=1;i<NN;i++)
      {
	 if(nodes[i].flag==0)
	 {
	    T[i].flag=0;
	    T[i].i1=nodes[i].nxt1;
	    T[i].i2=nodes[i].nxt2;
	    T[i].q1=0.5;
	    T[i].q2=0.5;
	    dc[j]=i;
	    j++;
	 }
	 else
	 {
	    T[i].flag=1;
	    T[i].i1=nodes[i].nxt1;
	    T[i].i2=nodes[i].nxt2;
	    T[i].q1=0.75;
	    T[i].q2=0.25;
	 }
      }
}



/***********************************************************************
  This program is to calculate the transition probability matrix of
  Manhattan Street Network
***********************************************************************/
void MSN()
{
   int            m, n, rfr, cfr, r, c, rnxt, cnxt, r1, c1, i, j;
   struct snode   nodes[1024];

   printf("Please input paremeter n(n must be even 2-32):");
   scanf("%d",&n);
   while(n<2||n>32||n%2==1)
   {
      printf("The n inputed is wrong, please input again n(n must be even 2-32):");
      scanf("%d",&n);
   }

      dcn=0;
      m=n;
      NN=m*n;
      for(cfr=0;cfr<n;cfr++)
      {
	 for(rfr=0;rfr<m;rfr++)
	 {
	    if(rfr!=0||cfr!=0)
	    {
	       if(cfr%2==0)
		  rnxt=mod(rfr+1,m);
	       else
		  rnxt=mod(rfr-1,m);
	       if(rfr%2==0)
		  cnxt=mod(cfr+1,n);
	       else
		  cnxt=mod(cfr-1,n);
	       i=cfr*m+rfr;
	       r=m/2-mod(m/2+rfr,m);
	       c=n/2-mod(n/2+cfr,n);
	       r1=m/2-mod(m/2+rnxt,m);
	       c1=n/2-mod(n/2+cnxt,n);
//             judge the quadrant

	       if(r>0&&c>0)       //in Q1 quadrant
	       {
		  if(r==m/2&&c==n/2)
		  {
		     nodes[i].flag=0;
		     dcn++;
		     nodes[i].nxt1=cfr*m+rnxt;
		     nodes[i].nxt2=cnxt*m+rfr;
		  }
		  if(r==m/2&&c<n/2)
		  {
		     if(c1-c<0)
		     {
			nodes[i].flag=0;
			dcn++;
			nodes[i].nxt1=cfr*m+rnxt;
			nodes[i].nxt2=cnxt*m+rfr;
		     }
		     else
		     {
			nodes[i].flag=1;
			nodes[i].nxt1=cfr*m+rnxt;
			nodes[i].nxt2=cnxt*m+rfr;
		     }
		  }
		  if(r<m/2&&c==n/2)
		  {
		     if(r1-r<0)
		     {
			nodes[i].flag=0;
			dcn++;
			nodes[i].nxt1=cfr*m+rnxt;
			nodes[i].nxt2=cnxt*m+rfr;
		     }
		     else
		     {
			nodes[i].flag=1;
			nodes[i].nxt1=cnxt*m+rfr;
			nodes[i].nxt2=cfr*m+rnxt;
		     }
		  }
		  if(r<m/2&&c<n/2)
		  {
		     if((r1-r<0&&c1-c<0)||(r1-r>0&&c1-c>0))
		     {
			nodes[i].flag=0;
			dcn++;
			nodes[i].nxt1=cfr*m+rnxt;
			nodes[i].nxt2=cnxt*m+rfr;
		     }
		     else
		     {
			if(r1-r<0)
			{
			   nodes[i].flag=1;
			   nodes[i].nxt1=cfr*m+rnxt;
			   nodes[i].nxt2=cnxt*m+rfr;
			}
			if(c1-c<0)
			{
			   nodes[i].flag=1;
			   nodes[i].nxt1=cnxt*m+rfr;
			   nodes[i].nxt2=cfr*m+rnxt;
			}
		     }
		  }
	       }