engee.clear() # clearing the workspace
is_dinamic_plot = true # choosing the graph type:true (plotlyjs()) - interactive, false (gr()) - static
is_dinamic_plot ? plotlyjs() : gr()

let
    installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
    list_packages = ["JLD2","Hungarian"]
    for pack in list_packages
        pack in installed_packages || Pkg.add(pack)
    end
end

using LinearAlgebra,Hungarian,DataFrames
import JLD2

include("utils/plot_sig_and_spec.jl") # the function of building an oscillogram and a spectrogram
include("utils/viz_scene.jl") # the function of building a radar operation scenario
include("utils/calc_plot_rcs.jl") # backscatter diagram construction function
include("utils/sort_ang_mse.jl") # the function of sorting bearing estimates
include("utils/plot_spec_azel.jl") # the function of constructing a spatial spectrum
include("utils/calc_table_ang_analysis.jl") # table calculation function for accuracy analysis

# Fundamental radar parameters
fc = 1e9 # [Hz], carrier frequency
c = physconst("LightSpeed") # [m/s], the propagation speed of the signal
lambda = c / fc; # wavelength, m

println("Carrier frequency: $(FC*1E-9)GHz")
println("Speed of light: $c m/s")
println("Wavelength: $lambda m")

# --- Signal parameters ---
bw = 100e6; # LFM band, Hz
T = 10e-6; # Pulse duration, s
Q = 2 # pulse continuity
prf = round(1 / (Q * T);sigdigits = 15); # Pulse repetition rate (rounding if necessary)
type_dir = "Up" # frequency change direction ["Up","Down"]
type_sweep = "Positive"  # frequency change type ["Positive","Symmetrical"]
n_pulses = 2 # number of accumulation pulses
n_groups_imps = 2 # number of pulse bursts

fs = 2*bw; # Sampling rate, Hz
NBand = 64; # number of frequency sub-bands

T_coh = n_pulses/prf # the time of the coregent accumulation
ΔR = c/(2*bw) # range resolution, m
ΔV = lambda/(2*T_coh) # speed resolution, m/s
Rmin = c*T/2 # minimum one-digit range, m
Rmax = c/(2*prf) # maximum one-digit range, m
Vmax = lambda*prf/4 # maximum single-digit speed, m/s

# formation of the system object of the LFM signal generator
waveform = EngeePhased.LinearFMWaveform( 
    SampleRate = fs, # Sampling rate, Hz
    PulseWidth = T, # Pulse duration
    PRF = prf, # Pulse repetition rate
    SweepBandwidth = bw,  # The LCHM band
    SweepDirection = "Up", # direction of frequency change
    SweepInterval = "Positive", # saw type, non-symmetrical or symmetrical
    OutputFormat = "Pulses", # output format type []
    NumPulses = n_pulses # number of pulses
);

# calculation of the LFM signal
sig_lfm = waveform();

# construction of an oscillogram and a spectrogram
plot_sig_and_spec(sig_lfm;fs=fs,name_sig = "LCHM");

println("Sampling rate: $(round(fs*1e-6;sigdigits=5)) MHz")
println("Pulse duration: $(round(T*1e6;sigdigits=5)) microseconds")
println("Follow-up period: $(round(1e6/prf;sigdigits=5)) mks")
println("Repetition rate: $(round(prf*1e-3;sigdigits=5)) kHz")
println("Borehole: $(round(Q;sigdigits=5))")
println("Coherent accumulation time: $(round(T_coh*1e3;sigdigits=5)) ms")
println("Number of pulses per packet: $n_pulses")
println("Number of packs: $n_groups_imps")
println("Range resolution : $(round(ΔR;sigdigits=3)) m") 
println("Speed resolution : $(round(ΔV;sigdigits=3)) m/s") 
println("Maximum one-digit range : $(round(Rmax;sigdigits=3)) m") 
println("Minimum one-digit range : $(round(Rmin;sigdigits=3)) m") 
println("Maximum single-digit speed: $(round(Vmax;sigdigits=3)) m/s") 

