netlim1.cpp

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

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

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

	       if(r<=0&&c>0)    //in Q4 quadrant
	       {
		  if(r==0)
		  {
		     nodes[i].flag=1;
		     nodes[i].nxt1=cnxt*m+rfr;
		     nodes[i].nxt2=cfr*m+rnxt;
		  }
		  if(c==1&&r!=0)
		  {
		     nodes[i].flag=1;
		     nodes[i].nxt1=cfr*m+rnxt;
		     nodes[i].nxt2=cnxt*m+rfr;
		  }
		  if(r!=0&c!=1)
		  {
		     if(r1==m/2)  r1=r1-m;
		     if(c1==-n/2+1)  c1=c1+n;
		     if((r1-r>0&&c1-c<0)||(r1-r<0&&c1-c>0))
		     {
			dcn++;
			nodes[i].flag=0;
			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+cnxt;
			}
		     }
		  }
	       }
	    }
	 }
      }
      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 P1=T*P0
  where T is the transition matrix with two nonzero element in each column
  and row
*************************************************************************/
void TP(double * P1, double * P0)
{
   int i, j;

   for(i=0;i<NN;i++)
      P1[i]=0.0;
   for(i=0;i<NN;i++)
   {
      j=T[i].i1;
      P1[j]=P1[j]+T[i].q1*P0[i];
      j=T[i].i2;
      P1[j]=P1[j]+T[i].q2*P0[i];
   }
}



int mod(int x, int y)
{
   if(x<0)
      x=x%y+y;
   else
      x=x%y;
   return x;
}



double calcupc(int n,  double pu)
{
   int      i;
   double   **a, *b, ep, pcm;

   a=TwoArrayAlloc(n,n);
   b=(double *)calloc(n,sizeof(double));
   if(b==NULL)
   {
      printf("\nMemory calloc fail!");
      exit(1);
   }

   ep=1e-12;

   a[0][0]=1.0-pu*pu-1.0; a[0][1]=1.0-pu;             a[0][2]=1-pu;
   a[1][0]=0.5*pu*pu;     a[1][1]=pu-0.25*pu*pu-1.0;  a[1][2]=0.25*pu*pu;
   a[2][0]=1.0;           a[2][1]=1.0;                a[2][2]=1.0;

   b[0]=0.0; b[1]=0.0; b[2]=1.0;

   Gs2(n, a, b, ep);
//   for(i=0;i<=2;i++)
//      printf("%f\n",b[i]);
   pcm=b[n-1];
   TwoArrayFree(a);
   free(b);
   return(pcm);
}



/************************************************************************
   This program is to solve the linear equations;
************************************************************************/
void Gs2(int n, double **a, double *b, double ep)
{
   int     i, j, k, l;
   double  c, t;

   for(k=1;k<=n;k++)
   {
      c=0.0;
      for(i=k;i<=n;i++)
      {
	 if(fabs(a[i-1][k-1])>fabs(c))
	 {
	    c=a[i-1][k-1];
	    l=i;
	 }
      }
      if(fabs(c)<=ep)
      {
	 printf("\nThe equation has no solution");
	 exit(1);
      }
      if(l!=k)
      {
	 for(j=k;j<=n;j++)
	 {
	    t=a[k-1][j-1];
	    a[k-1][j-1]=a[l-1][j-1];
	    a[l-1][j-1]=t;
	 }
	 t=b[k-1];
	 b[k-1]=b[l-1];
	 b[l-1]=t;
      }
      c=1.0/c;
      for(j=k+1;j<=n;j++)
      {
	 a[k-1][j-1]=a[k-1][j-1]*c;
	 for(i=k+1;i<=n;i++)
	    a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1];
      }
      b[k-1]*=c;
      for(i=k+1;i<=n;i++)
	 b[i-1]-=b[k-1]*a[i-1][k-1];
   }
   for(i=n;i>=1;i--)
      for(j=i+1;j<=n;j++)
	 b[i-1]-=a[i-1][j-1]*b[j-1];
}


