串行源程序.txt

串行运算的原代码,单机实现计算,程序写的很不错,

    for(i=0;i<=14;i++)
     for(nl=6;nl>=1;nl--)
     { i1=i+tx0;
       nl0=nl+m0;
       zez2[k][i1][nl0]=zez2[k][i1][nl0-1];
       zex2[k][i1][nl0]=zex2[k][i1][nl0-1]; }
   for(k=0;k<=Z1;k++)
    for(i=0;i<=14;i++)
    { i1=i+tx;
      if(yy==0) i2=yy+i;
      else i2=yy-i;
      zey1[k][i1][m]=ey1[xx][i2][k];
      zez1[k][i1][m]=ez1[xx][i2][k];}
   for(k=0;k<=Z1;k++)
    for(i=0;i<=14;i++)
    { i1=i+tx0;
      if(xx==0) i2=xx+i;
      else i2=xx-i;
      zez2[k][i1][m0]=ez1[i2][yy][k];
      zex2[k][i1][m0]=ex1[i2][yy][k]; }
 }

 void liaox(int xs,int xf,int yy,int zz)
 { int m,m0,tx,tx0,nl,nl0,ml,ml0,lm,l,i,j,i1,i2;
   float r,temp,tempx1[X1+1],tempx2[X1+1],tempy[X1+1],tempz[X1+1];
   if (yy==0) {m=0;tx0=0;}
   else {m=7;tx0=15;}
   if (zz==Z1) {m0=0;tx=0;}
   else {m0=7;tx=15;}
   for (i=xs;i<=xf;i++)
    {if(zz==Z2) i1=i+30;
     else i1=i;
     tempy[i1]=0;
     nl=m;
     for(l=1;l<=7;l++)
      {r=0;
       lm=2*l;
       for(ml=0;ml<=lm;ml++)
        {ml0=ml+tx;
         r=r+tt[l][ml]*xey1[i1][ml0][nl];}
       temp=pow(-1,l+1)*c7[l]*r;
       tempy[i1]=tempy[i1]+temp;
       nl++;} }
   for(i=xs;i<=xf;i++)
     { if(zz==Z2) ey2[i][yy][zz]=tempy[i+30];
       else ey1[i][yy][zz]=tempy[i]; }
   for (i=xs;i<=xf;i++)                   
    {if(zz==Z2) i1=i+30;
     else i1=i;
     tempz[i1]=0;
     nl=m0;
     for(l=1;l<=7;l++)
      {r=0;
       lm=2*l;
       for(ml=0;ml<=lm;ml++)
        {ml0=ml+tx0;
         r=r+tt[l][ml]*xez2[i1][ml0][nl];}
       temp=pow(-1,l+1)*c7[l]*r;
       tempz[i1]=tempz[i1]+temp;
       nl++;} }
   for(i=xs;i<=xf;i++)
     { if(zz==Z2) ez2[i][yy][zz]=tempz[i+30];
       else ez1[i][yy][zz]=tempz[i]; }
   for (i=xs;i<=xf;i++)
    {if(zz==Z2) i1=i+30;
     else i1=i;
     tempx1[i1]=0;
     nl=m;
     for(l=1;l<=7;l++)
      {r=0;
       lm=2*l;
       for(ml=0;ml<=lm;ml++)
        {ml0=ml+tx;
         r=r+tt[l][ml]*xex1[i1][ml0][nl];}
       temp=pow(-1,l+1)*c7[l]*r;
       tempx1[i1]=tempx1[i1]+temp;
       nl++;} }
   for (i=xs;i<=xf;i++)
    {if(zz==Z2) i1=i+30;
     else i1=i;
     tempx2[i1]=0;
     nl=m0;
     for(l=1;l<=7;l++)
      {r=0;
       lm=2*l;
       for(ml=0;ml<=lm;ml++)
        {ml0=ml+tx0;
         r=r+tt[l][ml]*xex1[i1][ml0][nl];}
       temp=pow(-1,l+1)*c7[l]*r;
       tempx2[i1]=tempx2[i1]+temp;
       nl++;} }
   for(i=xs;i<=xf;i++)
     { if(zz==Z2) ex2[i][yy][zz]=(tempx1[i+30]+tempx2[i+30])/2;
       else ex1[i][yy][zz]=(tempx1[i]+tempx2[i])/2; }
   for(i=0;i<=X1;i++)
    for(j=0;j<=14;j++)
     for(nl=6;nl>=1;nl--)
     { i1=j+tx;
       nl0=nl+m;
       xey1[i][i1][nl0]=xey1[i][i1][nl0-1];
       xex1[i][i1][nl0]=xex1[i][i1][nl0-1]; }
   for(i=0;i<=X1;i++)
    for(j=0;j<=14;j++)
     for(nl=6;nl>=1;nl--)
     { i1=j+tx0;
       nl0=nl+m0;
       xez2[i][i1][nl0]=xez2[i][i1][nl0-1];
       xex2[i][i1][nl0]=xex2[i][i1][nl0-1]; }
   for(i=xs;i<=xf;i++)
    for(j=0;j<=14;j++)
    { i1=j+tx;
      if(yy==0) i2=yy+i;
      else i2=yy-i;
      if (zz==Z2){xey1[i+30][i1][m]=ey2[i][i2][zz];
                  xex1[i+30][i1][m]=ex2[i][i2][zz];}
      else {xey1[i][i1][m]=ey1[i][i2][zz];
            xex1[i][i1][m]=ex1[i][i2][zz];}    }
   for(i=xs;i<=xf;i++)
    for(j=0;j<=14;j++)
    { i1=j+tx0;
      if(yy==0) i2=yy+i;
      else i2=yy-i;
      if (zz==Z2){xez2[i+30][i1][m0]=ez2[i][i2][zz];
                  xex2[i+30][i1][m0]=ex2[i][i2][zz];}
      else {xez2[i][i1][m0]=ez1[i][i2][zz];
            xex2[i][i1][m0]=ex1[i][i2][zz];}    }

 }

 void  fourier(float e[X2+1][Y+1][rn+1])
  {  int    n,k,i,j;
     float  gk,gre[X2+1][Y+1][4],gim[X2+1][Y+1][4],a[rn+1];
     gk=8.5/df;
     for (k=0;k<=3;k++)
      for(i=0;i<=X2;i++)
       for(j=0;j<=Y;j++)
          {gre[i][j][k]=0;
           gim[i][j][k]=0;}
     for(k=0;k<=3;k++)
      { for(n=0;n<=rn;n++)
        { a[n]=(2*pi*gk*n)/rn; }
        for(i=0;i<=X2;i++)
         for(j=0;j<=Y;j++)
          for(n=0;n<=rn;n++)
            {gre[i][j][k]=gre[i][j][k]+e[i][j][n]*cos(a[n]);
             gim[i][j][k]=gim[i][j][k]+e[i][j][n]*sin(a[n]); }
        gk=gk+1/df; }
     if (e==rex)
        {  for (i=0;i<=X2;i++)
            for (j=0;j<=Y;j++)
                { fex[i][j][0]=gre[i][j][1];
                  fex[i][j][1]=gim[i][j][1];}    }
     if (e==rey)
        {  for (i=0;i<=X2;i++)
            for (j=0;j<=Y;j++)
                { fey[i][j][0]=gre[i][j][1];
                  fey[i][j][1]=gim[i][j][1];}    }
     if (e==rhx)
        {  for (i=0;i<=X2;i++)
            for (j=0;j<=Y;j++)
                { fhx[i][j][0]=gre[i][j][1];
                  fhx[i][j][1]=gim[i][j][1];}    }
     if (e==rhy)
        {  for (i=0;i<=X2;i++)
            for (j=0;j<=Y;j++)
                { fhy[i][j][0]=gre[i][j][1];
                  fhy[i][j][1]=gim[i][j][1];}    }
  }

 void sss()
{ int a,b,c0;
  float k,gr1,gi1,gr2,gi2,temp;
  float gr[X2+1][Y+1][46],gi[X2+1][Y+1][46];
  k=2*pi/ld;
  for(c0=-45;c0<=-1;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=fey[a][b][0]-ita*fhx[a][b][0]*cos(c0*dc);
       gi1=fey[a][b][1]-ita*fhx[a][b][1]*cos(c0*dc);
       temp=k*((-1)*(b-50)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0+45]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0+45]=(gr1*gi2+gi1*gr2);
       Er1[a][b][c0+45]=gi[a][b][c0+45]*(k/(4*pi));
       Ei1[a][b][c0+45]=(-1)*gr[a][b][c0+45]*(k/(4*pi));  }
  for(c0=0;c0<=45;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=-fey[a][b][0]+ita*fhx[a][b][0]*cos(c0*dc);
       gi1=-fey[a][b][1]+ita*fhx[a][b][1]*cos(c0*dc);
       temp=k*((b-50)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0]=(gr1*gi2+gi1*gr2);
       Er1[a][b][c0+45]=gi[a][b][c0]*(k/(4*pi));
       Ei1[a][b][c0+45]=(-1)*gr[a][b][c0]*(k/(4*pi));}
  for(c0=0;c0<=90;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { se1[c0][0]=se1[c0][0]+Er1[a][b][c0];
       se1[c0][1]=se1[c0][1]+Ei1[a][b][c0];}
  for(c0=-45;c0<=-1;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=fex[a][b][0]*cos(c0*dc)+ita*fhy[a][b][0];
       gi1=fex[a][b][1]*cos(c0*dc)+ita*fhy[a][b][1];
       temp=k*((-1)*(b-50)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0+45]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0+45]=(gr1*gi2+gi1*gr2);
       Er2[a][b][c0+45]=(-1)*gi[a][b][c0+45]*(k/(4*pi));
       Ei2[a][b][c0+45]=gr[a][b][c0+45]*(k/(4*pi));    }
  for(c0=0;c0<=45;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=-fex[a][b][0]*cos(c0*dc)-ita*fhy[a][b][0];
       gi1=-fex[a][b][1]*cos(c0*dc)-ita*fhy[a][b][1];
       temp=k*((b-50)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0]=(gr1*gi2+gi1*gr2);
       Er2[a][b][c0+45]=(-1)*gi[a][b][c0]*(k/(4*pi));
       Ei2[a][b][c0+45]=gr[a][b][c0]*(k/(4*pi));       }
  for(c0=0;c0<=90;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { se2[c0][0]=se2[c0][0]+Er2[a][b][c0];
       se2[c0][1]=se2[c0][1]+Ei2[a][b][c0];}
  for(c0=0;c0<=90;c0++)
   {  Me1[c0]=sqrt(pow(se1[c0][0],2)+pow(se1[c0][1],2));
      Me2[c0]=sqrt(pow(se2[c0][0],2)+pow(se2[c0][1],2)); }
}

 void sss1()
{ int a,b,c0;
  float k,gr1,gi1,gr2,gi2,temp;
  float gr[X2+1][Y+1][46],gi[X2+1][Y+1][46];
  k=2*pi/ld;
  for(c0=-45;c0<=-1;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=-fex[a][b][0]-ita*fhy[a][b][0]*cos(c0*dc);
       gi1=-fex[a][b][1]-ita*fhy[a][b][1]*cos(c0*dc);
       temp=k*((-1)*(a-10)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0+45]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0+45]=(gr1*gi2+gi1*gr2);
       Er3[a][b][c0+45]=gi[a][b][c0+45]*(k/(4*pi));
       Ei3[a][b][c0+45]=(-1)*gr[a][b][c0+45]*(k/(4*pi));  }
  for(c0=0;c0<=45;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=fex[a][b][0]+ita*fhy[a][b][0]*cos(c0*dc);
       gi1=fex[a][b][1]+ita*fhy[a][b][1]*cos(c0*dc);
       temp=k*((a-10)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0]=(gr1*gi2+gi1*gr2);
       Er3[a][b][c0+45]=gi[a][b][c0]*(k/(4*pi));
       Ei3[a][b][c0+45]=(-1)*gr[a][b][c0]*(k/(4*pi));}
  for(c0=0;c0<=90;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { se3[c0][0]=se3[c0][0]+Er3[a][b][c0];
       se3[c0][1]=se3[c0][1]+Ei3[a][b][c0];}
  for(c0=-45;c0<=-1;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=fey[a][b][0]*cos(c0*dc)-ita*fhx[a][b][0];
       gi1=fey[a][b][1]*cos(c0*dc)-ita*fhx[a][b][1];
       temp=k*((-1)*(a-10)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0+45]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0+45]=(gr1*gi2+gi1*gr2);
       Er4[a][b][c0+45]=(-1)*gi[a][b][c0+45]*(k/(4*pi));
       Ei4[a][b][c0+45]=gr[a][b][c0+45]*(k/(4*pi));    }
  for(c0=0;c0<=45;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { gr1=-fey[a][b][0]*cos(c0*dc)+ita*fhx[a][b][0];
       gi1=-fey[a][b][1]*cos(c0*dc)+ita*fhx[a][b][1];
       temp=k*((a-10)*dr*sin(c0*dc)+z0*cos(c0*dc));
       gr2=cos(temp);
       gi2=sin(temp);
       gr[a][b][c0]=(gr1*gr2-gi1*gi2);
       gi[a][b][c0]=(gr1*gi2+gi1*gr2);
       Er4[a][b][c0+45]=(-1)*gi[a][b][c0]*(k/(4*pi));
       Ei4[a][b][c0+45]=gr[a][b][c0]*(k/(4*pi));       }
  for(c0=0;c0<=90;c0++)
   for(a=0;a<=X2;a++)
    for(b=0;b<=Y;b++)
     { se4[c0][0]=se4[c0][0]+Er4[a][b][c0];
       se4[c0][1]=se4[c0][1]+Ei4[a][b][c0];}
  for(c0=0;c0<=90;c0++)
   {  Me3[c0]=sqrt(pow(se3[c0][0],2)+pow(se3[c0][1],2));
      Me4[c0]=sqrt(pow(se4[c0][0],2)+pow(se4[c0][1],2)); }
}

