Displays an FDTD electrodynamics simulation as an animation

FDTD.c

//
//  main.c
//  scruffy4mm
//
//  Created by Frederick Benjamin Woodruff on 19/03/2017.
//  Copyright © 2017 Frederick Benjamin Woodruff. All rights reserved.
//
//
//  main.c
//  5mmWaveplate
//
//  Created by Frederick Benjamin Woodruff on 28/02/2017.
//  Copyright © 2017 Frederick Benjamin Woodruff. All rights reserved.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <pthread.h>

#define DIMS 3
#define NTHREADS 5
#define strM 1023
#define speedC 299792458.0 // m/s
#define VAC_PMBY 4e-7*M_PI
#define VAC_PTVY 1./(VAC_PMBY*speedC*speedC)
#define FLOAT float
#define DOUBLE double // for extension to complex

#define max(a,b) \
({ __typeof__ (a) _a = (a); \
__typeof__ (b) _b = (b); \
_a > _b ? _a : _b; })
#define min(a,b) \
({ __typeof__ (a) _a = (a); \
__typeof__ (b) _b = (b); \
_a < _b ? _a : _b; })


char ExReferenceFilename[strM] = "a4mmDATAscrappy/ExPaperRef.csv"; // free space
char EyReferenceFilename[strM] = "a4mmDATAscrappy/EyPaperRef.csv";
char EzReferenceFilename[strM] = "a4mmDATAscrappy/EzPaperRef.csv";

char ExSampleFilename[strM] = "a4mmDATAscrappy/ExPaperSample.csv"; // pillar array
char EySampleFilename[strM] = "a4mmDATAscrappy/EyPaperSample.csv";
char EzSampleFilename[strM] = "a4mmDATAscrappy/EzPaperSample.csv";

char materialFile[strM] = "a4mmDATAscrappy/PaperMaterial.csv";

//char ExPtReferenceFilename[strM] = "ExPtFieldFSPaperf.csv";
//char EyPtReferenceFilename[strM] = "EyPtFieldFSPaperf.csv";
//char EzPtReferenceFilename[strM] = "EzPtFieldFSPaperf.csv";

//char ExPtSampleFilename[strM] = "ExPtFieldPAPaperf.csv";
//char EyPtSampleFilename[strM] = "EyPtFieldPAPaperf.csv";
//char EzPtSampleFilename[strM] = "EzPtFieldPAPaperf.csv";

char dExReferenceFileName[strM]= "a4mmDATAscrappy/dExPaper4mmReference.csv";
char dEyReferenceFileName[strM]= "a4mmDATAscrappy/dEyPaper4mmReference.csv";
char dEzReferenceFileName[strM]= "a4mmDATAscrappy/dEzPaper4mmReference.csv";

char dHxReferenceFileName[strM]= "a4mmDATAscrappy/dHxPaper4mmReference.csv";
char dHyReferenceFileName[strM]= "a4mmDATAscrappy/dHyPaper4mmReference.csv";
char dHzReferenceFileName[strM]= "a4mmDATAscrappy/dHzPaper4mmReference.csv";

char dExSampleFileName[strM]= "a4mmDATAscrappy/dExPaper4mmSample.csv";
char dEySampleFileName[strM]= "a4mmDATAscrappy/dEyPaper4mmSample.csv";
char dEzSampleFileName[strM]= "a4mmDATAscrappy/dEzPaper4mmSample.csv";

char dHxSampleFileName[strM]= "a4mmDATAscrappy/dHxPaper4mmSample.csv";
char dHySampleFileName[strM]= "a4mmDATAscrappy/dHyPaper4mmSample.csv";
char dHzSampleFileName[strM]= "a4mmDATAscrappy/dHzPaper4mmSample.csv";


// Domain
const DOUBLE CouS = 1.8; // 1/Courant stability factor
const FLOAT deltaS = 60e-6/2;
const DOUBLE domainX = 0.006; // metres
const DOUBLE domainY = 0.006; // metres
const DOUBLE domainZ = 0.006; // metres Should be 0.010, smaller value is an approximation
const DOUBLE timePeriod = 5e-11; // seconds



// Current Source
DOUBLE t0=1.2e-12; // seconds
DOUBLE FWHM = 1e-12; //  seconds (pulse's full width half maximum)
DOUBLE frq = 0.8e12; // Hz
DOUBLE phase = 0.4; // radians
const FLOAT sourceZ = domainZ; // m (z)
const FLOAT sourceY = domainY; // m (y)
const FLOAT sourceY0 = domainY/2; // m (y)
const FLOAT sourceX0 = 0;


// Waveplate
const DOUBLE nPaper = 1.41; // refractive index
const DOUBLE alphaPaper = 1460;
const DOUBLE nAir = 1.00027;
FLOAT thPaper = 120e-6; // thickness of a sheet

//float nVac = 2.4;
FLOAT thVac = 180e-6;

FLOAT width = 0.004; // x-direction width in metres
FLOAT xStart = 0.001+domainX-domainY; // metres



FILE* ExReferenceFile;
FILE* EyReferenceFile;
FILE* EzReferenceFile;