# antenna parameters
d = lambda /2 # the distance between the elements
freq_range = [fc-bw fc+bw] # frequency range of radiation
Nx = 24 # Number of elements on the X-axis (azimuth resolution)
Ny = 24 # Number of elements on the Y-axis (resolution by angle of location)
Δaz = rad2deg(lambda/(d*(Nx-1))) # azimuth resolution, deg
Δel = rad2deg(lambda/(d*(Ny-1))) # seat angle resolution, deg
# Total number of elements N = Nx * Ny

# WITH an isotropic antenna element
ant_elem = EngeePhased.IsotropicAntennaElement( 
    FrequencyRange = freq_range, # frequency range
    BackBaffled = true # removing the backscatter of the antenna bottom
) 

# creating a rectangular headlight
array = EngeePhased.URA(
    Element = ant_elem, # antenna element
    Size = [Nx Ny], # grid size
    ElementSpacing = [d d], # the distance between the elements
);

pattern(array,fc) |> display;

println("Distance between elements : $(round(d;sigdigits=3)) m") 
println("Rectangular HEADLIGHT Dimensions: $(Nx)x$Ny [Nx, Ny]") 
println("Frequency range of radiation: $(freq_range*1e-9) GHz")
println("Azimuth resolution: $(round(Δaz;sigdigits=3)) °") 
println("Location angle resolution: $(round(Δel;sigdigits=3)) °") 

# Parameters of the receiving and transmitting path
Pt = 1e3 # W, peak power
Gtx = 20 # dB, gain of the transmitting antenna
LFtx = 1.5 # dB, losses in the transmission path

Grx = 20 # dB, gain of the receiving antenna
NF = 2 # dB, the noise coefficient of the receiving path
LFrx = 1.5 # dB, losses in the receiving path
TN = 290 # Noise temperature
NP = physconst("boltzmann") * bw * TN; # W, noise power

# Creatures with a transmitter
transmitter = EngeePhased.Transmitter( 
    Gain = Gtx, # transmitter gain
    PeakPower = Pt, # Peak power
    LossFactor=LFtx # losses in the transmission path
)

# Devices WITH broadband transmission lights
radiator = EngeePhased.WidebandRadiator(
    Sensor = array, # The bottom of the antenna element
    SampleRate = fs, # sampling rate
    CarrierFrequency = fc, # carrier frequency
    NumSubbands = NBand # number of frequency sub-bands
)

# devices WITH broadband receiving headlights
collector = EngeePhased.WidebandCollector(
    Sensor = array, # The bottom of the antenna element
    SampleRate = fs, # sampling rate
    CarrierFrequency = fc, # central frequency
    NumSubbands = NBand # number of frequency sub-bands
)

# creation of a receiving preamplifier (LNA)
receiver = EngeePhased.ReceiverPreamp(
    NoiseMethod="Noise power", # noise setting method ["Noise power", "Noise temperature"]
    Gain = Grx, # receiver gain, dB
    NoisePower = NP, # noise power
    LossFactor=LFrx, # losses in the receiving tract
    NoiseFigure = NF, # the noise factor in the path
    SampleRate = fs # sampling rate
);

println("----------------- The transmitter-----------------------------")
println("Transmitter power: $(round(Pt*1e-3;sigdigits=5)) kW")
println("Transmitter gain: $(round(Gtx;sigdigits=5)) dB")
println("Transmission path loss: $(round(LFtx;sigdigits=5)) dB")
println("--------------------------------------------------------")
println("-----------------The receiver-----------------------------")
println("Receiver gain: $(round(Gtx;sigdigits=5)) dB")
println("Losses in the receiving path: $(round(LFtx;sigdigits=5)) dB")
println("Noise temperature: $(round(LFrx;sigdigits=5)) dB")
println("Noise factor: $(round(NF;sigdigits=5)) dB")
println("Noise temperature: $(round(TN;sigdigits=5)) °C")
println("Noise power: $(round(NP;sigdigits=5)) W")

