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

# Copyright 2018 Lumerical Solutions, Inc.

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


filename = "wg_coarsely_sampled.dat"; # Phase data from coarsely sampled component level simulation
filename_GD =  "wg_spar_coarsely_sampled.dat"; # S-Parameter file with group delay option

L = 100e-6;
n_g = 4;
tau_g = L*n_g/c; # group delay evaluated at 193.1 THz

# Read data
data = readdata(filename); 
f = data(:,1); # Frequency
amp = data(:,2); # Amplitude
phase = data(:,3); # Phase
fc = (f(1)+f(end))/2; # Center frequency

# Calculate phase for S-Parameter file
phase_temp1 = phase-2*pi*(f-fc)*tau_g; # Subtract linear (group delay) phase contribution
phase_temp2 = mod(phase_temp1,2*pi); 
phase_GD = unwrap(phase_temp2);


# Write data to file
dataToWrite = [f, amp, phase_GD];
numsamp = length(f);
header1 = "(\"port 2\",\"TE\",1,\"port 1\",1,\"transmission\","+num2str(tau_g)+")";
header2 = "(\"port 1\",\"TE\",1,\"port 2\",1,\"transmission\","+num2str(tau_g)+")";
if(fileexists(filename_GD)) {rm(filename_GD);}
write(filename_GD, header1);
write(filename_GD,"("+num2str(numsamp)+",3)");
write(filename_GD,num2str(dataToWrite));
write(filename_GD, header2);
write(filename_GD,"("+num2str(numsamp)+",3)");
write(filename_GD,num2str(dataToWrite));

# Plot
plot(f/1e12,phase,phase_temp1,phase_temp2,phase_GD);
setplot("x label","Frequency (THz)");
setplot("y label","Phase (rad)");
legend("Original","After group delay phase subtraction", "After modulo 2 pi", "Final");