FILE* ExSampleFile;
FILE* EySampleFile;
FILE* EzSampleFile;

FILE* materialFieldOut;

//FILE* ExPtReferenceFile;
//FILE* EyPtReferenceFile;
//FILE* EzPtReferenceFile;

//FILE* ExPtSampleFile;
//FILE* EyPtSampleFile;
//FILE* EzPtSampleFile;

FILE* dExReferenceFile;
FILE* dEyReferenceFile;
FILE* dEzReferenceFile;

FILE* dExSampleFile;
FILE* dEySampleFile;
FILE* dEzSampleFile;

FILE* dHxReferenceFile;
FILE* dHyReferenceFile;
FILE* dHzReferenceFile;

FILE* dHxSampleFile;
FILE* dHySampleFile;
FILE* dHzSampleFile;

const FLOAT deltaT = deltaS/(speedC*CouS); // metres
const int GRIDX = domainX/deltaS;
const int GRIDY = domainY/deltaS;
const int GRIDZ = domainZ/deltaS;
const int GRIDYZ = GRIDY*GRIDZ;


const int MEM = GRIDX * GRIDY * GRIDZ;
int bounds[NTHREADS+1];



const int TSTEPS = timePeriod/deltaT;
DOUBLE sourceJz[TSTEPS];

const DOUBLE sigmaPaper = (alphaPaper*nPaper)/(speedC*VAC_PMBY);
const DOUBLE permPaper = nPaper*nPaper;
const DOUBLE permAir =  nAir*nAir;

// dielectric, conductivity, mu, mcond.
DOUBLE vac[4] = {permAir,0,1,0};
DOUBLE paper[4] = {permPaper,sigmaPaper,1,0};
FLOAT C_vac[4];
FLOAT C_paper[4];


FLOAT getC(char material,char type) {
    if(material=='p') {
        return C_paper[type];
    }
    if (material=='v') {
        return C_vac[type];
    }
    printf("error");
    return 0;
}

FLOAT Ex[MEM];
FLOAT Ey[MEM];
FLOAT Ez[MEM];
FLOAT Hx[MEM];
FLOAT Hy[MEM];
FLOAT Hz[MEM];
char materialType[MEM];


static FLOAT EyAtX0[GRIDY*GRIDZ];
static FLOAT HyAtXmax[GRIDY*GRIDZ];
static FLOAT EzAtX0[GRIDY*GRIDZ];
static FLOAT HzAtXmax[GRIDY*GRIDZ];
static FLOAT ExAtY0[GRIDX*GRIDZ];
static FLOAT HxAtYmax[GRIDX*GRIDZ];
static FLOAT EzAtY0[GRIDX*GRIDZ];
static FLOAT HzAtYmax[GRIDX*GRIDZ];
static FLOAT ExAtZ0[GRIDX*GRIDY];
static FLOAT HxAtZmax[GRIDX*GRIDY];
static FLOAT EyAtZ0[GRIDX*GRIDY];
static FLOAT HyAtZmax[GRIDX*GRIDY];

int IJKf(int,int,int);
DOUBLE current_timestamp();
void writeToFile(FILE*, FLOAT*);
void writePointToFile(FILE*, FLOAT*);
void getBounds();

int IJKf(int i, int j, int k) {
    return i*GRIDY*GRIDZ + j*GRIDZ + k;
}

double current_timestamp() {
    struct timeval te;
    gettimeofday(&te, NULL); // get current time
    double milliseconds = te.tv_sec*1000LL + te.tv_usec/1000; // caculate milliseconds
    // printf("milliseconds: %lld\n", milliseconds);
    return milliseconds;
}

void writeToFile(FILE* file, FLOAT *field) {
    for (int i=GRIDX-GRIDY+1; i<GRIDX-1; i++) {
        for (int j=1; j<GRIDY-1; j++) {
            int IJK = IJKf(i,j,GRIDZ/2);
            fprintf(file, "%.8f; ", field[IJK]);
        }
        fprintf(file, "\n");
    }
    fprintf(file, "\n");
}

void writeMaterialFile(FILE* file, char *field) {
    for (int i=GRIDX-GRIDY+1; i<GRIDX-1; i++) {
        for (int j=1; j<GRIDY-1; j++) {
            int IJK = IJKf(i,j,GRIDZ/2);
            fprintf(file, "%c;", field[IJK]);
        }
        fprintf(file, "\n");
    }
    fprintf(file, "\n");
}

void writePointToFile(FILE* file, FLOAT *field) {
    int IJK = IJKf(GRIDX-2,GRIDY/2,GRIDZ/2);
    fprintf(file, "%.8f; ", field[IJK]);
}

void takeWallToSingleValueAndWriteToFile(FILE* file, FLOAT *field) {
    double total=0;
    for (int j=0; j<GRIDY; j++) {
        for (int k=0; k<GRIDZ; k++) {
            int IJK = IJKf(GRIDX-2,j,k);
            total+=field[IJK];
        }
    }
    fprintf(file, "%.8f; ", total);
}