az_scan = -20:0.5:20 # azimuth scanning range
el_scan = 0:0.5:30 # scanning range by location angle
n_sig = 4 # number of ratings (number of goals)

algorithm = "MUSIC" # @param ["Beamscan","MVDR","MUSIC"]

 beamscan = EngeePhased.BeamscanEstimator2D( 
    SensorArray = array, # The antenna
    OperatingFrequency = fc, # carrier frequency
    AzimuthScanAngles = az_scan,  # azimuth scanning range
    ElevationScanAngles = el_scan,  # scanning range by location angle
    DOAOutputPort = true, # enabling the output of bearing estimates
    NumSignals = n_sig # Number of ratings
);

 mvdr = EngeePhased.MVDREstimator2D( 
    SensorArray = array, # The antenna
    OperatingFrequency = fc, # carrier frequency
    AzimuthScanAngles = az_scan,  # azimuth scanning range
    ElevationScanAngles = el_scan,  # scanning range by location angle
    DOAOutputPort = true, # enabling the output of bearing estimates
    NumSignals = n_sig # Number of ratings
);

music =  EngeePhased.MUSICEstimator2D(
    SensorArray = array, # The antenna
    OperatingFrequency = fc, # carrier frequency
    AzimuthScanAngles = az_scan,  # azimuth scanning range
    ElevationScanAngles = el_scan,  # scanning range by location angle
    DOAOutputPort = true, # enabling the output of bearing estimates
    NumSignalsSource="Property",
    NumSignals = n_sig, # number of ratings
    ForwardBackwardAveraging=false # accounting for forward and reverse averaging
);

doa_alg =  if algorithm == "Beamscan"
    println("Algorithm: Beamscan")
    deepcopy(beamscan)
elseif algorithm == "MVDR"
    println("Algorithm: MVDR");
    deepcopy(mvdr)
elseif algorithm == "MUSIC"
    println("Algorithm: MUSIC");  
    deepcopy(music)
end;

# distribution environment parameters
is_atmosphere = true # consideration of atmospheric influence on the propagation of electromagnetic waves
T_ch = 15 # atmospheric temperature, °C
press_air = 101325.0 # atmospheric pressure, PA
vap_density = 7.5 # density of wet steam
liq_density = 0.0 # condensate density
rain_rate = 0.0 # Precipitation intensity

channel_fwd = EngeePhased.WidebandLOSChannel( 
    PropagationSpeed = c,
    SampleRate = fs,
    OperatingFrequency = fc,
    NumSubbands = NBand,
    TwoWayPropagation = false, # accounting for bidirectional signal propagation
    SpecifyAtmosphere = is_atmosphere,
    Temperature = T_ch,
    DryAirPressure = press_air,
    WaterVapourDensity = vap_density,
    LiquidWaterDensity = liq_density,
    RainRate = rain_rate
);

channel_bwd = deepcopy(channel_fwd); # distribution channel from targets to radar (identical)

data = JLD2.load("quad_rcs_pattern.jld2")

num_tgts = n_sig # number of targets (set an equal number of direction finder ratings)
rcs_az = data["az"] # Azimuth EPR accounting range (-180:180), deg
rcs_el = data["el"] # EPR accounting range by seat angle (-90:90), deg
rcs_pattern = data["rcs"] # target DOOR matrix [181x361], m^2

type_tgt = "Nonfluctuating" # target type ["Non-Fluctuating","Swerling1","Swerling2","Swerling3","Swerling4"]

# formation of broadband SDR targets
rcs_targets = EngeePhased.WidebandBackscatterRadarTarget( 
    AzimuthAngles = rcs_az, # Azimuth-based EPR accounting range
    ElevationAngles = rcs_el, # EPR accounting range by seat angle
    RCSPattern = rcs_pattern, # DOR goals
    PropagationSpeed = c, # propagation velocity, m/s
    OperatingFrequency = fc, # carrier frequency, Hz
    NumSubbands = NBand, # number of frequency sub-bands
    Model = type_tgt # The type of fluctuation of the target's EPR
);

