#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "mpi.h"

#define IE 1000
#define JE 1000

int main(int argc, char *argv[] )
{
 int xd;
 int x;
 int p;
 /*int status;*/

 float cc=2.9979e8;
 float pi=3.14159;
 float muz=4.0*pi*1.0e-7;
 float epsz=1.0/(cc*cc*muz);
 float freq=5.0e9;
 float lambda=cc/freq;
 float omega=2.0*pi*freq;

 float dx=1.0e-8;
 float dt=dx/(2.0*cc);
 int nsteps=10000;
 int w=200000;

 int i;
 int j;
 int n;

 float pulse;

 int IS=IE/2;
 int JS=150;

 float CA[IE][JE],CB[IE][JE],CP[IE][JE],CQ[IE][JE];
 float EX[IE][JE-1],EY[IE-1][JE],HZ[IE][JE];

 float sigma=0.0;
 float sigmam=0.0;

 
 FILE *fp1,*fp2;        /*定義文件指針*/

 float eaf;
 float eafm;

 float LD,RD,RU,LU;

 float bjc=(cc*dt-dx)/(cc*dt+dx);
 float jbjc=(cc*dt-1.4142*dx)/(cc*dt+1.4142*dx);
 float bjc2=(cc*cc*(8.85e-12)*dt)/(2*cc*dt+2*dx);
 float jbjc2=(cc*dt)/(1.4142*dx);
 
 float ZCHZ[IE][JE];     /*用來暫存上一時間步的磁場值*/

 float CDHZ[JE];         /*用來暫存傳遞邊界上的HZ和EY的值*/
 float CDEY[JE];

 eaf=(sigma*dt)/(2*epsz);
 eafm=(sigmam*dt)/(2*muz);

 for(i=0;i<IE;i++)
 {
  for(i=0;j<JE;j++)
  {
   CA[i][j]=(1-eaf)/(1+eaf);
   CB[i][j]=(dt/epsz)/(1+eaf);
   CP[i][j]=(1-eafm)/(1+eafm);
   CQ[i][j]=(dt/muz)/(1+eafm);

   EX[i][j]=0.0;
   EY[i][j]=0.0;
   HZ[i][j]=0.0;
  }
 }

 for(j=0;j<JE;j++)
 {
  CDHZ[j]=0.0;
  CDEY[j]=0.0;
 }

  
 MPI_Init(&argc,&argv);                 /*初始化MPI*/
 MPI_Comm_rank(MPI_COMM_WORLD,&xd);     /*得到各個進程號*/
 MPI_Comm_size(MPI_COMM_WORLD,&p);      /*得到總的進程數*/
 MPI_Status status;

 /*主循環開始*/

 for(n=0;n<10000;n++)
 {
  if(xd==0)           /*0號進程*/
  {
   /*接收進程1的EY的值*/
   MPI_Recv(&CDEY[0],1000,MPI_FLOAT,1,j,MPI_COMM_WORLD,&status);
   
   for(i=1;i<IE/2;i++)
   {
	for(j=1;j<JE-1;j++)
	{
	 ZCHZ[i][j]=HZ[i][j];
    }
   }
   LD=HZ[1][1];
   LU=HZ[1][JE-2];
   
   for(i=1;i<IE/2;i++)
   {
	for(j=1;j<JE-1;j++)
	{
	 HZ[i][j]=CP[i][j]*HZ[i][j]-(CQ[i][j]/dx)*(EY[i][j]-EY[i-1][j]-EX[i][j]+EX[i][j-1]);
    }
   }

   for(j=1;j<JE-1;j++)
   {
	HZ[0][j]=ZCHZ[i][j]+bjc*(HZ[1][j]-HZ[0][j])+bjc2*(EX[0][j]-EX[0][j-1]+EX[1][j]-EX[1][j-1]);
   }
   for(i=1;i<IE/2;i++)
   {
	HZ[i][0]=ZCHZ[i][1]+bjc*(HZ[i][1]-HZ[i][0])-bjc2*(EY[i][0]-EY[i-1][0]+EY[i][1]-EY[i-1][1]);
	HZ[i][JE-1]=ZCHZ[i][JE-1]+bjc*(HZ[i][JE-2]-HZ[i][JE-1])-bjc2*(EY[i][JE-1]-EY[i-1][JE-1]+EY[i][JE-2]-EY[i-1][JE-2]);
   }
   HZ[0][0]=LD+jbjc*(HZ[1][1]-HZ[0][0]);
   HZ[0][JE-1]=LU+jbjc*(HZ[1][JE-2]-HZ[1][JE-1]);

   /*加源*/
   pulse=sin((2*pi*cc/632.8e-9)*n*dt);
   HZ[IS][JS]=pulse;

   for(i=0;i<IE/2;i++)
   {
	for(j=0;j<JE-1;j++)
	{
	 EX[i][j]=CA[i][j]*EX[i][j]+(CB[i][j]/dx)*(HZ[i][j+1]-HZ[i][j]);
	}
   }
   for(i=0;i<IE/2-1;i++)
   {
	for(j=0;j<JE;j++)
	{
	 EY[i][j]=CA[i][j]-(CB[i][j]/dx)/(HZ[i+1][j]-HZ[i][j]);
	}
   }

   for(j=0;j<JE;j++)
   {
	CDHZ[j]=HZ[IE/2-1][j];
   }
   
   /*傳遞HZ的值到進程1*/
   MPI_Send(&CDHZ[0],1000,MPI_FLOAT,1,j,MPI_COMM_WORLD);

   if(n==7000)
   {
	if((fp1=fopen("/public/home/zhouqing/c1","w"))==NULL)
	{
	 printf("cannot opening file\n");
	 exit(1);
	}
   }
   for(i=0;i<IE/2;i++)
   {
    for(j=0;j<JE;j++)
	{
	 printf(fp1,"%f\n",HZ[i][j]);
	}
   }
   fclose(fp1);
  }
 

 else if(xd==1)
  {
   
   for(j=0;j<JE;j++)
   {
   MPI_Recv(&CDHZ[0],1000,MPI_FLOAT,0,j,MPI_COMM_WORLD,&status);
   MPI_Send(&CDEY[0],1000,MPI_FLOAT,0,j,MPI_COMM_WORLD);
   }
   for(i=IE/2;i<IE;i++)
   {
	for(j=0;j<JE;j++)
	{
	 ZCHZ[i][j]=HZ[i][j];
	}
   }
   RD=HZ[IE-2][1];
   RU=HZ[IE-2][JE-2];

   for(i=IE/2;i<IE;i++)
   {
	for(j=0;j<JE;j++)
	{
	 HZ[i][j]=CP[i][j]*HZ[i][j]-(CQ[i][j]/dx)*(EY[i][j]-EY[i-1][j]-EX[i][j]+EX[i][j-1]);
	}
   }
   for(j=1;j<JE-1;j++)
   {
	HZ[IE-1][j]=ZCHZ[IE-2][j]+bjc*(HZ[IE-2][j]-HZ[IE-1][j])+bjc2*(EX[IE-1][j]-EX[IE-1][j-1]+EX[IE-2][j]-EX[IE-2][j-1]);
   }
   for(i=IE/2;i<IE-1;i++)
   {
    HZ[i][0]=ZCHZ[i][1]+bjc*(HZ[i][1]-HZ[i][0])-bjc2*(EY[i][0]-EY[i-1][0]+EY[i][1]-EY[i-1][1]);
	HZ[i][JE-1]=ZCHZ[i][JE-1]+bjc*(HZ[i][JE-2]-HZ[i][JE-1])-bjc2*(EY[i][JE-1]-EY[i-1][JE-1]+EY[i][JE-2]-EY[i-1][JE-2]);
   }
   HZ[IE-1][0]=RD+jbjc*(HZ[IE-2][1]-HZ[IE-1][0]);
   HZ[IE-1][JE-1]=RU+jbjc*(HZ[IE-2][JE-2]-HZ[IE-1][JE-1]);
   
   for(i=IE/2;i<IE;i++)
   {
	for(j=0;j<JE-1;j++)
	{
	 EX[i][j]=CA[i][j]*EX[i][j]+(CB[i][j]/dx)*(HZ[i][j+1]-HZ[i][j]);
	}
   }
   for(i=IE/2-1;i<IE-1;i++)
   {
	for(j=0;j<JE;j++)
	{
	 EY[i][j]=CA[i][j]*EY[i][j]-(CB[i][j]/dx)*(HZ[i+1][j]-HZ[i][j]);
	}
   }

   for(j=0;j<JE;j++)
   {
	CDEY[j]=EY[IE/2-1][j];
   }

   if(n==7000)
   {
	if((fp2=fopen("/public/home/zhouqing/c2","w"))==NULL)
	{
	 printf("cannot opening file\n");
	 exit(1);
	}
   }
   for(i=IE/2;i<IE;i++)
   {
    for(j=0;j<JE;j++)
	{
	 printf(fp2,"%f\n",HZ[i][j]);
	}
   }
   fclose(fp2);
 }
 }

 MPI_Finalize();
 return 0;
}