double  **TwoArrayAlloc(int r, int c)
{
   double  *x, **y;
   int     n;

   x=(double *)calloc(r*c, sizeof(double));
   y=(double **)calloc(r, sizeof(double *));
   if(x==NULL||y==NULL)
   {
      printf("\nMemory calloc fail");
      exit(1);
   }
   for(n=0;n<=r-1;++n)
      y[n]=&x[c*n];
   return(y);
}


void TwoArrayFree(double **x)
{
   free(x[0]);
   free(x);
}


/***********************************************************************
  This program is to calculate the conditional probability p(e/n)
************************************************************************/
double conber(double R, double l, int n, double m1, int flag)
{
   int     i;
   double  D, belta, wl, Aeff, alpha, Ne, G, G0, Gf, T, tw, tau, sigma, t1, t2;
   double  dt1, dt2, sum, n1, p, pt, pe, E, N0, m, sita;

   n1=(double)n;
   T=1000.0/R;   //T is the pulse interval in ps;
   tau=T/m1;     //tau is the FWHA(full width half amplitude) of the pulse;
   tw=0.4*T;     //2*tw is the width of the detecting width;
   Ne=1.5;       //Ne is the excess noise parameter;
   wl=1.55;      //wl is the wavelengh in um;
   D=1.0;        //D is the dispersion coefficient in ps/nm/km;
   belta=0.53*wl*wl*D;   //belta in psps/km;
   Aeff=35;      //Aeff is the effective fiber core areas in um*um;
   alpha=0.25/4.343;  //alpha is the fiber loss in /km;
   G0=pow(10.0,1.5);  //G0 is the gain compensate for the node loss;
   Gf=exp(alpha*l);   //Gf is the distributed amplifier gain for fiber;
   G=G0*Gf;           //G is the total gain;
   m=8.0;             //2m=2BT+1 is the signal dimensionality;

   sita=3.295*wl*wl*D*n1*l*exp(-0.8814*m1)/(tau*tau);
   dt2=0.2836*tau*(log(2.0)-log(1.0+cos(2.0*sita)));
   //dt2 is the time shift induced by short-range interaction;
   t1=1.25e-5*Ne*(G-1.0)*D*l*l*(n1*n1*n1/3.0+n1*n1/2.0+n1/6.0)/(tau*Aeff);
   //t1 is the Gordon-Haus jitter variance;
   dt1=1.4e-2*belta*belta/pow(tau,4.0)*(l*n-(1.0-exp(-4.0*alpha*l*n))/(4.0*alpha))/(4.0*alpha);
   //dt1 is the time shift induced by soliton self-frequency shift
   sum=0.0;
   for(i=0;i<n;i++)
      sum=sum+(l*(n-i)-(1.0-exp(-4.0*alpha*l*(n-i)))/(4.0*alpha))/(4.0*alpha);
   t2=1.02e-8*pow(wl,4.0)*pow(D,3.0)*Ne*(G-1.0)*sum*sum/(pow(tau,7.0)*Aeff);
   //t2 is the time jitter variance induced by self-frequency shift and ASE
   E=9.30e-3*pow(wl,3.0)*Aeff*D/tau;        //signal energy;
   N0=1.99e-7*n1*Ne*(G-1.0)/wl;

   if(flag==1)
   {
      sigma=sqrt(t1+t2);
      pt=0.5*(qfunc((tw-dt1)/sigma)+qfunc((tw+dt1)/sigma));
      pt=pt+0.25*(qfunc((tw-dt1-dt2)/sigma)+qfunc((tw+dt1+dt2)/sigma)+qfunc((tw-dt1+dt2)/sigma)+qfunc((tw+dt1-dt2)/sigma));
      p=0.5*pt;
   }
   if(flag==2)
   {
      sigma=sqrt(t1);
      pt=qfunc(tw/sigma)+0.5*qfunc((tw-dt2)/sigma)+0.5*qfunc((tw+dt2)/sigma);
      p=0.5*pt;
   }
   if(flag==3)
   {
      sigma=sqrt(t2);
      pt=0.5*(qfunc((tw-dt1)/sigma)+qfunc((tw+dt1)/sigma));
      pt=pt+0.25*(qfunc((tw-dt1-dt2)/sigma)+qfunc((tw+dt1+dt2)/sigma)+qfunc((tw-dt1+dt2)/sigma)+qfunc((tw+dt1-dt2)/sigma));
      p=0.5*pt;
   }
   if(flag==4)
      p=qfunc(E/N0/(sqrt(m)+sqrt(m+2.0*E/N0)));
   return(p);
}
/***********************************************************************/