// Ca, Cb, Da, Db
void getConsts(FLOAT *Cvalues,DOUBLE *material) {
    Cvalues[0] =(1-material[1]*deltaT/(2*material[0]*VAC_PTVY))/(1+material[1]*deltaT/(2*material[0]*VAC_PTVY));
    Cvalues[1] = (deltaT/(material[0]*VAC_PTVY*deltaS))/(1+material[1]*deltaT/(2*material[0]*VAC_PTVY));
    Cvalues[2] = (1-(material[3]*deltaT/(2*material[2]*VAC_PMBY)))/(1+(material[3]*deltaT/(2*material[2]*VAC_PMBY)));
    Cvalues[3] = (deltaT/(material[2]*VAC_PMBY*deltaS))/(1+(material[3]*deltaT/(2*material[2]*VAC_PMBY)));
}



// calculated parameters

// ____________________________________________________
// MEMORY INIT.

int initABC() {
    for (int m = 0; m<GRIDX*GRIDZ; m++) {
        ExAtY0[m]=0;
        HxAtYmax[m]=0;
        EzAtY0[m]=0;
        HzAtYmax[m]=0;
    }
    for (int m = 0; m<GRIDY*GRIDZ; m++) {
        EyAtX0[m]=0;
        HyAtXmax[m]=0;
        EzAtX0[m]=0;
        HzAtXmax[m]=0;
    }
    for (int m = 0; m<GRIDX*GRIDY; m++) {
        ExAtZ0[m]=0;
        HxAtZmax[m]=0;
        EyAtZ0[m]=0;
        HyAtZmax[m]=0;
    }
    return 1;
}
void getBounds() {
    for (int tid=0; tid<NTHREADS+1; tid++) {
        bounds[tid]=tid*(GRIDX-2)/(NTHREADS)+1;
    }
}
int ABC(int X0Interface, int XmaxInterface, int Y0Interface, int YmaxInterface, int ZmaxInterface, int Z0Interface) {
    
    FLOAT abccoef = ((CouS )-1)/((CouS)+1);
    
    if(X0Interface) {
        //back
        for(int j=0;j<(GRIDY-1);j++) {
            for(int k=0;k<GRIDZ;k++) {
                Ey[0*(GRIDY*GRIDZ)+j*GRIDZ+k] = EyAtX0[j*GRIDZ+k]-abccoef*(Ey[1*GRIDY*GRIDZ+j*GRIDZ+k]-Ey[0*(GRIDY*GRIDZ)+j*GRIDZ+k]);
                EyAtX0[j*GRIDZ+k] = Ey[1*(GRIDY*GRIDZ)+j*GRIDZ+k];
            }
            for(int j=0;j<GRIDY;j++) {
                for(int k=0;k<(GRIDZ-1);k++) {
                    Ez[0*(GRIDY*GRIDZ)+j*GRIDZ+k] = EzAtX0[j*GRIDZ+k]-abccoef*(Ez[1*(GRIDY*GRIDZ)+j*GRIDZ+k]-Ez[0*(GRIDY*GRIDZ)+j*GRIDZ+k]);
                    EzAtX0[j*GRIDZ+k] = Ez[1*(GRIDY*GRIDZ)+j*GRIDZ+k];
                }
            }
        }
    }
    
    if(XmaxInterface) {
        //front
        for(int j=0;j<(GRIDY-1);j++) {
            for(int k=0;k<GRIDZ;k++) {
                Hy[(GRIDX-1)*(GRIDY*GRIDZ)+j*GRIDZ+k] = HyAtXmax[j*GRIDZ+k]-abccoef*(Hy[(GRIDX-2)*(GRIDY*GRIDZ)+j*GRIDZ+k]-Hy[(GRIDX-1)*(GRIDY*GRIDZ)+j*GRIDZ+k]);
                HyAtXmax[j*GRIDZ+k] = Hy[(GRIDX-2)*(GRIDY*GRIDZ)+j*GRIDZ+k];
            }
        }
        for(int j=0;j<GRIDY;j++) {
            for(int k=0;k<(GRIDZ-1);k++) {
                Hz[(GRIDX-1)*GRIDYZ+j*GRIDZ+k] = HzAtXmax[j*GRIDZ+k]-abccoef*(Hz[(GRIDX-2)*GRIDYZ+j*GRIDZ+k]-Hz[(GRIDX-1)*(GRIDY*GRIDZ)+j*GRIDZ+k]);
                HzAtXmax[j*GRIDZ+k] = Hz[(GRIDX-2)*(GRIDY*GRIDZ)+j*GRIDZ+k];
            }
        }
    }
    
    if(Y0Interface) {
        //left
        for(int i=0;i<(GRIDX-1);i++) {
            for(int k=0;k<GRIDZ;k++) {
                Ex[i*GRIDYZ+0*GRIDZ+k] = ExAtY0[i*GRIDZ+k]-abccoef*(Ex[i*GRIDYZ+1*GRIDZ+k]-Ex[i*GRIDYZ+0*GRIDZ+k]);
                ExAtY0[i*GRIDZ+k] = Ex[i*GRIDYZ+1*GRIDZ+k];
            }
        }
        for(int i=0;i<GRIDX;i++) {
            for(int k=0;k<(GRIDZ-1);k++) {
                Ez[i*GRIDYZ+0*GRIDZ+k] = EzAtY0[i*GRIDZ+k]-abccoef*(Ez[i*GRIDYZ+1*GRIDZ+k]-Ez[i*GRIDYZ+0*GRIDZ+k]);
                EzAtY0[i*GRIDZ+k] = Ez[i*GRIDYZ+1*GRIDZ+k];
            }
        }
    }
    
    if(YmaxInterface) {
        //right
        for(int i=0;i<(GRIDX-1);i++) {
            for(int k=0;k<GRIDZ;k++) {
                Hx[i*GRIDYZ+(GRIDY-1)*GRIDZ+k] = HxAtYmax[i*GRIDZ+k]-abccoef*(Hx[i*GRIDYZ+(GRIDY-2)*GRIDZ+k]-Hx[i*GRIDYZ+(GRIDY-1)*GRIDZ+k]);
                HxAtYmax[i*GRIDZ+k] = Hx[i*GRIDYZ+(GRIDY-2)*GRIDZ+k];
            }
        }
        for(int i=0;i<GRIDX;i++) {
            for(int k=0;k<(GRIDZ-1);k++) {
                Hz[i*GRIDYZ+(GRIDY-1)*GRIDZ+k] = HzAtYmax[i*GRIDZ+k]-abccoef*(Hz[i*GRIDYZ+(GRIDY-2)*GRIDZ+k]-Hz[i*GRIDYZ+(GRIDY-1)*GRIDZ+k]);
                HzAtYmax[i*GRIDZ+k] = Hz[i*GRIDYZ+(GRIDY-2)*GRIDZ+k];
            }
        }
    }
    
    if(Z0Interface) {
        //bottom
        
        for(int i=0;i<(GRIDX-1);i++) {
            for(int j=0;j<GRIDY;j++) {
                Ex[i*GRIDYZ+j*GRIDZ+0] = ExAtZ0[i*GRIDY+j]-abccoef*(Ex[i*GRIDYZ+j*GRIDZ+1]-Ex[i*GRIDYZ+j*GRIDZ+0]);
                ExAtZ0[i*GRIDY+j] = Ex[i*GRIDYZ+j*GRIDZ+1];
            }
        }
        
        for(int i=0;i<GRIDX;i++) {
            for(int j=0;j<(GRIDY-1);j++) {
                Ey[i*GRIDYZ+j*GRIDZ+0] = EyAtZ0[i*GRIDY+j]-abccoef*(Ey[i*GRIDYZ+j*GRIDZ+1]-Ey[i*GRIDYZ+j*GRIDZ+0]);
                EyAtZ0[i*GRIDY+j] = Ey[i*GRIDYZ+j*GRIDZ+1];
            }
        }
        
        
    }
    
    if(ZmaxInterface) {
        //top
        
        for(int i=0;i<(GRIDX-1);i++) {
            for(int j=0;j<GRIDY;j++) {
                Hx[i*GRIDYZ+j*GRIDZ+GRIDZ-1] = HxAtZmax[i*GRIDY+j]-abccoef*(Hx[i*GRIDYZ+j*GRIDZ+GRIDZ-2]-Hx[i*GRIDYZ+j*GRIDZ+(GRIDZ-1)]);
                HxAtZmax[i*GRIDY+j] = Hx[i*GRIDYZ+j*GRIDZ+(GRIDZ-2)];
            }
        }
        
        for(int i=0;i<GRIDX;i++) {
            for(int j=0;j<(GRIDY-1);j++) {
                Hy[i*GRIDYZ+j*GRIDZ+GRIDZ-1] = HyAtZmax[i*GRIDY+j]-abccoef*(Hy[i*GRIDYZ+j*GRIDZ+GRIDZ-2]-Hy[i*GRIDYZ+j*GRIDZ+(GRIDZ-1)]);
                HyAtZmax[i*GRIDY+j] = Hy[i*GRIDYZ+j*GRIDZ+(GRIDZ-2)];
            }
        }
    }
    
    return 1;
}