text_title = "Visualization and DOR goals"
k_norm  = 0.6 # DOOR size scaling factor
calc_plot_rcs(rcs_pattern,rcs_az,rcs_el;title = text_title,k = k_norm)

# Goal movement model
targets_pos = [
            [1000; 100; 200];; # 1 goal
            [1200; -10; 200];; # 2 the goal
            [1100; -150; 150];; # 3 the goal
            [950; 50; 320];; # 4 the goal
] # the initial position of the target, m
targets_vel = [
            [-130; 0; 0];; # The speed of the 1st target (on the radar parallel to the x-axis)
            [-150; 0; 0];; # the speed of the 2nd target (on the radar parallel to the x-axis)
            [-200; 0; 0];; # the speed of the 3rd target (on the radar parallel to the x-axis)
            [-230; 0; 0];; # the speed of the 4th target (on the radar parallel to the x-axis)
] # target speed, m/s
type_traj = "Velocity" # The motion model is uniform ["Velocity","Acceleration","Custom"]

targets_platform = EngeePhased.Platform(
    MotionModel = type_traj,
    InitialPosition = targets_pos, # starting position
    Velocity = targets_vel # target speed
)

# Radar motion model
radar_pos = [0; 0; 0] # Radar position
radar_vel = [0; 0; 0] # Radar speed

radar_platform = EngeePhased.Platform(
    InitialPosition = radar_pos, # Radar position
    Velocity = radar_vel # Radar speed
);

abs_dist,ang_true = rangeangle(targets_pos,radar_pos)

for i in 1:num_tgts
    println("----------------Goal # $(i)----------------------------")
    println("Azimuth: $(round(ang_true[1,i];sigdigits=3))°, elevation angle: $(round(ang_true[2,i];sigdigits=3))°")
    println("Radar range to target no.$(i): $(round(abs_dist[i]*1e-3;sigdigits=5)) km")
end

viz_scene(targets_pos, radar_pos) |> display;

println("-----------------RADAR station-----------------------------")
println("Carrier frequency: $(round(fc*1e-9;sigdigits=3)) GHz")
println("Speed of light: $c m/s")
println("Wavelength: $(round(lambda;sigdigits=3)) m") 
println("Number of frequency sub-bands: $NBand m")
println("Radar position: $radar_pos [Px,Py,Pz], m")
println("Radar speed: $radar_vel [Vx,Vyc,Vz],m/s")
println("--------------------------------------------------------")
println("-----------------Generator-----------------------------")
println("Sampling rate: $(round(fs*1e-6;sigdigits=5)) MHz")
println("Pulse duration: $(round(T*1e6;sigdigits=5)) microseconds")
println("Follow-up period: $(round(1e6/prf;sigdigits=5)) mks")
println("Repetition rate: $(round(prf*1e-3;sigdigits=5)) kHz")
println("Borehole: $(round(Q;sigdigits=5))")
println("Coherent accumulation time: $(round(T_coh*1e3;sigdigits=5)) ms")
println("Number of pulses per packet: $n_pulses")
println("Number of packs: $n_groups_imps")
println("Range resolution : $(round(ΔR;sigdigits=3)) m") 
println("Speed resolution : $(round(ΔV;sigdigits=3)) m/s") 
println("Maximum one-digit range : $(round(Rmax;sigdigits=3)) m") 
println("Minimum one-digit range : $(round(Rmin;sigdigits=3)) m") 
println("Maximum single-digit speed: $(round(Vmax;sigdigits=3)) m/s") 
println("--------------------------------------------------------")
println("----------------- The transmitter-----------------------------")
println("Transmitter power: $(round(Pt*1e-3;sigdigits=5)) kW")
println("Transmitter gain: $(round(Gtx;sigdigits=5)) dB")
println("Transmission path loss: $(round(LFtx;sigdigits=5)) dB")
println("--------------------------------------------------------")
println("-----------------The receiver-----------------------------")
println("Receiver gain: $(round(Gtx;sigdigits=5)) dB")
println("Losses in the receiving path: $(round(LFtx;sigdigits=5)) dB")
println("Noise temperature: $(round(LFrx;sigdigits=5)) dB")
println("Noise factor: $(round(NF;sigdigits=5)) dB")
println("Noise temperature: $(round(TN;sigdigits=5)) °C")
println("Noise power: $(round(NP;sigdigits=5)) W")
println("--------------------------------------------------------")
println("-----------------The antenna-------------------------------")
println("Distance between elements : $(round(d;sigdigits=3)) m") 
println("Rectangular HEADLIGHT Dimensions: $(Nx)x$Ny [Nx, Ny]") 
println("Frequency range of radiation: $(freq_range*1e-9) GHz")
println("Azimuth resolution: $(round(Δaz;sigdigits=3)) °") 
println("Location angle resolution: $(round(Δel;sigdigits=3)) °") 
println("--------------------------------------------------------")
println("-----------------Direction finder-----------------------------")
println("Azimuth scan range: $(az_scan) deg.")
println("Scanning range by location angle: $(el_scan) deg.")
println("Number of ratings: $n_sig")
println("--------------------------------------------------------")

