/*引入總場域，引入均勻平面波*/
#include <stdio.h>
#include <math.h>
float f=1.20352e9;  /*定義頻率2G*/
#define Pi 3.141592653
#define pi 3.141592653
#define c  300000000
#define Np 4000   //fft變換的點數
float dt=1.5766e-12; /*定義時間單位步長1*10^(-8)*/
#define NN 4000/*定義跌代次數*/
const float ee=1/(36*Pi)*0.000000001;
const float u=4*Pi*0.0000001;
float er=10.2;    //介質介電常數
const float sig=0;    //介質導電率
float e=ee;

const float a=34.0/1000;   //天線尺寸(波長數)
const float b=34.0/1000;
const float d=1.0/1000;     //天線厚度
const float A=3*a;      //介質板電尺寸
const float B=3*b;      //介質板電尺寸
const float D=5.0/1000;   //介質板厚度

const int DD=20;       //微帶片xy面剖分段數10*10
 float dx=a/DD;
 float dy=b/DD;
const float dz=1.0/1000;   //常數為介質剖分段數

const int Nx=80;  /*定義剖分網格數*/
const int Ny=80;  /*定義剖分網格數*/
const int Nz=60;  /*定義剖分網格數*/
const int DDD=3;  //外推數據邊界到PML邊界距離


const int kx=5;   //潰源位置(相對天線正中)
 int ky=10;   //潰源位置(相對天線正中)
const float r0=0.635/1000;

#define NNN  10 //定義PML層厚度
float sgmmax=4.47;


 float Ezx[Nx+1][Ny+1][Nz];/*定義全局z向電場*/
 float Exy[Nx][Ny+1][Nz+1];/*定義全局x向電場*/
 float Eyx[Nx+1][Ny][Nz+1];/*定義全局y向電場*/

 float Hzx[Nx][Ny][Nz+1];/*定義全局z向磁場*/
 float Hxy[Nx+1][Ny][Nz];/*定義全局x向磁場*/
 float Hyx[Nx][Ny+1][Nz];/*定義全局y向磁場*/
             /*1代表前一時刻,2代表后一時刻*/

 float Ezy[Nx+1][Ny+1][Nz];/*定義全局z向電場*/
 float Exz[Nx][Ny+1][Nz+1];/*定義全局x向電場*/
 float Eyz[Nx+1][Ny][Nz+1];/*定義全局y向電場*/

 float Hzy[Nx][Ny][Nz+1];/*定義全局z向磁場*/
 float Hxz[Nx+1][Ny][Nz];/*定義全局x向磁場*/
 float Hyz[Nx][Ny+1][Nz];/*定義全局y向磁場*/
             /*1代表前一時刻,2代表后一時刻*/
 //float E1[NN];      /*E場值隨時間變化*/
     
 //float H1[NN];      /*H場值隨時間變化*/

 

float sgmx,sgmy,sgmz,sgmmx,sgmmy,sgmmz;
float canshu1,canshu2,canshu3,canshu4,canshu5,canshu6;
float Iz=0,Uz=0;

 int flag=1;  //定義傳值時間
 void fuzhi();
 int chuanbo1(int M);
 int chuanbo2(int M);
 int bianjie();
 void distill(int M);
 void waitui(int M);
 void zukang();
 void zhubobi();
 void ffft(float *p1,float *p2);

 //int shuchu(int M);
float rou[NN/2];   //遠區電場,駐波比
float si[180],co[180];
 int flag3=1; /*定義傳值時刻*/
 struct bb
 {
   float shi;
   float xu;
 }pp1[Np],pp2[Np],zk[Np/2],fxhs[180];

 struct bb chu(struct bb p,struct bb q);
 struct bb pr[300];

struct aa
 {
   struct bb sita;
   struct bb fai;
 }fxhs1;
struct cc
 {
     float zhengfu;
     float xiangwei;
 };
struct cc Jxu[Nx-2*NNN-2*DDD][Ny-2*NNN-2*DDD+1],Jxd[Nx-2*NNN-2*DDD][Ny-2*NNN-2*DDD+1],Jmxu[Nx-2*NNN-2*DDD][Ny-2*NNN-2*DDD+1],Jmxd[Nx-2*NNN-2*DDD][Ny-2*NNN-2*DDD+1];
struct cc Jxr[Nx-2*NNN-2*DDD][Nz-2*NNN-2*DDD+1],Jxl[Nx-2*NNN-2*DDD][Nz-2*NNN-2*DDD+1],Jmxr[Nx-2*NNN-2*DDD][Nz-2*NNN-2*DDD+1],Jmxl[Nx-2*NNN-2*DDD][Nz-2*NNN-2*DDD+1];
struct cc Jyu[Nx-2*NNN-2*DDD+1][Ny-2*NNN-2*DDD],Jyd[Nx-2*NNN-2*DDD+1][Ny-2*NNN-2*DDD],Jmyu[Nx-2*NNN-2*DDD+1][Ny-2*NNN-2*DDD],Jmyd[Nx-2*NNN-2*DDD+1][Ny-2*NNN-2*DDD];
struct cc Jyf[Ny-2*NNN-2*DDD+1][Nz-2*NNN-2*DDD],Jyb[Ny-2*NNN-2*DDD+1][Nz-2*NNN-2*DDD],Jmyf[Ny-2*NNN-2*DDD+1][Nz-2*NNN-2*DDD],Jmyb[Ny-2*NNN-2*DDD+1][Nz-2*NNN-2*DDD];
struct cc Jzr[Nx-2*NNN-2*DDD][Nz-2*NNN-2*DDD+1],Jzl[Nx-2*NNN-2*DDD][Nz-2*NNN-2*DDD+1],Jmzr[Nx-2*NNN-2*DDD][Nz-2*NNN-2*DDD+1],Jmzl[Nx-2*NNN-2*DDD][Nz-2*NNN-2*DDD+1];
struct cc Jzf[Ny-2*NNN-2*DDD+1][Nz-2*NNN-2*DDD],Jzb[Ny-2*NNN-2*DDD+1][Nz-2*NNN-2*DDD],Jmzf[Ny-2*NNN-2*DDD+1][Nz-2*NNN-2*DDD],Jmzb[Ny-2*NNN-2*DDD+1][Nz-2*NNN-2*DDD];