void* Ekernel(void *x) {
    int tid;
    tid = *((int *) x);
    
    int start = bounds[tid];
    int end = bounds[tid+1];
    for (int i=start; i<end; i++) {
        for (int j=1; j<GRIDY-1; j++) {
            int IJ = i*GRIDYZ + j*GRIDZ;
            for (int k=1; k<GRIDZ-1; k++) {
                int IJK = IJ + k;
                int Ijk = IJK + GRIDYZ;;
                int iJk = IJK + GRIDZ;
                int ijK = IJK + 1;
                Ez[IJK] = getC(materialType[IJK], 0)*Ez[IJK] + getC(materialType[IJK], 1)*(Hy[Ijk]-Hy[IJK]-Hx[iJk]+Hx[IJK]);
                Ex[IJK] = getC(materialType[IJK], 0)*Ex[IJK] + getC(materialType[IJK], 1)*(Hz[iJk]-Hz[IJK]-Hy[ijK]+Hy[IJK]);
                Ey[IJK] = getC(materialType[IJK], 0)*Ey[IJK] + getC(materialType[IJK], 1)*(Hx[ijK]-Hx[IJK]-Hz[Ijk]+Hz[IJK]);
            }
        }
    }
    return NULL;
}
void* Hkernel(void *x) {
    int tid;
    tid = *((int *) x);
    int start = bounds[tid];
    int end = bounds[tid+1];
    for (int i=start; i<end; i++) {
        for (int j=1; j<GRIDY-1; j++) {
            int IJ = i*GRIDYZ + j*GRIDZ;
            for (int k=1; k<GRIDZ-1; k++) {
                int IJK = IJ + k;
                int iJK = IJK - GRIDYZ;
                int IjK = IJK - GRIDZ;
                int IJk = IJK - 1;
                Hx[IJK] = getC(materialType[IJK], 2)*Hx[IJK] + getC(materialType[IJK], 3)*(Ey[IJK]-Ey[IJk]-Ez[IJK]+Ez[IjK]);
                Hy[IJK] = getC(materialType[IJK], 2)*Hy[IJK] + getC(materialType[IJK], 3)*(Ez[IJK]-Ez[iJK]-Ex[IJK]+Ex[IJk]);
                Hz[IJK] = getC(materialType[IJK], 2)*Hz[IJK] + getC(materialType[IJK], 3)*(Ex[IJK]-Ex[IjK]-Ey[IJK]+Ey[iJK]);
            }
        }
    }
    return NULL;
}

