#####################################################
# microstrip_rlc_port.lsf
#
# Reflected power for a microstrip transmission line
# with a lumped RLC element 
#
# Copyright 2016 Lumerical Solutions Inc
#####################################################

# Parameters:
Z0 = 40.4236; # characteristic impedance of the microstrip from FDE simulation

R = 1.5*Z0; # resistance
L = 6.5e-9; # inductance 
C = 2e-12; # capacitance 

# Flags to enable the desired components of the RLC element:
# Example: for a pure resistance use: test_R = true; test_L = false; test_C = false;
test_R = true;
test_L = true;
test_C = true;

delta = 0.15; #tolerance for upper and lower bounds.

#-----------------------------

# set up RLC properties
switchtolayout;
setnamed("load","material","<RLC>");
setnamed("load","c enabled",test_C);
setnamed("load","l enabled",test_L);
setnamed("load","r enabled",test_R);

if (test_R){
setnamed("load","r enabled",true);
setnamed("load","r",R);
}

if (test_L){
setnamed("load","l enabled",true);
setnamed("load","l",L);
}

if (test_C){
setnamed("load","c enabled",true);
setnamed("load","c",C);
}


#run simulation with resistor:
save;
run;


#get reflection from mode expansion monitor:
S_result = getresult("FDTD::ports::port 1","S");
reflection = abs(S_result.S)^2;

f = S_result.f;


#analytical calculation: 
j = 0.0 + 1i;
omega = 2.0*pi*f;

ZR_inv = matrix(length(f));
ZL_inv = matrix(length(f));
ZC_inv = matrix(length(f));

ZR_inv_pdelta = matrix(length(f));
ZL_inv_pdelta = matrix(length(f));
ZC_inv_pdelta = matrix(length(f));

ZR_inv_ndelta = matrix(length(f));
ZL_inv_ndelta = matrix(length(f));
ZC_inv_ndelta = matrix(length(f));

if (test_R)
{ZR_inv = ZR_inv + 1.0/R;
dR = delta*R;
ZR_inv_pdelta = ZR_inv_pdelta + 1.0/(R+dR);
ZR_inv_ndelta = ZR_inv_ndelta + 1.0/(R-dR);
}

if (test_L)
{ZL_inv = 1.0/(j*omega*L);
dL = delta*L;
ZL_inv_pdelta = 1.0/(j*omega*(L+dL));
ZL_inv_ndelta = 1.0/(j*omega*(L-dL));
}

if (test_C) 
{ZC_inv = j*omega*C;
dC = delta*C;
ZC_inv_pdelta = j*omega*(C+dC);
ZC_inv_ndelta = j*omega*(C-dC);
}

Zload = Z0 + 1.0/(ZR_inv + ZL_inv + ZC_inv); #Load impedance
Zload_pdelta = Z0 + 1.0/(ZR_inv_pdelta + ZL_inv_pdelta + ZC_inv_pdelta); #Upper boundary of load impedance
Zload_ndelta = Z0 + 1.0/(ZR_inv_ndelta + ZL_inv_ndelta + ZC_inv_ndelta); #Lower boundary of load impedance

Gamma2_load = abs((Zload-Z0)/(Zload+Z0))^2; #Reflected power
Gamma2_load_pdelta = abs((Zload_pdelta-Z0)/(Zload_pdelta+Z0))^2; #Upper boundary for reflected power
Gamma2_load_ndelta = abs((Zload_ndelta-Z0)/(Zload_ndelta+Z0))^2; #Lower boundary for reflected power

#plots:
plot(f*1e-9,reflection,Gamma2_load, Gamma2_load_pdelta, Gamma2_load_ndelta,"f (GHz)","power","Reflected power");
legend("numerical","analytical", "upper bound", "lower bound");

# Smith chart:
ncent=find(f,1e9);
Zcent=Zload(ncent);