using Plots, Latexify

ϵ0 = 8.854187813*(10^-12)
μ0 = 1.256637062*(10^-6)
j = im

#questão 1
d2_1 = 35.21997862/1000
#

#questão 2
d2 = 0.049808571430478205
d3 = 0.024904285715239102
#
global f0 = 1.064*(10^9)

function getN(μr,ϵr)
    return √((μr*μ0)/(ϵ0*ϵr))
end

function getk(μr,ϵr, f)
    return 2*pi*f*√(μr*μ0*ϵ0*ϵr)
end

function getλ(f, μr, ϵr)
    return 1/(f*√(μr*μ0*ϵ0*ϵr))
end

function ρN(N1, N0)
    return ((N1-N0)/(N1+N0))    
end

function ρ_phase(z, ρ, k)
    return ρ*exp(-j*2*k*z)
end

function getN_ρ(ρ, N)
    return N*((ρ+1)/(1-ρ))
end

n1 = getN(1,1)
n2 = getN(1,4)
n3 = getN(1,16)
n4 = getN(1,16)

#println("d2 = ", getλ(f0, 1, 2)/4, " || d3 = ", getλ(f0, 1, 8)/4)

ρ3 = ρN(n4,n3)



global fi = f0 - 500*(10^6)

global x = []
global y = []


global bw_y = []



global fc_l = 0
global fc_h = 0
global count_f0 = 1

while (fi <= (f0 + 500*(10^6)))
    #=
    local k3 = getk(1,8,fi)
    local ρ2_p = ρ_phase(d3,ρ3,k3) 
    local n3_p = getN_ρ(ρ2_p, n3)
    =#
    
    local ρ2 = ρN(n3,n2)
    local k2 = getk(1,4,fi)
    local ρ1_p = ρ_phase(d2_1,ρ2,k2)
    local n2_p = getN_ρ(ρ1_p,n2)
    local ρ1 = ρN(n2_p,n1)
    
    local ab2_ρ1 = abs2(ρ1)

    append!(x, fi)
    append!(y, ab2_ρ1)
    
    if (ab2_ρ1 <= 0.001)
        append!(bw_y, ab2_ρ1)
        if(fc_l == 0)
            global fc_l = fi
        end
    else
        if((fc_h == 0) && (fc_l != 0))
            global fc_h = fi
        end
        append!(bw_y, 1)
    end
    
    if (fi != f0) && (fi < f0)
        global count_f0+= 1
    end
    global fi += 10*10^(3)
end

p = plot(x,y, label="", title=((latexify("abs2(ρ_1)"))))
plot!(x, bw_y, fillrange = 1, fillalpha = 0.125, fillcolor = :red, label = "largura da banda", legend = :outerbottom)
xlims!(p, (f0 - 500 *10^6), (f0 + 500*10^6))
ylims!(p, 0, 0.11)
xticks!(p, (f0-500*(10^6)):(150*(10^6)):(f0+500*(10^6)))
yticks!(p,(0:0.01:0.11))
vline!([f0], label="Frequência central(1.064 GHz)", legend = :outerbottom)
println(x[1]," Hz: ", y[1], ";\n",x[length(x)]," Hz: ", y[length(y)])
println(x[count_f0], " Hza: ", y[count_f0])
display(p)

println("fc_l: ", fc_l, "  ||  fc_h: ", fc_h, " || BW: ", (fc_h-fc_l))
savefig(p, "c3")
#println(bw_y)