import cmath

s11 = cmath.rect(0.2, -60)
s12 = cmath.rect(0.36, 10)
s21 = cmath.rect(5.8, -60)
s22 = cmath.rect(0.2, -60)

#######  Stability check - d<1 , k>1  ######
d = abs((s11*s22 - s12*s21))
k= (1 + (d*d) - abs(s11*s11) - abs(s22*s22))/(2 * abs(s12*s21))
print("d: " ,d)
print("k: " ,k)
if d<1 and k>1 :
    print("Unconditionally Stable")
else:
    print ("Not Stable")


#######  Stability Circles  ######
CL= ((s22 - d * s11.conjugate()).conjugate())/(abs(s22)**2 - abs(d)**2)
RL= abs((s12 * s21)/(abs(s22)**2 - abs(d)**2))
CS= ((s11 - d * s22.conjugate()).conjugate())/(abs(s11)**2 - abs(d)**2)
RS= abs((s12 * s21)/(abs(s11)**2 - abs(d)**2))
print("CL: ", CL)
print("RL: ", RL)
print("CS: ", CS)
print("RS: ", RS)


# if device isnt unilateral and s12 isnt small enough: bilateral , 2 answers for each Gama
B1 = 1 + abs(s11)**2 - abs(s22)**2 - abs(d)**2
B2 = 1 + abs(s22)**2 - abs(s11)**2 - abs(d)**2
C1 = s11 - d * s22.conjugate()
C2 = s22 - d * s11.conjugate()

GamaS1= (B1 + ((B1)**2 - 4*abs(C1)**2))/(2*C1)
print("GamaS1= ", GamaS1)
GamaS2= (B1 - ((B1)**2 - 4*abs(C1)**2))/(2*C1)
print("GamaS2= ", GamaS2)
GamaL1= (B2 + ((B2)**2 - 4*abs(C2)**2))/(2*C2)
print("GamaL1= ", GamaL1)
GamaL2= (B2 - ((B2)**2 - 4*abs(C2)**2))/(2*C2)
print("GamaL2= ", GamaL2)
Gamain1= s11 + (s12 * s21 * GamaL1)/(1-s22 * GamaL1)
print("Gamain1= ", Gamain1)
Gamain2= s11 + (s12 * s21 * GamaL2)/(1-s22 * GamaL2)
print("Gamain2= ", Gamain2)
Gamaout1= s22 + (s12 * s21 * GamaS1)/(1-s11 * GamaS1)
print("Gamaout1= ", Gamaout1)
Gamaout2= s22 + (s12 * s21 * GamaS2)/(1-s11 * GamaS2)
print("Gamaout2= ", Gamaout2)


####### Matching with Gama ########
Z0= 50

Zin1= Z0*(1+Gamain1)/(1-Gamain1)
print("Zin1= ", Zin1)
Zin2= Z0*(1+Gamain2)/(1-Gamain2)
print("Zin2= ", Zin2)

Zout1= Z0*(1+Gamaout1)/(1-Gamaout1)
print("Zout1= ", Zout1)
Zout2= Z0*(1+Gamaout2)/(1-Gamaout2)
print("Zout2= ", Zout2)


######## Power Gains  ########
#gain of Utilateral device  
GTUmax = abs(s21)**2 / ((1 - abs(s11)**2)*(1 - abs(s22)**2))
print("GTUmax: ", GTUmax)

###### Operating Power Gain : Pl / Pin
G = (abs(s21)**2 * (1-abs(GamaL1)**2)) / ((1-abs(Gamain1)**2)* abs(1-s22 * GamaL1)**2)
print("G= ", G)
##### Available Gain : Pavn / Pavs
GA = (abs(s21)**2 * (1-abs(GamaS1)**2)) / ((1-abs(Gamaout1)**2) * abs(1- s11 * GamaS1)**2)
print("GA= ", GA)
##### Trancducer Gain : Pl / Pavs
GT = (abs(s21)**2 * (1-abs(GamaL1)**2) * (1-abs(GamaS1)**2)) / ((abs(1- GamaS1 * Gamain1)**2) * abs(1-s22 * GamaL1)**2)
print("GT= ", GT)
#Maximum Transducing power gain of bilateral device
GTmax = (abs(s21) / abs(s12)) * (k - (k**2 - 1)**0.5)
print("GTmax= ", GTmax)

######## Divided Max Gains ######
Gsmax= 1/(1-abs(s11)**2)
Glmax= 1/(1-abs(s22)**2)
G0= abs(s21)**2
print ("Gsmax= ", Gsmax)
print ("Glmax= ", Glmax)
print ("G0= ", G0)


###### Unilateral figure of Merit , (M = U) ###### GT/GTUmax
M= (abs(s11)*abs(s12)*abs(s21)*abs(s22))/((1-abs(s11)**2)*(1-abs(s22)**2))
print("M: ", M)
print("check Maximum error: ", 1/(1+M)**2 , ", ", 1/(1-M)**2)
if 1/(1+M)**2 < GT/GTUmax and GT/GTUmax < 1/(1-M)**2 :
    print ("Device is Unilateral")
else:
    print("Out of Range, Couldnot be Unilateral")

# if device is Unilateral: [GamaS = s11.conjugate()] , [GamaL = s22.conjugate()]
GamaSmax= s11.conjugate()
print("GamaSmax= ", GamaSmax)
GamaLmax= s22.conjugate()
print("GamaLmax= ", GamaLmax)



######## Matching Circuits #########
Rs= 50 
Xs= 0
RL= 40
XL= 60
'''
# RL > Rs :  Rs -> series -> parallel -> RL
X1= Xs + (Rs*(RL-Rs)+(Rs/RL)*(XL)**2)**0.5
X2= Xs - (Rs*(RL-Rs)+(Rs/RL)*(XL)**2)**0.5
B1= (Rs-RL)/(RL*Xs+Rs*XL-RL*X1)
B2= (Rs-RL)/(RL*Xs+Rs*XL-RL*X2)
print("X1= ", X1) 
print("X2= ", X2) 
print("B1= ", B1) 
print("B2= ", B2) 

'''
# Rs > RL :  Rs -> parallel -> series -> RL
X1= -XL + (RL*(Rs-RL)+(RL/Rs)*(Xs)**2)**0.5
X2= -XL - (RL*(Rs-RL)+(RL/Rs)*(Xs)**2)**0.5
B1= (Rs-RL)/(Rs*XL+RL*Xs-Rs*X1)
B2= (Rs-RL)/(Rs*XL+RL*Xs-Rs*X2)
print("X1= ", X1) 
print("X2= ", X2) 
print("B1= ", B1) 
print("B2= ", B2)

'''
####### Noise Figures #######
F= 2.24
Fmin= 2
Rn= 4
GamaOpt= cmath.rect(0.45, 180)
Z0= 50

N= ((F - Fmin) / (4*Rn/Z0)) * abs((1+GamaOpt)**2)
print("GamaOpt: ", GamaOpt)
print("N: ", N)
print ("Cf: ", GamaOpt/(N+1))
print("Rf: ", (N*(N+1-abs(GamaOpt)**2))**0.5/(N+1))
'''