using LinearAlgebra
using Plots

function unit_cell(Mse::Array{Complex{Float64},2},
		Mtl::Array{Complex{Float64},2},
		Msh::Array{Complex{Float64},2})
	# izracuna ABCD parametre ene π celice sukalnika
	# Mse - ABCD matrika serijske impedance
	# Mtl - ABCD matrika linije
	# Msh - ABCD matrika vzporedne impedance (shunt)
	Mse * Mtl * Msh * Mtl * Mse
end

function AtoS(A::Array{Complex{Float64},2})
	# konvertira ABCD matriko v S matriko
	imenovalec = A[1,1] + A[1,2]/50. + A[2,1]*50. + A[2,2]
	S11 = A[1,1] + A[1,2]/50. - A[2,1]*50. - A[2,2]
	S22 = -A[1,1] + A[1,2]/50. - A[2,1]*50. + A[2,2]
	S12 = 2(A[1,1]*A[2,2] - A[1,2]*A[2,1])
	S21 = 2
	[S11 S12; S21 S22]/imenovalec
end

struct varactor
	R::Float64
	L::Float64
	C::Float64
end

function Zd(d::varactor, f::Float64)
	d.R + im*2*pi*f*d.L - im/2/pi/f/d.C
end

function Yin(Z::Complex{Float64}, β::Float64, l::Float64)
	# admitanca l dolge linije zakljucene na Z
	# karakteristicna impedanca linije je 50/√2 = 35Ω
	√2/50. * (50/√2 + im*Z*tan(β*l)) / (Z+im*50/√2*tan(β*l))
end

function Aseries(Z::Complex{Float64})
	[1. Z; 0. 1.]
end

function Atline(β::Float64, l::Float64)
	[cos(β*l) im*50*sin(β*l); im*sin(β*l)/50. cos(β*l)]
end

function Ashunt(Y::Complex{Float64})
	[1. 0.; Y 1.]
end

struct model
	l1::Float64
	ϵ1::Float64
	l2::Float64
	ϵ2::Float64
end

function response1(f::Float64, m::model, d::varactor)
	# izracunaj vse vmesne velicine in vrni ABCD parametre 
	Z = Zd(d, f)
	Mse = Aseries(Z)

	Mtr = Atline(beta(f, m.ϵ1), m.l1)

	Y = Yin(Z, beta(f, m.ϵ2), m.l2)
	Msh = Ashunt(Y)

	unit_cell(Mse,Mtr,Msh)
end

function S2row(S::Array{Complex{Float64},2})
	mag = 10*log10.(abs2.(S))
	ang = angle.(S)
	[mag[1,1] ang[1,1] mag[2,1] ang[2,1] mag[1,2] ang[1,2] mag[2,2] ang[2,2]]
end

function celo_obmocje(f::Array{Float64,1}, m::model, d::varactor)
	n = length(f)
	out = Array{Float64,2}(undef, n, 9)

	for (i, fi) in enumerate(f)
		abcd = response1(fi,m,d)
		S = AtoS(abcd)
		out[i,1] = fi
		out[i,2:end] = S2row(S)
	end
	out
end

function beta(f::Float64, ϵ::Float64)
	2*pi*f*√ϵ/3e8
end

#dioda = varactor(.5, 0.6e-9, 1.3e-12)
#m = model(2e-3, 2.788, 18e-3, 2.937)
#f = collect(range(1e9, 4e9, length=801))
#
#@time test = celo_obmocje(f, m, dioda)
#p = plot(test[:,1], test[:,2],label="S11")
#plot!(p, test[:,1], test[:,4],label="S21")

function cela_simulacija(f, C, l; mapa="./", R=.5, L=.6e-9)
	for l1 in l
		p11 = plot(title="S₁₁ @ l = $l1 mm",xlabel="Frekvenca (GHz)",
							 ylabel="S11 (dB)")
		p21 = plot(title="S₂₁ @ l = $l1 mm",xlabel="Frekvenca (GHz)",
							 ylabel="S21 (dB)")

		m = model(l1*1e-3, 2.8504, 18e-3, 2.992)
		for Cd in C
			d = varactor(R, L, Cd*1e-12) 
			s2p = celo_obmocje(f, m, d)
			plot!(p11, s2p[:,1], s2p[:,2],label="$Cd pF")
			plot!(p21, s2p[:,1], s2p[:,4],label="$Cd pF")
		end
		savefig(p11, "s11-$l1.pdf")
		savefig(p21, "s21-$l1.pdf")
	end
end



f = collect(range(1e9, 4e9, length=801))
l = range(2.,10., step=2.)
C = [8., 4., 2.]

cela_simulacija(f,C,l)