import java.io.*;
import java.util.*;


public class CircPolarztn{
  public static double Exini( int tim, double x,double phx) {//Ex all x and time t
    return Math.cos(tim-2*Math.PI*x/200.0 +phx ); 
 }
 public static double Eyini( int tim, double x,double phy) {//Ey all x and time t
    return Math.cos(tim-2*Math.PI*x/200.0 +phy ); 
 }
 public static double Hxini( int tim, double x,double phy) {//Hx all x and time t
    return Math.cos(tim-2*Math.PI*x/200.0 +phy+Math.PI); 
 }
 public static double Hyini( int tim, double x,double phx) {//Hy all x and time t
    return Math.cos(tim-2*Math.PI*x/200.0 +phx); 
 }
public static void main(String[] argv) throws IOException,
    FileNotFoundException{
double dx,dt,c,c0;
double phx,phy;
int time,max,i,n,j,k;
phx=0.5*Math.PI;//phases
phy=0.0;//phase difference pi/2
// c= (c0/dz)*dt
c = 0.1;  //Courant stability condition
time = 100; //change this time, stable up to 1700
max = 201;            
double Ex[][] = new double[max+2][2];// Ex and Hy components
double Hy[][] = new double[max+2][2];     
double Ey[][] = new double[max+2][2];// Ey and Hx components
double Hx[][] = new double[max+2][2];     
//Initialize variables
for(i=0;i<=max;i++){
  for (j=0;j<=1;j++){
     Ex[i][j]=0;
     Ey[i][j]=0;
     Hx[i][j]=0;
     Hy[i][j]=0;
   }
}
//Initial values for E and H fields
 for(k = 0;k<max;k++){//x position
   Ex[k][0]=Exini(0,(double)k,phx);
   Ey[k][0]=Eyini(0,(double)k,phy);
   Hx[k][0]=Hxini(0,(double)k,phy);
   Hy[k][0]=Hyini(0,(double)k,phx);   
 }
PrintWriter w =new PrintWriter(new FileOutputStream("Efield.dat"), true);
PrintWriter q =new PrintWriter(new FileOutputStream("Hfield.dat"), true);
for(n = 0;n<time;n++){     
  for(k = 1;k<max-1;k++){//New Ex, Ey components
    Ex[k][1] = Ex[k][0] + c*(Hy[k-1][0]-Hy[k+1][0]);  
    Ey[k][1] = Ey[k][0] + c*(Hx[k+1][0]-Hx[k-1][0]);
    Hx[k][1] = Hx[k][0] + c*(Ey[k+1][0]-Ey[k-1][0]);
    Hy[k][1] = Hy[k][0] + c*(Ex[k-1][0]-Ex[k+1][0]);    
  }     
  Ex[0][1]  =Ex[0][0] + c*(Hy[200-1][0]-Hy[1][0]);    //periodic boundary
  Ex[200][1]=Ex[200][0]+c*(Hy[200-1][0]-Hy[1][0]);    //conditions for first 
  Ey[0][1]  =Ey[0][0]+  c*(Hx[1][0]- Hx[200-1][0]);   //and last points
  Ey[200][1]=Ey[200][0]+c*(Hx[1][0]- Hx[200-1][0]);
  Hx[0][1]  =Hx[0][0]+  c*(Ey[1][0]- Ey[200-1][0]);
  Hx[200][1]=Hx[200][0]+c*(Ey[1][0]- Ey[200-1][0]);
  Hy[0][1]  =Hy[0][0]+  c*(Ex[200-1][0]-Ex[1][0]);
  Hy[200][1]=Hy[200][0]+c*(Ex[200-1][0]-Ex[1][0]); 	
  for (k=0;k<max;k++){         //update fields old=new
    Ex[k][0]=Ex[k][1];
    Ey[k][0]=Ey[k][1];
    Hx[k][0]=Hx[k][1];
    Hy[k][0]=Hy[k][1];
  }         	 
}//time
for(k=0;k<max;k=k+4){	//to plot every 4 points
    w.println(""+0.0+" "+0.0+" "+k+"  "+Ex[k][1]+" "+Ey[k][1]+"  "+0.0);
    q.println(""+0.0+" "+0.0+" "+k+"  "+Hx[k][1]+" "+Hy[k][1]+"  "+0.0);
   }   
 }
}