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


#define KE 200    /*KE is the number of cells to be used */

            
void main ()
{
float dx[KE],ex[KE],hy[KE],ix[KE],sx[KE],sxm2[KE],sxm1[KE],vc,omega;

int n,k,kc,ke,kstart,nsteps;
float ddx,dt,T,epsz,epsilon,sigma;
float t0,spread,pi,freq_in,arg,pulse;
FILE *fp;//, *fopen();
float ex_low_m1=0,ex_low_m2=0,ex_high_m1=0,ex_high_m2=0,ga[KE];

pi=3.14159;
kc=KE/2;                 //center of the problem space
freq_in=1.0e14;         //incident frequency
ddx=3e8/freq_in/10;     //cell size
dt=ddx/6e8;             //time steps
epsz=8.8e-12;
vc=5.7e13;
omega=2.*pi*1.e15;


  /*initialize to free space*/
for (k=0; k<KE;k++)
{

    ex[k]=0.;
	dx[k]=0.;
    hy[k]=0.;
    ga[k]=0.;
	sx[k]=0.;
	sxm2[k]=0.;
	sxm1[k]=0.;
 }
printf("dielectric starts at-->");
scanf("%d",&kstart);

//these parameters specify the input pulse

t0=50.0;
spread=20;
T=0;
nsteps=1;

//main part of the program

while(nsteps>0)
{
printf("nsteps-->");
scanf("%d",&nsteps);
printf("%d \n",nsteps);

for(n=1;n<=nsteps;n++)
{
T=T+1;
//calculate the dx field
for(k=1;k<KE;k++)
{dx[k]=dx[k]+0.5*(hy[k-1]-hy[k]);}

//put a gaussion pulse at the low endsin(2*pi*freq_in*dt*T)*
//pulse=exp(-0.5*(pow((t0-T)/spread,2)))*sin(2*pi*freq_in*dt*T);
pulse=sin(2*pi*freq_in*dt*T);
//printf("%f",pow((t0-T)/spread,2));
dx[3]=dx[3]+pulse;

for(k=kstart;k<=KE-1;k++){
ga[k]=1;}

for(k=0;k<=KE-1;k++)
{
ex[k]=dx[k]-sx[k];
sx[k]=ga[k]*((1+exp(-vc*dt))*sxm1[k]-exp(-vc*dt)*sxm2[k]+(pow(omega,2.)*dt/vc)*(1-exp(-vc*dt))*ex[k]);
sxm2[k]=sxm1[k];
sxm1[k]=sx[k];
//printf("%f",ex[k]);
}

//boundary conditions
ex[0]=ex_low_m2;
ex_low_m2=ex_low_m1;
ex_low_m1=ex[1];

ex[KE-1]=ex_high_m2;
ex_high_m2=ex_high_m1;
ex_high_m1=ex[KE-2];

//calculate the hy field
for(k=0;k<KE-1;k++)
{hy[k]=hy[k]+0.5*(ex[k]-ex[k+1]);}
}



//write the e field out to a file ex
fp=fopen("Ex.txt","w");
for(k=0;k<KE;k++)
{fprintf(fp,"%6.3f \n",ex[k]);}
fclose(fp);
printf("%5.1f\n",T);

//write the e field out to a file ex
fp=fopen("Hy.txt","w");
for(k=0;k<KE;k++)
{fprintf(fp,"%6.3f \n",hy[k]);}
fclose(fp);
printf("%5.1f\n",T);

//write the e field out to a file ex
fp=fopen("Dx.txt","w");
for(k=0;k<KE;k++)
{fprintf(fp,"%6.3f \n",dx[k]);}
fclose(fp);
printf("%5.1f\n",T);

}
}