void E_CPU_para() {
    // Calculates E and H fields at the next timestep
    pthread_t threads[NTHREADS];
    int rc,ti;
    int thread_args[NTHREADS];
    
    
    for (ti=0; ti<NTHREADS; ++ti)
    {
        thread_args[ti] = ti;
        rc = pthread_create(&threads[ti], NULL, &Ekernel,(void *)&thread_args[ti]);
    }
    for (ti=0; ti<NTHREADS; ++ti) {
        rc = pthread_join(threads[ti], NULL);
    }
}
void H_CPU_para() {
    // Calculates H fields at the next timestep
    pthread_t threads[NTHREADS];
    int rc,ti;
    int thread_args[NTHREADS];
    
    
    for (ti=0; ti<NTHREADS; ++ti)
    {
        thread_args[ti] = ti;
        rc = pthread_create(&threads[ti], NULL, &Hkernel,(void *)&thread_args[ti]);
    }
    for (ti=0; ti<NTHREADS; ++ti) {
        rc = pthread_join(threads[ti], NULL);
    }
}


void initArraysToZero() {
    for (int ijk=0; ijk<MEM; ijk++) {
        
        materialType[ijk]='v';
        Ex[ijk]=0;
        Ey[ijk]=0;
        Ez[ijk]=0;
        Hx[ijk]=0;
        Hy[ijk]=0;
        Hz[ijk]=0;
    }
}


void makesource() {
    DOUBLE twoCSq = 2*(FWHM/2.35482)*(FWHM/2.35482);
    printf("\nsource waveform:\n");
    for (int n=0; n<TSTEPS; n++) {
        DOUBLE t = n*deltaT;
        DOUBLE tx = (t-t0);
        sourceJz[n] = pow(2,-tx*tx/(twoCSq)) * sin(t*frq*2*M_PI + phase);
        if (t<4e-12 && n%3==0) {
            printf("%f; ",sourceJz[n]);
        }
    }
    printf("\n");
}
void source(int n) {
    DOUBLE source = sourceJz[n];
    
    int starty = max(((sourceY0/domainY)*GRIDX*(1-sourceY/domainY)),0);
    int endy = min(((sourceY0/domainY)*GRIDX*(1+sourceY/domainY)),GRIDY);
    int startz = max((0.5*GRIDZ*(1-sourceZ/domainZ)),0);
    int endz = min((0.5*GRIDZ*(1+sourceZ/domainZ)),GRIDZ);
    
    DOUBLE theta = 45; // degrees in yz relative to +z
    
    DOUBLE sintheta = sin(theta * M_PI/180);
    DOUBLE costheta = cos(theta);
    
    for (int i=(sourceX0/domainX)*GRIDX +1; i<(sourceX0/domainX)*GRIDX +3; i++) {
        for (int j=starty; j<endy; j++) {
            for (int k=startz; k<endz; k++) {
                int IJK = IJKf(i,j,k);
                int rsq = ((GRIDY/2)-j)*((GRIDY/2)-j) + ((GRIDZ/2)-k)*((GRIDZ/2)-k);
                int norm = (GRIDY/2)*(GRIDY/2) + (GRIDZ/2)*(GRIDZ/2);
                double mult=pow(2,-6*rsq/norm);
                Ez[IJK]+=source*costheta*mult;
                Ey[IJK]+=source*sintheta*mult;
            }
        }
    }
}