# Memory allocation for the output signal and bearing estimates
data_y = zeros(ComplexF64,length(sig_lfm),Nx*Ny,n_groups_imps)
data_ang_true = zeros(2, num_tgts, n_groups_imps); # true directions for goals
data_ang_est = zeros(2, num_tgts, n_groups_imps); # true directions for goals
data_spec = zeros(length(el_scan),length(az_scan), n_groups_imps); # spectral data from the viewing angles at the direction finder output

len_one_imps = round(Int,length(sig_lfm)/n_pulses) # number of counts of 1 pulse
sig_lfm_one_imps = sig_lfm[1:len_one_imps] # allocation of 1 pulse
sig_buffer_one_groups = similar(data_y[:,:,1]) # buffer for accumulating a burst of pulses

println("The calculation of the scenario is running...")
@inbounds for i in 1:n_groups_imps
    println("Calculation of bundle number $(i)...")

    @inbounds for k in 1:n_pulses
        println("Calculation of the pulse number $(k)...")
        # Updating the position of the radar and targets in increments of the pulse period 1/prf
        radar_pos, radar_vel = radar_platform(1/prf)
        tgts_pos, tgts_vel = targets_platform(1/prf)

        # Calculating and maintaining the target bearing
        _,ang_doa = rangeangle(tgts_pos,radar_pos)
        data_ang_true[:,:,i] .= deepcopy(ang_doa)

        # Broadband signal emission on Wednesday
        sig_tx = transmitter(sig_lfm_one_imps) # passage of the transmission path
        sig_rad = radiator(sig_tx,ang_doa) # radiating a signal into space

        # Propagation from radar to targets
        sig_fwd = channel_fwd(sig_rad,radar_pos,tgts_pos,radar_vel,tgts_vel)
        # Reflection from goals
        sig_refl = rcs_targets(sig_fwd,ang_doa)

        # Propagation from targets to radar
        sig_ret = channel_bwd(sig_refl,tgts_pos,radar_pos,tgts_vel,radar_vel)

        # the receiving tract
        y = collector(sig_ret,ang_doa) # receiving lights
        sig_buffer_one_groups[(1:len_one_imps).+(k-1)*len_one_imps,:] .= receiver(y) # pre-amplification based on receiver thermal noise (LNA)
    end

    # bearing calculation
    data_spec[:,:,i],data_ang_est[:,:,i] = doa_alg(sig_buffer_one_groups); # assessment of the range and directions of arrival of goals
    data_y[:,:,i] .= sig_buffer_one_groups; # saving the received signal
end

println("Scenario calculation is completed!")

