sname = "length1";
length1 = getsweepdata(sname,"length");
width1 = getsweepdata(sname,"width");
Ex_sweep = getsweepdata(sname,"Ex");
Ey_sweep = getsweepdata(sname,"Ey");
n_wavelength = 6;
n_length = 56;
n_width = 57;

Ex = matrix(n_length,n_width,n_wavelength);
Ey = matrix(n_length,n_width,n_wavelength);

for (i=1:n_length){
    for(j=1:n_width){
        for(k=1:n_wavelength){
            Ex(i,j,k) = Ex_sweep(k,j,i);
            Ey(i,j,k) = Ey_sweep(k,j,i);
        }
    }
}

E = getresult("Phase","E");
lambda = E.lambda;
f = E.f;

E_LCP = Ex+1i*Ey;
E_RCP = Ex-1i*Ey;


phi_RCP = angle(E_RCP);
phi_LCP = angle(E_LCP);


for (i=1:n_length){
    for(j=1:n_width){
            phi_LCP(i,j,1:n_wavelength) = unwrap(phi_LCP(i,j,1:1:n_wavelength));
            phi_RCP(i,j,1:1:n_wavelength) = unwrap(phi_RCP(i,j,1:1:n_wavelength));
    }
}

delta_phi_lw1 = phi_LCP(21,6,1:1:n_wavelength);
delta_phi_lw2 = phi_LCP(26,26,1:1:n_wavelength);

delta_phi = phi_LCP(1:56,1:29,n_wavelength)-phi_LCP(1:56,1:29,1);
delta_phi = delta_phi/pi*180;


phi_lambda_max = phi_LCP(1:56,1:57,6);
phi_lambda_max = phi_lambda_max/pi*180;
phi_lambda_max = mod(phi_lambda_max,360);

matlabsave('Unitcell_sname = "length1";
length1 = getsweepdata(sname,"length");
width1 = getsweepdata(sname,"width");
Ex_sweep = getsweepdata(sname,"Ex");
Ey_sweep = getsweepdata(sname,"Ey");
n_wavelength = 6;
n_length = 56;
n_width = 57;

Ex = matrix(n_length,n_width,n_wavelength);
Ey = matrix(n_length,n_width,n_wavelength);

for (i=1:n_length){
    for(j=1:n_width){
        for(k=1:n_wavelength){
            Ex(i,j,k) = Ex_sweep(k,j,i);
            Ey(i,j,k) = Ey_sweep(k,j,i);
        }
    }
}

E = getresult("Phase","E");
lambda = E.lambda;
f = E.f;

E_LCP = Ex+1i*Ey;
E_RCP = Ex-1i*Ey;


phi_RCP = angle(E_RCP);
phi_LCP = angle(E_LCP);


for (i=1:n_length){
    for(j=1:n_width){
            phi_LCP(i,j,1:n_wavelength) = unwrap(phi_LCP(i,j,1:1:n_wavelength));
            phi_RCP(i,j,1:1:n_wavelength) = unwrap(phi_RCP(i,j,1:1:n_wavelength));
    }
}

#Sigle unitcell phase vs frequency
delta_phi_lw1 = phi_LCP(21,6,1:1:n_wavelength);
#plot(f,delta_phi_lw1/pi*180,"f","Phase(°)","L =1.2 W=0.3","linewidth=2");
delta_phi_lw2 = phi_LCP(26,26,1:1:n_wavelength);
#plot(f,delta_phi_lw2/pi*180,"f","Phase(°)","L =1.3 W=0.7","linewidth=2");


# Phase distribution difference between max and min wavelength 
delta_phi = phi_LCP(1:56,1:29,n_wavelength)-phi_LCP(1:56,1:29,1);
delta_phi = delta_phi/pi*180;
#image(length1,width1,delta_phi,"length","width","different phase");


# Phase distribution in max wavelength
phi_lambda_max = phi_LCP(1:56,1:57,6);
phi_lambda_max = phi_lambda_max/pi*180;
phi_lambda_max = mod(phi_lambda_max,360);
image(length1*1e6,width1*1e6,phi_lambda_max,"length(um)","width(um)","max phase");

matlabsave('unitcell_PhaseMap_ExEy',Ex,Ey,phi_LCP,length1,width1,lambda,f);