void makeplanar(double theta_l, double phi_l, double x0_l,double y0_l,double z0_l,double d_l, double thickness) {
    double cth=cos(theta_l*M_PI/180);
    double sth=sin(theta_l*M_PI/180);
    double cph=cos(phi_l*M_PI/180);
    double sph=sin(phi_l*M_PI/180);
    
    double yAng = sth*cph;
    double xAng=sth*sph;
    double zAng=cth;
    
    for (int xi=0; xi<GRIDX; xi++) {
        for (int yi=0; yi<GRIDY; yi++) {
            for (int zi=0; zi<GRIDZ; zi++) {
                double xq= xi*domainX/GRIDX;
                double yq= yi*domainY/GRIDY;
                double zq= zi*domainZ/GRIDZ;
                
                int conditionP = (xq*xAng + yq*yAng + zq*zAng<(d_l/2.)+x0_l*xAng + y0_l*yAng + z0_l*zAng);
                int conditionM = (xq*xAng + yq*yAng + zq*zAng>(-d_l/2.)+x0_l*xAng + y0_l*yAng + z0_l*zAng);
                int conditionLR = (fabs(xq-x0_l)<thickness/2.);
                
                if(conditionP && conditionM && conditionLR) {
                    materialType[IJKf(xi,yi,zi)] = 'p';
                }
            }
        }
    }
    
}

void makeWaveplate() {
    
    makeplanar(90.,0,domainX/2,0,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0003,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0006,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0009,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0012,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0015,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0018,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0021,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0024,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0027,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0030,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0033,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0036,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0039,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0042,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0045,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0048,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0051,domainX/2,0.00012,0.004);
    makeplanar(90.,0,domainX/2,0.0054,domainX/2,0.00012,0.004);
}

//_______Near To Far Field______
FLOAT ExTminusXmax[GRIDYZ];
FLOAT EyTminusXmax[GRIDYZ];
FLOAT EzTminusXmax[GRIDYZ];

FLOAT HxTminusXmax[GRIDYZ];
FLOAT HyTminusXmax[GRIDYZ];
FLOAT HzTminusXmax[GRIDYZ];

FLOAT dEx[TSTEPS+10]={0};
FLOAT dEy[TSTEPS+10]={0};
FLOAT dEz[TSTEPS+10]={0};
FLOAT dHx[TSTEPS+10]={0};
FLOAT dHy[TSTEPS+10]={0};
FLOAT dHz[TSTEPS+10]={0};

void initdEH() {
    for (int n=0; n<TSTEPS+10; n++) {
        dEx[n]=0;
        dEy[n]=0;
        dEz[n]=0;
        dHx[n]=0;
        dHy[n]=0;
        dHz[n]=0;
    }
}

void addToPotentialFields(int n, double l) {
    double r = (domainX/2.)+l;
    double f = ((r*r) - domainX/2.)/(speedC*deltaT);
    double nnFloatE = n+0.5+f;
    int nnE = nnFloatE;
    double aE = nnFloatE-nnE;
    
    double nnFloatH = n+f;
    int nnH = nnFloatH;
    double aH = nnFloatH-nnH;
    
    double prefactor = speedC*deltaT*CouS/(4*M_PI*r);
    for (int j=0; j<GRIDY; j++) {
        for (int k=0; k<GRIDZ; k++) {
            int IJKmax=IJKf(GRIDX-2, j, k);
            int wallJK=IJKf(0, j, k);
            
            dEy[nnE-1]+=(1-aE)*prefactor*(Ey[IJKmax]-EyTminusXmax[wallJK]);
            dEy[nnE] += aE*prefactor*(Ey[IJKmax]-EyTminusXmax[wallJK]);
            dEz[nnE-1]+=(1-aE)*prefactor*(Ez[IJKmax]-EzTminusXmax[wallJK]);
            dEz[nnE] += aE*prefactor*(Ez[IJKmax]-EzTminusXmax[wallJK]);
            dEx[nnE-1]+=(1-aE)*prefactor*(Ex[IJKmax]-ExTminusXmax[wallJK]);
            dEx[nnE] += aE*prefactor*(Ex[IJKmax]-ExTminusXmax[wallJK]);
            
            dHy[nnH]-=(1-aH)*prefactor*(Hy[IJKmax]-HyTminusXmax[wallJK]);
            dHy[nnH-1] -= aH*prefactor*(Hy[IJKmax]-HyTminusXmax[wallJK]);
            dHz[nnH]-=(1-aH)*prefactor*(Hz[IJKmax]-HzTminusXmax[wallJK]);
            dHz[nnH-1] -= aH*prefactor*(Hz[IJKmax]-HzTminusXmax[wallJK]);
            dHx[nnH]-=(1-aH)*prefactor*(Hx[IJKmax]-HxTminusXmax[wallJK]);
            dHx[nnH-1] -= aH*prefactor*(Hx[IJKmax]-HxTminusXmax[wallJK]);
            
            ExTminusXmax[wallJK]=Ex[IJKmax];
            EyTminusXmax[wallJK]=Ey[IJKmax];
            EzTminusXmax[wallJK]=Ez[IJKmax];
            HxTminusXmax[wallJK]=Hx[IJKmax];
            HyTminusXmax[wallJK]=Hy[IJKmax];
            HzTminusXmax[wallJK]=Hz[IJKmax];
        }
    }
}

