串行源程序.txt

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

/*串行源程序*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <conio.h>
#include <time.h>
#include <dos.h>

#define pi  3.1415926
#define c   2.998e-1
#define E1  8.85e-3
#define E2  2.58*E1
#define dt  0.1374642e-2
#define f0  10.5
#define dr  0.713809e-3
#define t0  0.15
#define t   0.05
#define df  1/(218*dt)
#define u   1.25e+3
#define n0  0
#define N   258
#define rn  218
#define dc  pi/90
#define ld  c/f0
#define ita 120*pi
#define z0  35*dr
#define s   0.577
#define X1  80
#define Y   100
#define Z1  75
#define X2  20
#define Z2  100
float ex1[X1+1][Y+1][Z1+1],ey1[X1+1][Y+1][Z1+1],ez1[X1+1][Y+1][Z1+1];
float hx1[X1+1][Y+1][Z1+1],hy1[X1+1][Y+1][Z1+1],hz1[X1+1][Y+1][Z1+1];
float ex2[X2+1][Y+1][Z2+1],ey2[X2+1][Y+1][Z2+1],ez2[X2+1][Y+1][Z2+1];
float hx2[X2+1][Y+1][Z2+1],hy2[X2+1][Y+1][Z2+1],hz2[X2+1][Y+1][Z2+1];
float xey[Y+1][Z1+1][10],xez[Y+1][Z1+1][10];
float yexb[X1+1][Z1+1][10],yezb[X1+1][Z1+1][10];
float yexs[X2+1][Z2+2][10],yezs[X2+1][Z2+2][10];
float zex[X1+3][Y+1][10],zey[X1+3][Y+1][10];
float yey1[Y+1][30][14],yez1[Y+1][30][14];
float yex2[Y+1][30][14],yey2[Y+1][30][14];
float zey1[Z1+1][30][14],zez1[Z1+1][30][14];
float zex2[Z1+1][30][14],zez2[Z1+1][30][14];
float xex1[X1+1][30][14],xey1[X1+1][30][14];
float xex2[X1+1][30][14],xez2[X1+1][30][14];
float A,B,C;
float tt0[3][16],tt[8][16];
float rex[X2+1][Y+1][rn+1],rey[X2+1][Y+1][rn+1];
float rhx[X2+1][Y+1][rn+1],rhy[X2+1][Y+1][rn+1];
float fex[X2+1][Y+1][2],fey[X2+1][Y+1][2];
float fhx[X2+1][Y+1][2],fhy[X2+1][Y+1][2];
float Er1[X2+1][Y+1][91],Ei1[X2+1][Y+1][91];
float Er2[X2+1][Y+1][91],Ei2[X2+1][Y+1][91];
float Er3[X2+1][Y+1][91],Ei3[X2+1][Y+1][91];
float Er4[X2+1][Y+1][91],Ei4[X2+1][Y+1][91];
float mex[X2+1][Y+1],mey[X2+1][Y+1];
float mhx[X2+1][Y+1],mhy[X2+1][Y+1];
float se1[91][2],se2[91][2];
float Me1[91],Me2[91];
float se3[91][2],se4[91][2];
float Me3[91],Me4[91];
float Vs[rn+1],Is[rn+1];
float Vsw[2],Isw[2],pin[2],Mp;
float G1[91],G2[91],G3[91],G4[91];
int c7[8]={1,7,21,35,35,21,7,1};


void main()
{  
  FILE  *fp14,*fp15,*fp16,*fp17,*fp18;
  int l,m,m0,m1,m2,m3,i,j,k,i1,j1,k1,n;
  float tmp,e,y,x1,y1,z1;
  time_t first, second;
  void bigareah();
  void smallareah();
  void bigareae();
  void smallareae();
  void murx(int g);
  void mury1(int g);
  void mury2(int g);
  void murz(int g,int xs,int xf,int yf);
  void liaoy(int xx,int zz);
  void liaoz(int xx,int yy);
  void liaox(int xs,int xf,int yy,int zz);
  void fourier(float e[X2+1][Y+1][rn+1]);
  void fourier1(float e[rn+1]);
  void mo(float e[X2+1][Y+1][2]);
  void sss();
  void sss1();
//开始计时
first = time(NULL);

////////////////////////////////////////////////////////////////////////

//初始化定义数组
  

  tt[1][1]=(2-s)*(1-s)/2;
  tt[1][2]=s*(2-s);
  tt[1][3]=s*(s-1)/2;
  for (l=2;l<=7;l++)
   { m2=2*l+1;
     m1=l-1;
     for (j=1;j<=m2;j++)
      {if (j<(2*l-1)) tt0[1][j]=tt[m1][j];
       else tt0[1][j]=0; }
     m3=2*l;
     tt0[2][1]=0;
     for (j=2;j<=m3;j++)
      { j1=j-1;
        tt0[2][j]=tt[m1][j1]; }     tt0[2][m2]=0;
     tt0[3][1]=0;
     tt0[3][2]=0;
     for (j=3;j<=m2;j++)
      { j1=j-2;
        tt0[3][j]=tt[m1][j1]; }
     m0=2*l+1;
     for (i=1;i<=m0;i++)
      { tt[l][i]=0;
        for (j=1;j<=3;j++)
           { tmp=tt[1][j]*tt0[j][i];
             tt[l][i]=tt[l][i]+tmp; }    }   }
  A=(c*dt-dr)/(c*dt+dr);
  B=(2*dr)/(c*dt+dr);
  C=pow((c*dt),2)/(2*dr*(c*dt+dr));
  printf("End1.\n");
  
///////////////////////////////////////////////////////////////////////////
  
  //循环开始
  for(n=1;n<=N;n++)
  { y=(((n*dt-t0)/t)*((n*dt-t0)/t));
    ex1[30][50][13]=1000*exp((-1.0*y))*sin(2*pi*f0*(n-n0)*dt);

/////////////////////////////////////////////////////////////////////
    
    //计算大区域磁场
	bigareah();
   printf("endh1[%d]\n",n);
//////////////////////////////////////////////////////////////////

    ///这里要用到MPI并行
   
    for(i=0;i<=X2;i++)
     for(j=0;j<=Y;j++)
    
        {hx2[i][j][0]=hx1[i+30][j][0];
         hy2[i][j][0]=hy1[i+30][j][0];
         hz2[i][j][0]=hz1[i+30][j][0];}
         
//////////////////////////////////////////////////////////////////////      

   //计算小区域磁场
    smallareah();
    printf("endh2[%d]\n",n);
//////////////////////////////////////////////////////////////////////    
   
    if(n<=rn)
      { //记录Vs,Is
      	Vs[n]=ex1[30][50][13];
       Is[n]=fabs(hz1[30][51][13])+fabs(hz1[30][50][13])+fabs(hy1[30][50][14])+fabs(hy1[30][50][13]); }

///////////////////////////////////////////////////////////////////////////
 
 //计算大区域电场 
bigareae();
printf("ende1[%d]\n",n);
///////////////////////////////////////////////////////////////////////////

//计算小区域电场 
smallareae();
printf("ende2[%d]\n",n);
///////////////////////////////////////////////////////////////////////////
  ///电场传递
    for(i=30;i<=50;i++)
     for(j=0;j<=Y;j++)
        {ex1[i][j][0]=ex2[i-30][j][0];
         ey1[i][j][0]=ey2[i-30][j][0];
         ez1[i][j][0]=ez2[i-30][j][0];}
         

    murx(0); printf("endmur1\n");
    murx(X1);printf("endmur2\n");
    mury1(0);printf("endmur3\n");
    mury2(0);printf("endmur4\n");
    mury1(Y);printf("endmur5\n");
    mury2(Y);printf("endmur6\n");
    murz(Z1,1,X1-1,Y-1);printf("endmur7\n");
    murz(0,1,29,Y-1);printf("endmur8\n");
    murz(0,51,X1-1,Y-1);printf("endmur9\n");
    murz(Z2,1,X2-1,Y-1);printf("endmur10\n");
    liaoy(0,0);printf("endliao1\n");
    liaoy(0,Z1);printf("endliao2\n");
    liaoy(X1,0);printf("endliao3\n");
    liaoy(X1,Z1);printf("endliao4\n");
    liaoz(0,0);printf("endliao5\n");
    liaoz(0,Y);printf("endliao6\n");
    liaoz(X1,0);printf("endliao7\n");
    liaoz(X1,Y);printf("endliao8\n");
    liaox(0,X1,0,Z1);printf("endliao9\n");
    liaox(0,X1,Y,Z1);printf("endliao10\n");
    liaox(0,29,0,0);printf("endliao11\n");
    liaox(0,29,Y,0);printf("endliao12\n");
    liaox(51,X1,0,0);printf("endliao13\n");
    liaox(51,X1,Y,0);printf("endliao14\n");
    liaox(1,X2-1,0,Z2);printf("endliao15\n");
    liaox(1,X2-1,Y,Z2);printf("endliao16\n");

/////////////////////////////////////////////////////////////////////////
    
  
 //////////////////////////////////////////////////////////////////////
 
 //记录口面Ex,Ey,Hx,Hy
         
    if ((n>=36)&&(n<=36+rn) )
 	  {m=n-36;
      for (i=30;i<=50;i++)
 	    for (j=0;j<=Y;j++)
        { rex[i-30][j][m]=ex1[i][j][35];
          rey[i-30][j][m]=ey1[i][j][35];
          rhx[i-30][j][m]=hx1[i][j][35];
          rhy[i-30][j][m]=hy1[i][j][35]; }  }

 }

  
/////////////////////////////////////////////////////////////////////////////

//将记录数组傅立叶变换  

 fourier(rex);printf("endfourier1\n");
 fourier(rey);printf("endfourier2\n");
 fourier(rhx);printf("endfourier3\n");
 fourier(rhy);printf("endfourier4\n");
 
 //////////////////////////////////////////////////////////////
 //积分并求增益
 mo(fex);
 mo(fey);
 mo(fhx);
 mo(fhy);
 
 
 sss();
 sss1();printf("endsss\n");
 fourier1(Vs);printf("endVs.\n");
 fourier1(Is);printf("endIs.\n");
 fp14=fopen("pin.txt","w");
 pin[0]=Vsw[0]*Isw[0]+Vsw[1]*Isw[1];
 pin[1]=-Vsw[0]*Isw[1]+Vsw[1]*Isw[0];
 Mp=sqrt(pow(pin[0],2)+pow(pin[1],2));
 fprintf(fp14,"pin=%.12f+j%.12f\nMode=%.12f\n",pin[0],pin[1],Mp);
 fclose(fp14);
 for(n=0;n<=90;n++)
  { tmp=fabs(pow(Me1[n],2)*pow(dr,2)/(30*Mp));
    G1[n]=10*log10(tmp); }
 for(n=0;n<=90;n++)
  {tmp=fabs(pow(Me2[n],2)*pow(dr,2)/(30*Mp));
   G2[n]=10*log10(tmp); }
 for(n=0;n<=90;n++)
  {tmp=fabs(pow(Me3[n],2)*pow(dr,2)*cos((n0-45)*dc)/(30*Mp));
   G3[n]=10*log10(tmp);}
 for(n=0;n<=90;n++)
  {tmp=fabs(pow(Me4[n],2)*pow(dr,2)/(30*Mp));
   G4[n]=10*log10(tmp);}
 fp15=fopen("G1.m","w");
 for(n=0;n<=90;n++)
  fprintf(fp15,"y(%d)=%.12f;\n",n+1,G1[n]);
 fprintf(fp15,"plot(y);");
 fclose(fp15);
 fp16=fopen("G2.m","w");
 for(n=0;n<=90;n++)
  fprintf(fp16,"y(%d)=%.12f;\n",n+1,G2[n]);
 fprintf(fp16,"plot(y);");
 fclose(fp16);
 fp17=fopen("G3.m","w");
 for(n=0;n<=90;n++)
  fprintf(fp17,"y(%d)=%.12f;\n",n+1,G3[n]);
 fprintf(fp17,"plot(y);");
 fclose(fp17);
 fp18=fopen("G4.m","w");
 for(n=0;n<=90;n++)
  fprintf(fp18,"y(%d)=%.12f;\n",n+1,G4[n]);
 fprintf(fp18,"plot(y);");
 fclose(fp18);
 //结束计时
second = time(NULL);
printf("The running time is: %f seconds",difftime(second,first));
 getch();

 }

	void bigareah()

	{	int i,j,k,i1,j1,k1;
  float x1,y1,z1;


		for(i=0;i<=X1-1;i++)
     for(j=0;j<=Y-1;j++)
      for(k=0;k<=Z1-1;k++)
       {i1=i+1;
        j1=j+1;
        k1=k+1;
        x1=ey1[i][j][k1]-ey1[i][j][k]+ez1[i][j][k]-ez1[i][j1][k];
        hx1[i][j][k]=hx1[i][j][k]+dt/(u*dr)*x1;
        y1=ez1[i1][j][k]-ez1[i][j][k]+ex1[i][j][k]-ex1[i][j][k1];
	     hy1[i][j][k]=hy1[i][j][k]+dt/(u*dr)*y1;
        z1=ex1[i][j1][k]-ex1[i][j][k]+ey1[i][j][k]-ey1[i1][j][k];
        hz1[i][j][k]=hz1[i][j][k]+dt/(u*dr)*z1;    }
	  }

void smallareah()
   {   int i,j,k,i1,j1,k1;
  float x1,y1,z1;
	
	for(i=0;i<=X2-1;i++)
     for(j=0;j<=Y-1;j++)
      for(k=1;k<=Z2;k++)
       {i1=i+1;
        j1=j+1;
        k1=k-1;
        x1=ey2[i][j][k1]-ey2[i][j][k]+ez2[i][j][k]-ez2[i][j1][k];
        hx2[i][j][k]=hx2[i][j][k]+dt/(u*dr)*x1;
        y1=ez2[i1][j][k]-ez2[i][j][k]+ex2[i][j][k]-ex2[i][j][k1];
        hy2[i][j][k]=hy2[i][j][k]+dt/(u*dr)*y1;
        z1=ex2[i][j1][k]-ex2[i][j][k]+ey2[i][j][k]-ey2[i1][j][k];
        hz2[i][j][k]=hz2[i][j][k]+dt/(u*dr)*z1;  }
	  }

void bigareae()
	{ int i,j,k,i1,j1,k1;
  float e,x1,y1,z1;
	
	for(i=1;i<=X1-1;i++)
     for(j=1;j<=Y-1;j++)
      for(k=1;k<=Z1-1;k++)
       {if((i>30 && i<50)&&(j>32 && j<50)&&(k>2 && k<24)) e=E2;
        else e=E1;
        i1=i-1;
        j1=j-1;
        k1=k-1;
        x1=hz1[i][j][k]-hz1[i][j1][k]+hy1[i][j][k1]-hy1[i][j][k];
        ex1[i][j][k]=ex1[i][j][k]+dt/(e*dr)*x1;
        y1=hx1[i][j][k]-hx1[i][j][k1]+hz1[i1][j][k]-hz1[i][j][k];
        ey1[i][j][k]=ey1[i][j][k]+dt/(e*dr)*y1;
        z1=hy1[i][j][k]-hy1[i1][j][k]+hx1[i][j1][k]-hx1[i][j][k];
        ez1[i][j][k]=ez1[i][j][k]+dt/(e*dr)*z1;
        if ((i==30 || i==50) && k<=35 )
          {ey1[i][j][k]=0;
           ez1[i][j][k]=0; }
        if((i>30 && i<50)&&(j>=32 && j<=50)&&(k==2 || k==24))
            {ex1[i][j][k]=(ex1[i][j][k-1]+ex1[i][j][k+1])/2;
             ey1[i][j][k]=(ey1[i][j][k-1]+ey1[i][j][k+1])/2;
             ez1[i][j][k]=(ez1[i][j][k-1]+ez1[i][j][k+1])/2;}
        if((i>30 && i<50)&&(k>=2 && k<=24)&&(j==32 || j==50))
            {ex1[i][j][k]=(ex1[i][j-1][k]+ex1[i][j+1][k])/2;
             ey1[i][j][k]=(ey1[i][j-1][k]+ey1[i][j+1][k])/2;
             ez1[i][j][k]=(ez1[i][j-1][k]+ez1[i][j+1][k])/2;}
        if ((i>=31 && i<=34) && j==50 && k==13)
            { ex1[i][j][k]=0; }
       }
	 }

void smallareae()