A.Abdur Rahman Akib
Subscriber

clear;

load("unit_cell.fsp");

 

period = getnamed("::model","period");

 

# get data from sweep

sname= "radius";

radius = getsweepdata(sname,"radius");

height = getsweepdata(sname,"height");

S = getsweepresult(sname,"S");

phase = pinch(angle(S.S21_Gn),1,1); # the dimension of S21_Gn is [frequency,radius,height]

phase = pinch(unwrap(phase)) - phase(1);

height = height(1);

E = getsweepresult(sname,"E");

H = getsweepresult(sname,"H");

wavelength = E.lambda;

matlabsave("EH_and_phase_vs_radius_fdtd",wavelength,height,period,phase,radius,E,H);

 

#############################################################

# interpolate radius and fields on a denser phase data points

#############################################################

phase_interp = linspace(0,2*pi,361); # limit the phase range to 0 - 2*pi

radius_interp = spline(radius,phase,phase_interp); # phase sampled at an interval of 1 degree.

plot(phase_interp*180/pi,radius_interp*1e9,"phase (degree)","radius (nm)","","linewidth=2");

holdon;

plot(phase*180/pi,radius*1e9,"phase (degree)","radius (nm)","","plot type = point,marker type=x, color=red");

legend("radius_interpolated","radius from sweep");

holdoff;

 

E = getsweepresult(sname,"E");

Ex = pinch(E.Ex); Hx = pinch(H.Hx);

Ey = pinch(E.Ey); Hy = pinch(H.Hy);

Ez = pinch(E.Ez); Hz = pinch(H.Hz);

x = E.x;

y = E.y;

z = E.z;

f = E.f;

wavlength = c/f;

 

# interpolate over radius

Ex = interp(Ex,x,y,radius,x,y,radius_interp);

Ey = interp(Ey,x,y,radius,x,y,radius_interp);

Ez = interp(Ez,x,y,radius,x,y,radius_interp);

Hx = interp(Hx,x,y,radius,x,y,radius_interp);

Hy = interp(Hy,x,y,radius,x,y,radius_interp);

Hz = interp(Hz,x,y,radius,x,y,radius_interp);

 

# spatial sampling

######################

do_sampling = true;

ns = 3; # sampling per period

 

if (do_sampling) {

xs = linspace(min(x),max(x),ns);

ys = linspace(min(y),max(y),ns);

Ex = interp(Ex,x,y,radius_interp,xs,ys,radius_interp);

Ey = interp(Ey,x,y,radius_interp,xs,ys,radius_interp);

Ez = interp(Ez,x,y,radius_interp,xs,ys,radius_interp);

Hx = interp(Hx,x,y,radius_interp,xs,ys,radius_interp);

Hy = interp(Hy,x,y,radius_interp,xs,ys,radius_interp);

Hz = interp(Hz,x,y,radius_interp,xs,ys,radius_interp);

 

Ex = reshape(Ex,[ns,ns,1,1,length(phase_interp)]);

Ey = reshape(Ey,[ns,ns,1,1,length(phase_interp)]);

Ez = reshape(Ez,[ns,ns,1,1,length(phase_interp)]);

Hx = reshape(Hx,[ns,ns,1,1,length(phase_interp)]);

Hy = reshape(Hy,[ns,ns,1,1,length(phase_interp)]);

Hz = reshape(Hz,[ns,ns,1,1,length(phase_interp)]);

 

E = rectilineardataset("E",xs,ys,z);

E.addparameter("f",f);

E.addparameter("radius",radius_interp);

E.addattribute("E",Ex,Ey,Ez);

H = rectilineardataset("H",xs,ys,z);

H.addparameter("f",f);

H.addparameter("radius",radius_interp);

H.addattribute("H",Hx,Hy,Hz);

phase = phase_interp;

radius = radius_interp;

}

 

mat_pillar = getnamed("pillar","material");

index_pillar = getindex("pillar",c/wavelength);

mat_sub = getnamed("substrate","material");

index_sub = getnamed("substrate","index");

 

matlabsave("EH_and_phase_vs_radius_interp_fdtd",wavelength,height,period,mat_pillar,index_pillar,mat_sub,index_sub,phase,radius,E,H);

 

visualize(E);