void saveField(FILE * file, FLOAT * field) {
    for (int n=0; n<TSTEPS; n++) {
        fprintf(file, "%.8f; ", field[n]);
    }
}



int main(int argc, const char * argv[]) {
    double tOld = current_timestamp();
    
    ExReferenceFile = fopen(ExReferenceFilename,"w");
    EyReferenceFile = fopen(EyReferenceFilename,"w");
    EzReferenceFile = fopen(EzReferenceFilename,"w");
    ExSampleFile = fopen(ExSampleFilename,"w");
    EySampleFile = fopen(EySampleFilename,"w");
    EzSampleFile = fopen(EzSampleFilename,"w");
    
    materialFieldOut = fopen(materialFile, "w");
    
    //ExPtReferenceFile = fopen(ExPtReferenceFilename,"w");
    //EyPtReferenceFile = fopen(EyPtReferenceFilename,"w");
    //EzPtReferenceFile = fopen(EzPtReferenceFilename,"w");
    //ExPtSampleFile = fopen(ExPtSampleFilename,"w");
    //EyPtSampleFile = fopen(EyPtSampleFilename,"w");
    //EzPtSampleFile = fopen(EzPtSampleFilename,"w");
    
    dExReferenceFile = fopen(dExReferenceFileName,"w");
    dEyReferenceFile = fopen(dEyReferenceFileName,"w");
    dEzReferenceFile = fopen(dEzReferenceFileName,"w");
    dHxReferenceFile = fopen(dHxReferenceFileName,"w");
    dHyReferenceFile = fopen(dHyReferenceFileName,"w");
    dHzReferenceFile = fopen(dHzReferenceFileName,"w");
    
    dExSampleFile = fopen(dExSampleFileName,"w");
    dEySampleFile = fopen(dEySampleFileName,"w");
    dEzSampleFile = fopen(dEzSampleFileName,"w");
    dHxSampleFile = fopen(dHxSampleFileName,"w");
    dHySampleFile = fopen(dHySampleFileName,"w");
    dHzSampleFile = fopen(dHzSampleFileName,"w");
    
    getBounds();
    getConsts(C_vac, vac);
    getConsts(C_paper, paper);
    
    // REFERENCE RUN
    
    initArraysToZero(); initABC(); makesource(); initdEH();
    
    printf("size of (t; x, y, z) = (%i; %i, %i, %i)\n", TSTEPS,GRIDX,GRIDY,GRIDZ);
    for (int n=0; n<TSTEPS; n++) {
        if(n%10==0) { printf("Reference: %i/%i\n",n/10,TSTEPS/10); }
        source(n);
        E_CPU_para();
        H_CPU_para();
        addToPotentialFields(n, 5*domainX);
        //E_CPU_seq();
        //H_CPU_seq();
        ABC(1,1,1,1,1,1);
        writeToFile(ExReferenceFile,Ex);
        writeToFile(EyReferenceFile,Ey);
        writeToFile(EzReferenceFile,Ez);
        
        //writePointToFile(ExPtReferenceFile,Ex);
        //writePointToFile(EyPtReferenceFile,Ey);
        //writePointToFile(EzPtReferenceFile,Ez);
        
        //takeWallToSingleValueAndWriteToFile(ExPtFSOut,Ex);
        //takeWallToSingleValueAndWriteToFile(EyPtFSOut,Ey);
        //takeWallToSingleValueAndWriteToFile(EzPtFSOut,Ez);
        
        
    }
    
    saveField(dExReferenceFile, dEx);
    saveField(dEyReferenceFile, dEy);
    saveField(dEzReferenceFile, dEz);
    saveField(dHxReferenceFile, dHx);
    saveField(dHyReferenceFile, dHy);
    saveField(dHzReferenceFile, dHz);
    
    // SAMPLE RUN
    
    initArraysToZero(); initABC(); makesource(); initdEH();
    makeWaveplate();
    writeMaterialFile(materialFieldOut, materialType);
    
    
    
    
    for (int n=0; n<TSTEPS; n++) {
        if(n%10==0) { printf("Sample: %i/%i\n",n/10,TSTEPS/10); }
        source(n);
        E_CPU_para();
        H_CPU_para();
        addToPotentialFields(n, 5*domainX);
        //E_CPU_seq();
        //H_CPU_seq();
        ABC(1,1,1,1,1,1);
        
        writeToFile(ExSampleFile,Ex);
        writeToFile(EySampleFile,Ey);
        writeToFile(EzSampleFile,Ez);
        
        
        //writePointToFile(ExPtSampleFile,Ex);
        //writePointToFile(EyPtSampleFile,Ey);
        //writePointToFile(EzPtSampleFile,Ez);
        
        //takeWallToSingleValueAndWriteToFile(ExPtWPOut,Ex);
        //takeWallToSingleValueAndWriteToFile(EyPtWPOut,Ey);
        //takeWallToSingleValueAndWriteToFile(EzPtWPOut,Ez);
    }
    
    saveField(dExSampleFile, dEx);
    saveField(dEySampleFile, dEy);
    saveField(dEzSampleFile, dEz);
    saveField(dHxSampleFile, dHx);
    saveField(dHySampleFile, dHy);
    saveField(dHzSampleFile, dHz);
    
    double tNew = current_timestamp();
    printf("%.2fs\n", (tNew-tOld)*0.001);
    
    return 0;
}

