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);