/***********************************************************************
  This program is to calculate the conditional intersymbol crosstalk
************************************************************************/
double conint(double R, double l, int n)
{
   int     i;
   double  D, wl, tau0, tau, T, tw, n1, p, belta, s, x;

   n1=(double)n;
   T=1000.0/R;   //tau is the FWHA(full width half amplitude) of the pulse;
   tau0=T/2.0;   //2*tw is the width of the detecting width;
   tw=0.4*T;
   wl=1.55;      //wl is the wavelengh in um;
   D=1.0;        //D is the dispersion coefficient in ps/nm/km;
   belta=0.53*wl*wl*D;
   s=4.0*log(2.0)*belta*n1*l/(tau0*tau0);
   tau=tau0*sqrt(1.0+s*s);
   x=2.0*sqrt(log(2.0))*(T-tw)/tau;
   p=erfc(x);
   return(p);
}
/***********************************************************************/




/***********************************************************************
     The following programs are to calculate the Q function
************************************************************************/
double qfunc(double x)
{
  if(x>=30)
     return(0.0);
  else
     return(simp2(fun,x,30,1.0e-8));
}
/***********************************************************************/



/***********************************************************************
     The following program is to calculate the erfc function
************************************************************************/
double erfc(double x)
{
  if(x>=30)
     return(0.0);
  else
     return(simp2(fun1,x,30,1.0e-8));
}
/**********************************************************************/



double simp2(double (*funcp)(double), double a, double b, double eps)
{
   int    i, j, k, n;
   double h, t1, t2, s1, s2, p, x, judge;

   n=1;
   h=b-a;
   t1=h*(funcp(a)+funcp(b))/2.0;
   s1=t1;
   while(1)
   {
      p=0;
      for(k=0;k<=n-1;k++)
      {
	 x=a+(k+0.5)*h;
	 p+=funcp(x);
      }
      t2=(t1+h*p)/2.0;
      s2=(4*t2-t1)/3.0;
      if(s1!=0)
	 judge=fabs((s2-s1)/s1);
      else
	 judge=fabs(s2-s1);
      if(judge>=eps)
      {
	 t1=t2;
	 n=n+n;
	 h=h/2;
	 s1=s2;
	 continue;
      }
      break;
   }
   return(s2);
}


double fun(double x)
{
   return(exp(-x*x/2)/sqrt(2.0*PI));
}



double fun1(double x)
{
   return(2.0*exp(-x*x)/sqrt(PI));
}
/***********************************************************************/



/************************************************************************
   This program is to determine the soliton pulse width, that is T/tau
************************************************************************/
double range(double R, double z)
{
   double  m, T, a, h, fx;

   a=5.0;
   h=1.0;
   T=1000.0/R;
   fx=funx(a,T,z);
   if(fx>=0.0)
   {
      m=a;
      goto end;
   }
   else
   {
      while(1)
      {
	 do
	 {
	    a=a+h;
	    fx=funx(a,T,z);
	 }while(fx<0.0);
	 while(1)
	 {
	    h=0.5*h;
	    a=a-h;
	    fx=funx(a,T,z);
	    if(fx>=0)
	    {
	       if(h<0.01)
	       {
		  m=a;
		  goto end;
	       }
	    }
	    else
	       break;
	 }
	 h=0.5*h;
      }
   }
   end:
   return(m);
}



double funx(double m, double T, double z)
{
  double wl, D;

  wl=1.55;
  D=1.0;
  return(0.1430*T*T*exp(0.8814*m)/(m*m*wl*wl*D)-z);
}
/**************************************************************************/