plotAnimation.py

import numpy
from matplotlib import pyplot as pyplot
from matplotlib import animation
import csv
from matplotlib import gridspec
from matplotlib.lines import Line2D
import matplotlib.text as text

# this plots the paper wave plate

ExFieldFSFile = 'Documents/SSProject/Debug/a4mmDATAscrappy/ExPaperRef.csv'
EyFieldFSFile = 'Documents/SSProject/Debug/a4mmDATAscrappy/EyPaperRef.csv'
EzFieldFSFile = 'Documents/SSProject/Debug/a4mmDATAscrappy/EzPaperRef.csv'

ExFieldPAFile = 'Documents/SSProject/Debug/a4mmDATAscrappy/ExPaperSample.csv'
EyFieldPAFile = 'Documents/SSProject/Debug/a4mmDATAscrappy/EyPaperSample.csv'
EzFieldPAFile = 'Documents/SSProject/Debug/a4mmDATAscrappy/EzPaperSample.csv'

materialField = 'Documents/SSProject/Debug/a4mmDATAscrappy/PaperMaterial.csv'

arrowDim = 32
#cutFactor = 13
start = 0
end = -2
hMax=0.6
colorMax=0.5

#maxSteps=100

deltaT = 8e-14
speedC = 299792458.0
deltaS = 60e-6/2

deltaT = deltaS/(speedC*1.8)

domainX = 0.006 # metres
domainY = domainX # metres
domainZ = domainX # metres
timePeriod =5e-11 # seconds

GRIDX = int(domainX/deltaS)-2
GRIDY = int(domainY/deltaS)-2
GRIDZ = int(domainZ/deltaS)-2
TSTEPS = int(timePeriod/deltaT) +1
print GRIDX
print GRIDY
print TSTEPS

cutFactor = GRIDX/arrowDim

def charCsvToNumpy(filename):
    field = numpy.zeros((GRIDX,GRIDY))
    with open(filename, 'rb') as material:
        Reader = csv.reader(material, delimiter=';', quotechar='|')
        rows=[]
        
        for row in Reader:
            thisRow = row[:-1]
            thisRowAsFloat = numpy.zeros(len(thisRow))
            if row!=[]:
                for xi in range(len(thisRow)):
                    #print ('char: %c,,,'%thisRow[xi])
                    if (thisRow[xi]=='v'):
                        thisRowAsFloat[xi]=0
                    else:
                        thisRowAsFloat[xi]=1
                        
                rows.append(thisRowAsFloat)
            else:
                field=numpy.array(rows).astype(numpy.float)
    return field

def csv2numpy(filename):
    field = numpy.zeros((TSTEPS,GRIDX,GRIDY))
    
    with open(filename, 'rb') as material:
        Reader = csv.reader(material, delimiter=';', quotechar='|')
        rows=[]
        tstep=0
        for row in Reader:
            #if tstep>maxSteps:
                #break
            thisRow = row[:-1]
            if row!=[]:
                rows.append(thisRow)
            else:
                field[tstep]=numpy.array(rows).astype(numpy.float)
                tstep+=1
                rows=[]
    return field
    
def cropField(field,cells):
    field=field[:,cells:-cells,cells:-cells]
    
    
# get fields
EyPAField = csv2numpy(EyFieldPAFile)
EzPAField = csv2numpy(EzFieldPAFile)

materialV = charCsvToNumpy(materialField)



imShowFieldPA = numpy.zeros((TSTEPS,GRIDX,GRIDY))
for i in range(TSTEPS):
    imShowFieldPA[i] = EyPAField[i] - materialV*hMax/2 # may need to explicitly for loop




unitScale=1000
domain = [0, domainX*unitScale, 0, domainY*unitScale]

fig = pyplot.figure(figsize =(12,7.2))

pyplot.xlabel('x/mm')
pyplot.ylabel('y/mm')
pyplot.axis(domain)



objects=[]
xvals=numpy.linspace(0,domainX*unitScale,GRIDX)

for tim in xrange(0,TSTEPS,1):
    timeStamp=(tim*deltaT * 1e12)-2.2
    
    im1 = pyplot.imshow(imShowFieldPA[tim,:,:].T, cmap='gist_heat', alpha=1,interpolation = 'nearest',\
    extent=domain,origin='lower',vmin=-hMax,vmax=hMax)

    
    timer = pyplot.text(0.18,3.58,'%.1fps'%timeStamp)
    objects.append([timer,im1])



ani = animation.ArtistAnimation(fig,objects, interval=10,blit=False)
pyplot.show()