#@DESCRIPTION{
# Testing nonlinear effects in waveguide (TPA and FCA) with coupled TE and TM modes
# Daniel 2010 Fig9
#
# Parameters from: Brian A. Daniel and Govind P. Agrawal, "Vectorial nonlinear propagation in silicon nanowire waveguides: polarization effects," J. Opt. Soc. Am. B 27, 956-965 (2010)
#}

clear;
closeall;

load("NLSE_WGD_Daniel2010_Fig9.icp");
switchtodesign;

##
run;

##
#Get results
myr=getresult("OOSC_2","mode 1/signal");
x=myr.getparameter(myr.getparameter);
Ax=myr.getattribute(myr.getattribute);

myr=getresult("OOSC_2","mode 2/signal");
x=myr.getparameter(myr.getparameter);
Ay=myr.getattribute(myr.getattribute);


m=1;
S0 = abs(Ax)^2 + abs(Ay)^2;
S1=   abs(Ax)^2 - abs(Ay)^2;
S2 = real(2*Ax * conj(Ay)*exp(1i*2*pi*m));
S3 = imag(2*Ax * conj(Ay)*exp(1i*2*pi*m));

loaddata("NLSE_WGD_Daniel2010_ref.dat");    

## Compare results against reference paper and save benchmark
## Data in Daniel2010 Fig9

##Read data from manuscript for comparison
mydata =readdata("NLSE_WGD_Daniel2010_Fig9b_S0.csv");
sz = size(mydata);
numsamples = sz(1);
x_S0 = mydata(1:numsamples(1),1);
y_S0 = mydata(1:numsamples,2);

mydata =readdata("NLSE_WGD_Daniel2010_Fig9_S1.csv");
sz = size(mydata);
numsamples = sz(1);
x_S1 = mydata(1:numsamples(1),1);
y_S1 = mydata(1:numsamples,2);

mydata =readdata("NLSE_WGD_Daniel2010_Fig9_S2.csv");
sz = size(mydata);"NLSE_WGD_Daniel2010_ref.dat";
numsamples = sz(1);
x_S2 = mydata(1:numsamples(1),1);
y_S2 = mydata(1:numsamples,2);

mydata =readdata("NLSE_WGD_Daniel2010_Fig9_S3.csv");
sz = size(mydata);
numsamples = sz(1);
x_S3 = mydata(1:numsamples(1),1);
y_S3 = mydata(1:numsamples,2);


##
# plot results 
tc= ( x(1) + x(length(x)) ) /2.0;
T=(x-tc);
plot(x_S0,y_S0);
holdon;
plot(T*1e12,S0);
setplot("title","S0");
setplot("x min",-30);
setplot("x max",20);
setplot("y min",0);
setplot("y max",0.4);
setplot("x label","time (ps)");
setplot("y label", "S0 (Watt)");
legend("ref[1]_fig9","INTC" );
holdoff;
##
plot(x_S1,y_S1);
holdon;
plot(x_S2,y_S2);
plot(x_S3,y_S3);
plot(T*1e12,S1/S0);
plot(T*1e12,S2/S0);
plot(T*1e12,S3/S0);
setplot("x min",-25);
setplot("x max",17);
setplot("y min",-1.0);
setplot("y max",1.0);
setplot("x label","time (ps)");
setplot("y label", "Normalized Stokes parameters");
legend("S1 - ref[1]_fig9","S2 - ref[1]_fig9","S3 - ref[1]_fig9","S1 - INTC", "S2 - INTC", "S3 - INTC" );
holdoff;