n_imp = 1 # the number of the analyzed pulse
plot_size = (600,600) # permission
plot_spec_azel(az_scan,el_scan,data_spec[:,:,n_imp];
    title = "The power spectrum at the direction finder output (algorithm $(algorithm))",
    size = plot_size
)

sort_ang_est =  sort_ang_mse(data_ang_true,data_ang_est) # sorted angle estimation data

for k in 1:n_groups_imps
    println("--------------- Pulse pack No. $(k) -----------------------------")
    for i in 1:num_tgts
        println("   ---------------- Goal #$(i) ----------------------------")
        println("       Target: azimuth = $(round(data_ang_true[1,i,k];sigdigits=4)) °, elevation = $(round(data_ang_true[2,i,k];sigdigits=4))°");
        println("       Estimate: azimuth = $(round(sort_ang_est[1,i,k];sigdigits=4))°, seat angle = $(round(sort_ang_est[2,i,k];sigdigits=4))°");
        println("       Error: azimuth = $(round(abs(sort_ang_est[1,i,k] - data_ang_true[1,i,k]);sigdigits=4))°, seat angle = $(round(abs(sort_ang_est[2,i,k] - data_ang_true[2,i,k]);sigdigits=4))°");
    end
    println("Maximum error: $(round(maximum(abs.(data_ang_true[:,:,k].-sort_ang_est[:,:,k])); sigdigits=4) )°")
    println("-------------------------------------------------------------")
end

tgtrng = maximum(abs_dist) # target range
# calculation of the SNR with an accuracy of 3 digits
snr = round(radareqsnr(
	lambda, # wavelength
	tgtrng, # target range
	Pt, # pulse power
	T, # pulse duration
	RCS = 0.01, # EPR of the goal (choose the average value)
	Gain = Grx, # antenna gain
	Loss = LFtx+LFrx,  # losses in the path
	Ts = TN # noise temperature
)[1];sigdigits = 3)

println("Current SNR: $snr db")

y = deepcopy(data_y[:,:,end]) # copying the signal data for the last pulse
ang_trues_end =  data_ang_true[:,:,end]; # copying the true directions of goals

noise_vec = collect(20:10:70) # the vector of noise power increase from 20 to 70 dB relative to NP
snr_vec = snr .- noise_vec # vector of SNR values

println("SNR vector: $(round.(snr_vec;sigdigits=3)) dB")

# allocation of memory for an incomparable spectrum and angle estimates for MUSIC
spec_music = zeros(length(el_scan),length(az_scan), length(noise_vec)) 
ang_music = zeros(2, num_tgts, length(noise_vec)) # true directions for goals

@inbounds for i in 1:length(noise_vec)
    println("Estimation calculation for SNR $(round(snr_vec[i];sigdigits=3)) dB")
    release!(music) # resetting the system object to its initial state
    noise_sig =  sqrt(NP*10^(noise_vec[i]/10))*randn(size(y)) # new data for the noise signal
    spec_music[:,:,i],ang_music[:,:,i] = music(y.+noise_sig); # assessment of the range and directions of arrival of goals
end

fig_spec_music_vec = [] # 
@inbounds for i in 1:length(noise_vec)
    push!(fig_spec_music_vec,
        plot_spec_azel(az_scan,el_scan,spec_music[:,:,i];
            title = "MUSIC (SNR = $(round(snr_vec[i];sigdigits=3)) db)"
        )
    )
end
plot(fig_spec_music_vec...)

sort_ang_music = similar(ang_music) # sorted angle estimation data
[sort_ang_music[:,:,i] .= sort_ang_mse(ang_trues_end,ang_music[:,:,i]) for i in 1:length(snr_vec)]

# building a table
calc_table_ang_analysis(ang_trues_end,sort_ang_music,snr_vec)

# allocation of memory for an incomparable spectrum and angle estimates for MVDR
spec_mvdr = zeros(length(el_scan),length(az_scan), length(noise_vec)) 
ang_mvdr = zeros(2, num_tgts, length(noise_vec)) # true directions for goals