void fourier1(float e[rn+1])
{ int n,k;
  float gk,gr[4],gi[4],a[rn+1];
  gk=8.5/df;
  for (k=0;k<=3;k++)
   { for(n=0;n<=rn;n++)
      { a[n]=(2*pi*gk*n)/rn;  }
     gr[k]=0;gi[k]=0;
     for(n=0;n<=rn;n++)
      { gr[k]=gr[k]+e[n]*cos(a[n]);
        gi[k]=gi[k]+e[n]*sin(a[n]);}
     gk=gk+1/df; }
  if (e==Vs) {Vsw[0]=gr[1];
              Vsw[1]=gi[1];}
  if (e==Is) {Isw[0]=gr[1];
              Isw[1]=gi[1];}
}
void mo(float e[X2+1][Y+1][2])
{FILE *fp;
 int i,j;
 float mod[X2+1][Y+1];
 for(i=0;i<=X2;i++)
  for(j=0;j<=Y;j++)
  { mod[i][j]=sqrt(pow(e[i][j][0],2)+pow(e[i][j][1],2)); }
 if (e==fex){ fp=fopen("mex.m","w");
             for(i=0;i<=X2;i++)
              for(j=0;j<=Y;j++)
               { mex[i][j]=mod[i][j];
                 fprintf(fp,"mex(%d,%d)=%.12f;\n",i+1,j+1,mex[i][j]);}
             fprintf(fp,"mesh(mex);");
             fclose(fp);}
 if (e==fey){ fp=fopen("mey.m","w");
             for(i=0;i<=X2;i++)
              for(j=0;j<=Y;j++)
               { mey[i][j]=mod[i][j];
                 fprintf(fp,"mey(%d,%d)=%.12f;\n",i+1,j+1,mey[i][j]);}
             fprintf(fp,"mesh(mey);");
             fclose(fp);}
 if (e==fhx){ fp=fopen("mhx.m","w");
             for(i=0;i<=X2;i++)
              for(j=0;j<=Y;j++)
               { mhx[i][j]=mod[i][j];
                 fprintf(fp,"mhx(%d,%d)=%.12f;\n",i+1,j+1,mhx[i][j]);}
             fprintf(fp,"mesh(mhx);");
             fclose(fp);}
 if (e==fhy){ fp=fopen("mhy.m","w");
             for(i=0;i<=X2;i++)
              for(j=0;j<=Y;j++)
               { mhy[i][j]=mod[i][j];
                 fprintf(fp,"mhy(%d,%d)=%.12f;\n",i+1,j+1,mhy[i][j]);}
             fprintf(fp,"mesh(mhy);");
             fclose(fp);}
}