/* fd3d_4.3c 3D Radiation from a dipole an FDTD program with
    a senven_point PML and total/scantteried field formulation*/


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

#define IE 100          /*IE is the number of cells to the X lable*/
#define JE 140        /*IE is the number of cells to the Y lable*/
#define KE 100        /*IE is the number of cells to the Z lable*/
#define ia 15
#define ja 15
#define ka 15
#define pi 3.1415926
#define epsz 8.85419e-12
#define TT 500
#define LL 7
main()
{
   static float dx[IE][JE][KE],dy[IE][JE][KE],dz[IE][JE][KE];
   static float ex[IE][JE][KE],ey[IE][JE][KE],ez[IE][JE][KE];
   static float hx[IE][JE][KE],hy[IE][JE][KE],hz[IE][JE][KE];
   static float gax[IE][JE][KE],gay[IE][JE][KE],gaz[IE][JE][KE];
   static int i,j,n,jc,ic,kc,k,NSTEPS,npml,n_pml,jmetal;
   static int ib,jb,kb,kzh,ixh,jyh,T;
   static float xn,xxn,curl_h,curl_e,curl_d;
   static float ddx,dt,epsilon,omega;
   static float t0,spread,pulse;
    FILE *fp,*fopen(),*fptime,*fpt1,*fpt2,*fpt3,*fpt4,*fpt5,*fpt6,*fpt7;
   static float ez_inc[JE],hx_inc[JE];
   static float ez_low_m1=0.0,ez_low_m2=0.0,ez_high_m1=0.0,ez_high_m2=0.0;
   static float gi1[IE],gi2[IE],gi3[IE];
   static float gj1[JE],gj2[JE],gj3[JE];
   static float gk1[KE],gk2[KE],gk3[KE];
   static float idxl[ia][JE][KE],idxh[ia][JE][KE];
   static float ihxl[ia][JE][KE],ihxh[ia][JE][KE];
   static float idyl[IE][ja][KE],idyh[IE][ja][KE];
   static float ihyl[IE][ja][KE],ihyh[IE][ja][KE];
   static float idzl[IE][JE][ka],idzh[IE][JE][ka];
   static float ihzl[IE][JE][ka],ihzh[IE][JE][ka];

   static int ipos[LL]={0,0,0,0,0,0,0},jpos[LL]={2,10,15,40,50,60,70},kpos[LL]={0,0,0,0,0,0,0};
   static int iwid,kwid,l,jgrid_end;
  static float R_dist[IE][KE][LL],r_y[IE][KE][LL],N_shift[IE][KE][LL];
   static float ezml[IE][KE],Eana[LL][TT],del_e[TT];

    ic=IE/2;
    jc=JE/2;
    kc=KE/2;
    ddx=.01;             /* Set the cell size to 0.01 m */
    dt=ddx/(2*3e8);   /* Calculate the time step  */
    ib=IE-ia-10;
    jb=JE-ja-10;
    kb=KE-ka-10;
    jmetal=50;
    iwid=2;
    kwid=4;

    /*Initialize the arrays */



    for(l=0;l<LL;l++){
        for(i=0;i<=TT;i++){
            Eana[l][i]=0.0;
        }
    }


    for(j=0;j<JE;j++){
        ez_inc[j]=0.0;
        hx_inc[j]=0.0;
        for(k=0;k<KE;k++){
         for(i=0;i<IE;i++){
          ex[i][j][k]=0.0;
          ey[i][j][k]=0.0;
          ez[i][j][k]=0.0;
          dx[i][j][k]=0.0;
          hz[i][j][k]=0.0;
          gax[i][j][k]=1.0;
          gay[i][j][k]=1.0;
          gaz[i][j][k]=1.0;
        }
        }
    }
   for(i=0;i<ia;i++){
        for(j=0;j<JE;j++){
            for(k=0;k<KE;k++){
                idxl[i][j][k]=0.0;
                idxh[i][j][k]=0.0;
                ihxl[i][j][k]=0.0;
                ihxh[i][j][k]=0.0;
            }
        }
    }
    for(i=0;i<IE;i++){
        for(j=0;j<ja;j++){
            for(k=0;k<KE;k++){
                idyl[i][j][k]=0.0;
                idyh[i][j][k]=0.0;
                ihyl[i][j][k]=0.0;
                ihyh[i][j][k]=0.0;
            }
        }
    }

    for(i=0;i<IE;i++){
        for(j=0;j<JE;j++){
            for(k=0;k<ka;k++){
                idzl[i][j][k]=0.0;
                idzh[i][j][k]=0.0;
                ihzl[i][j][k]=0.0;
                ihzh[i][j][k]=0.0;
            }
        }
    }

   for(l=0;l<LL;l++){
       for(i=ic-iwid;i<=ic+iwid;i++){
           for(k=kc-kwid;k<=kc+kwid;k++){
               R_dist[i][k][l]=sqrt(pow((ic+ipos[l]-i),2.0)
                    +pow((kc+kpos[l]-k),2.0)+pow((jpos[l]),2.0));
              r_y[i][k][l]=jpos[l]/R_dist[i][k][l];
              N_shift[i][k][l]=2*R_dist[i][k][l];
           }
       }
   }

   /* Bounday Conditions  */
   for(i=0;i<IE;i++){
          gi1[i]=0.0;
          gi2[i]=1.0;
          gi3[i]=1.0;
          fi1[i]=0.0;
          fi2[i]=1.0;
          fi3[i]=1.0;
      }
       for(j=0;j<JE;j++){
          gj1[j]=0.0;
          gj2[j]=1.0;
          gj3[j]=1.0;
          fj1[j]=0.0;
          fj2[j]=1.0;
          fj3[j]=1.0;
       }

       for(k=0;k<KE;k++){
          gk1[k]=0.0;
          gk2[k]=1.0;
          gk3[k]=1.0;
          fk1[k]=0.0;
          fk2[k]=1.0;
          fk3[k]=1.0;
       }

       printf("Number of PML cells--->");
       scanf("%d",&npml);
       printf("%d \n",npml);
       n_pml=npml;

       for(i=0;i<=n_pml;i++){
           xxn  =(npml-i)/npml;
           xn   =0.33*pow(xxn,3.0);
           printf(" %d  %7.4f  %7.4f  \n",i,xxn,xn);
              fi1[i]=xn;
              fi1[IE-i-1]=xn;
              gi2[i]      =1.0/(1.0+xn);
              gi2[IE-1-i] =1.0/(1.0+xn);
              gi3[i]      =(1.0-xn)/(1.0+xn);
              gi3[IE-i-1] =(1.0-xn)/(1.0+xn);
            xxn=(npml-i-0.5)/npml;
            xn =0.33*pow(xxn,3.0);
              gi1[i]=xn;
              gi1[IE-i-2]=xn;
              fi2[i]     =1.0/(1.0+xn);
              fi2[IE-2-i]=1.0/(1.0+xn);
              fi3[i]     =(1.0-xn)/(1.0+xn);
              fi3[IE-2-i] =(1.0-xn)/(1.0+xn);
       }

        for(j=0;j<=n_pml;j++){
           xxn  =(npml-j)/npml;
           xn   =0.33*pow(xxn,3.0);
           printf(" %d  %7.4f  %7.4f  \n",i,xxn,xn);
              fj1[j]=xn;
              fj1[JE-j-1]=xn;
              gj2[j]      =1.0/(1.0+xn);
              gj2[JE-1-j] =1.0/(1.0+xn);
              gj3[j]      =(1.0-xn)/(1.0+xn);
              gj3[JE-j-1] =(1.0-xn)/(1.0+xn);
            xxn=(npml-j-0.5)/npml;
            xn =0.33*pow(xxn,3.0);
              gj1[j]=xn;
              gj1[JE-j-2]=xn;
              fj2[j]     =1.0/(1.0+xn);
              fj2[JE-2-j]=1.0/(1.0+xn);
              fj3[j]     =(1.0-xn)/(1.0+xn);
              fj3[JE-2-j] =(1.0-xn)/(1.0+xn);
       }
        for(k=0;k<=n_pml;k++){
           xxn  =(npml-k)/npml;
           xn   =0.33*pow(xxn,3.0);
           printf(" %d  %7.4f  %7.4f  \n",i,xxn,xn);
              fk1[k]=xn;
              fk1[KE-k-1]=xn;
              gk2[k]      =1.0/(1.0+xn);
              gk2[KE-1-k] =1.0/(1.0+xn);
              gk3[k]      =(1.0-xn)/(1.0+xn);
              gk3[KE-k-1] =(1.0-xn)/(1.0+xn);
            xxn=(npml-k-0.5)/npml;
            xn =0.33*pow(xxn,3.0);
              gk1[k]=xn;
              gk1[KE-k-2]=xn;
              fk2[k]     =1.0/(1.0+xn);
              fk2[KE-2-k]=1.0/(1.0+xn);
              fk3[k]     =(1.0-xn)/(1.0+xn);
              fk3[KE-2-k] =(1.0-xn)/(1.0+xn);
       }

   /* add a metal plane at j=jgrid_end*/
    for(i=0;i<IE;i++){
         for(k=0;k<KE;k++){
             j=jmetal;
             gax[i][j][k]=0.0;
             gay[i][j][k]=0.0;
             gaz[i][j][k]=0.0;
         }
     }
  /* add a hole at ic,jc,kc*/
  for(i=ic-iwid;i<=ic+iwid;i++){
         for(k=kc-kwid;k<=kc+kwid;k++){
             j=jmetal;
             gax[i][j][k]=1.0;
             gay[i][j][k]=1.0;
             gaz[i][j][k]=1.0;
         }
     }
     t0=40;              /*Center of the incident pulse*/
     spread=6.0;          /*Width of the incident pulse*/
     T=0;
     NSTEPS=1;

   fptime=fopen("time.dat","w");
   fpt1=fopen("hou2.dat","w");
   fpt2=fopen("hou10.dat","w");
   fpt3=fopen("yuanchang15.dat","w");
   fpt4=fopen("yuanchang40.dat","w");
   fpt5=fopen("yuanchang50.dat","w");
   fpt6=fopen("yuanchang60.dat","w");
   fpt7=fopen("yuanchang70.dat","w");

     while (NSTEPS>0) {
       printf("NSTEPS-->");  /*NSTEPS is the number of times the*/
       scanf("%d",&NSTEPS);  /* main loop has executed*/
       printf("%d\n",NSTEPS);

     for (n=1; n<=NSTEPS; n++)
    {
      T=T+1;           /*T keeps track of the total number*/
      /* Start of the main FDTD Loop    */

   /* Calculate the incident buffer*/
   for(j=1;j<JE;j++){
            ez_inc[j]=ez_inc[j]+0.5*(hx_inc[j-1]-hx_inc[j]);
       }
     ez_inc[jmetal]=0.0;

      for(i=ic-iwid; i<=ic+iwid; i++){
         for(k=kc-kwid; k<=kc+kwid; k++){
             ezml[i][k]=ez[i][jmetal][k];
         }
     }
      /*  Source */

      pulse = sin(2*pi*1500*1e6*dt*T);
     pulse = exp(-0.5*(pow((t0-T)/spread,2.0)));
      ez_inc[3]=ez_inc[3]+pulse;

         /* ABC for the incident buffer*/
       ez_inc[0]       =ez_low_m2;
       ez_low_m2   =ez_low_m1;
       ez_low_m1   =ez_inc[1];

       ez_inc[JE-1]    =ez_high_m2;
       ez_high_m2  =ez_high_m1;
       ez_high_m1  =ez_inc[JE-2];

 /*   Calculate the Dx field  */
      for(i=1;i<ia;i++){
          for(j=jmetal+1;j<JE;j++){
            for(k=1;k<KE;k++){
                 curl_h=(hx[i][j][k]-hx[i][j-1][k]-hy[i][j][k]+hy[i][j][k-1]);
                  idxl[i][j][k]=idxl[i][j][k]+curl_h;
             dx[i][j][k] = dx[i][j][k]+gj2[j]*gk2[k]*(curl_h+
                           gi1[i]*idxl[i][j][k]);

            }
          }
      }

    for(i=ia;i<=ib;i++){
          for(j=jmetal+1;j<JE;j++){
            for(k=1;k<KE;k++){
                 curl_h=(hz[i][j][k]-hz[i][j-1][k]-hy[i][j][k]+hy[i][j][k-1]);
             dx[i][j][k] = gj3[j]dx[i][j][k]+gj2[j]*gk2[k]*0.5*curl_h;
            }
          }
      }



      /* Calculate the Dy field */
          for(i=1;i<IE;i++){
          for(j=1;j<ja;j++){
            for(k=1;k<KE;k++){
                 curl_h=(hx[i][j][k]-hx[i][j][k-1]-hz[i][j][k]+hz[i-1][j][k]);
                  idyl[i][j][k]=idyl[i][j][k]+curl_h;
             dy[i][j][k] = gi3[i]*gk3[k]*dy[i][j][k]+gi2[i]*gk2[k]*0.5*(curl_h+
                           gj1[j]*idyl[i][j][k]);

            }
          }
      }

    for(i=1;i<IE;i++){
          for(j=jmetal+1;j<=jb;j++){
            for(k=1;k<KE;k++){
                 curl_h=(hx[i][j][k]-hx[i][j][k-1]-hz[i][j][k]+hz[i-1][j][k]);
             dy[i][j][k] = gi3[j]*gk3[j]*dy[i][j][j]+gi2[i]*gk2[k]*0.5*curl_h;
            }
          }
      }

      for(i=1;i<IE;i++){
          for(j=jb+1;j<JE;j++){
              jyh=j-jb-1;
            for(k=1;k<KE;k++){
                 curl_h=(hx[i][j][k]-hx[i][j][k-1]-hz[i][j][k]+hz[i-1][j][k]);
                  idyh[i][jyh][k]=idyh[i][jyh][k]+curl_h;
             dy[i][j][k] = gi3[i]*gk3[k]*dy[i][j][k]+gi2[i]*gk2[k]*0.5*(curl_h+
             gj1[j]*idyh[i][jyh][k]);
            }
          }
      }

     /* Calculate the Incident Dy */
     for(i=ia;i<=ib;i++){
          for(j=ja;j<=jb-1;j++){
              dy[i][j][ka]=dy[i][j][ka]-0.5*hx_inc[j];

          }
     }
            /* Calculate the Dz field */
         for(i=1;i<IE;i++){
          for(j=1;j<JE;j++){
            for(k=0;k<ka;k++){
                 curl_h=(hy[i][j][k]-hy[i-1][j][k]-hx[i][j][k]+hx[i][j-1][k]);
                  idzl[i][j][k]=idzl[i][j][k]+curl_h;
             dz[i][j][k] = gi3[i]*gj3[j]*dz[i][j][k]+gi2[i]*gj2[j]*0.5*(curl_h+
                           gk1[k]*idzl[i][j][k]);

            }
          }
      }


      for(i=1;i<IE;i++){
          for(j=1;j<JE;j++){