@inbounds for i in 1:length(noise_vec)
    println("Estimation calculation for SNR $(round(snr_vec[i];sigdigits=3)) dB")
    release!(mvdr) # resetting the system object to its initial state
    noise_sig =  sqrt(NP*10^(noise_vec[i]/10))*randn(size(y)) # new data for the noise signal
    spec_mvdr[:,:,i],ang_mvdr[:,:,i] = mvdr(y.+noise_sig); # assessment of the range and directions of arrival of goals
end

fig_spec_mvdr_vec = [] # 
@inbounds for i in 1:length(noise_vec)
    push!(fig_spec_mvdr_vec,
        plot_spec_azel(az_scan,el_scan,spec_mvdr[:,:,i];
            title = "MVDR (SNR = $(round(snr_vec[i];sigdigits=3)) db)"
        )
    )
end
plot(fig_spec_mvdr_vec...)

sort_ang_mvdr = similar(ang_mvdr) # sorted angle estimation data
[sort_ang_mvdr[:,:,i] .= sort_ang_mse(ang_trues_end,ang_mvdr[:,:,i]) for i in 1:length(snr_vec)]

# building a table
calc_table_ang_analysis(ang_trues_end,sort_ang_mvdr,snr_vec)

# allocation of memory for an incomparable spectrum and angle estimates for beamscan
spec_beamscan = zeros(length(el_scan),length(az_scan), length(noise_vec)) 
ang_beamscan = zeros(2, num_tgts, length(noise_vec)) # true directions for goals

@inbounds for i in 1:length(noise_vec)
    println("Estimation calculation for SNR $(round(snr_vec[i];sigdigits=3)) dB")
    release!(beamscan) # resetting the system object to its initial state
    noise_sig =  sqrt(NP*10^(noise_vec[i]/10))*randn(size(y)) # new data for the noise signal
    spec_beamscan[:,:,i],ang_beamscan[:,:,i] = beamscan(y.+noise_sig); # assessment of the range and directions of arrival of goals
end

fig_spec_beamscan_vec = [] # 
@inbounds for i in 1:length(noise_vec)
    push!(fig_spec_beamscan_vec,
        plot_spec_azel(az_scan,el_scan,spec_beamscan[:,:,i];
            title = "beamscan (SNR = $(round(snr_vec[i];sigdigits=3)) db)"
        )
    )
end
plot(fig_spec_beamscan_vec...)

sort_ang_beamscan = similar(ang_beamscan) # sorted angle estimation data
[sort_ang_beamscan[:,:,i] .= sort_ang_mse(ang_trues_end,ang_beamscan[:,:,i]) for i in 1:length(snr_vec)]

# building a table
calc_table_ang_analysis(ang_trues_end,sort_ang_beamscan,snr_vec)

plot_spec_2d = [];
for snr_i in eachindex(snr_vec)
    spec_beamscan_snr_0db_norm = maximum(spec_beamscan[:,:,snr_i];dims=1)[:]./maximum(spec_beamscan[:,:,snr_i])
    spec_mvdr_snr_0db_norm = maximum(spec_mvdr[:,:,snr_i];dims=1)[:]./maximum(spec_mvdr[:,:,snr_i])
    spec_music_snr_0db_norm = maximum(spec_music[:,:,snr_i];dims=1)[:]./maximum(spec_music[:,:,snr_i])

    fig = plot(az_scan,spec_beamscan_snr_0db_norm,lab="BeamScan",xlab="Azimuth, degree",
        ylab="Rated power",lw=2,
        title="Comparison of direction finding algorithms for SNR = $(round(snr_vec[snr_i];sigdigits = 3)) dB"
    )
    plot!(az_scan,spec_mvdr_snr_0db_norm,lab="MVDR",lw=2)
    plot!(az_scan,spec_music_snr_0db_norm,lab="MUSIC",lw=2)
    push!(plot_spec_2d,fig)
end
plot(plot_spec_2d...;layout = (length(plot_spec_2d),1),size = (600,300*length(plot_spec_2d)))