R/SeaSondeRCS_MUSIC.R

In SeaSondeR: Radial Metrics from SeaSonde HF-Radar Data

#### Defaults ####

#' Default Parameters for MUSIC Algorithm
#'
#' This function returns the default parameters for the MUSIC algorithm used in the SeaSondeR package.
#'
#' @details
#' The default parameters are:
#' - 40: Threshold used in \code{seasonder_MUSICCheckEigenValueRatio}.
#' - 20: Threshold used in \code{seasonder_MUSICCheckSignalPowers}.
#' - 2: Threshold used in \code{seasonder_MUSICCheckSignalMatrix}.
#' - 20: Threshold used in \code{seasonder_MUSICCheckBearingDistance}.
#'
#' @return A numeric vector containing the default parameters for the MUSIC algorithm:
#' \code{c(40, 20, 2, 20)}.
#'
#' @seealso
#' \code{\link{seasonder_MUSICTestDualSolutions}} to understand the parameters in context.
#'
#' @examples
#' # Retrieve default parameters
#' params <- seasonder_defaultMUSIC_parameters()
#' print(params)
#' @export
seasonder_defaultMUSIC_parameters <- seasonder_defaultMUSICParameters <- function(){

  c(40,20,2,20)

}


#' Default Options for the MUSIC Algorithm
#'
#' This function returns a list of default options used in the MUSIC algorithm.
#'
#' The returned list includes:
#' \itemize{
#'   \item \code{PPMIN}: Lower threshold value (default is \code{NULL}).
#'   \item \code{PWMAX}: Upper threshold value (default is \code{NULL}).
#'   \item \code{smoothNoiseLevel}: Logical flag indicating whether the noise level should be smoothed (\code{FALSE} by default).
#'   \item \code{doppler_interpolation}: Doppler interpolation factor (default is \code{2}).
#'   \item \code{MUSIC_parameters}: A numeric vector of default parameters for the MUSIC algorithm, retrieved from \code{seasonder_defaultMUSIC_parameters()}.
#'   \item \code{discard_low_SNR}: Logical flag to discard solutions with low signal-to-noise ratio (\code{TRUE} by default).
#'   \item \code{discard_no_solution}: Logical flag to discard cases with no solution (\code{TRUE} by default).
#' }
#'
#' @return A list containing the default options for the MUSIC algorithm.
#'
#' @examples
#' # Retrieve the default options for the MUSIC algorithm
#' opts <- seasonder_defaultMUSICOptions()
#' print(opts)
#' @export
#' @aliases seasonder_defaultMUSIC_options
seasonder_defaultMUSICOptions <- seasonder_defaultMUSIC_options <- function(){
  
  list(PPMIN = NULL,
       PWMAX = NULL,
       smoothNoiseLevel =FALSE,
       doppler_interpolation = 2,
       MUSIC_parameters = seasonder_defaultMUSIC_parameters(),
       discard_low_SNR = TRUE,
       discard_no_solution = TRUE
       )
  
}


#' Initialize Covariance Matrix for MUSIC Algorithm
#'
#' This function initializes a covariance matrix for use in the MUSIC algorithm.
#'
#' @details
#' The covariance matrix is initialized as a 3 x 3 matrix filled with complex \code{NA} values.
#' This structure is specifically designed for three-channel antenna configurations commonly used
#' in SeaSondeR applications.
#'
#' @return A 3 x 3 matrix of complex values, each initialized to \code{NA_complex_}.
#'
#' @seealso
#' \code{\link{seasonder_defaultMUSIC_parameters}} for default MUSIC parameters.
#'
seasonder_MUSICInitCov <- function(){
  out <- matrix(rep(NA_complex_, 9), nrow = 3)

  return(out)
}



#' Initialize Projection Matrix for MUSIC Algorithm
#'
#' This function initializes a projection matrix for use in the MUSIC algorithm.
#'
#' @param bearings A numeric vector representing the bearings (in degrees) for which projections are initialized.
#'   Defaults to \code{0}.
#'
#' @details
#' The function creates a 2 x n complex matrix, where n is the number of bearings.
#' The matrix rows are labeled:
#' - \code{"single"}: For single projections.
#' - \code{"dual"}: For dual projections.
#'
#' An attribute \code{"bearings"} is attached to the matrix, storing the input bearings vector.
#'
#' @return A 2 x n matrix of complex values, each initialized to \code{NA_complex_}, with row names
#'   \code{"single"} and \code{"dual"}. The input bearings are stored as an attribute.
#'
#' @seealso
#' \code{\link{seasonder_MUSICInitCov}} for initializing covariance matrices.
#'
seasonder_MUSICInitProjections <- function(bearings = 0) {

  # Initialize a 2 x n matrix filled with NA_complex_,
  # where n is the length of the bearings vector.
  out <- matrix(rep(NA_complex_, 2 * length(bearings)), nrow = 2)

  # Assign row names to represent projection types.
  rownames(out) <- c("single", "dual")

  # Attach the bearings vector as an attribute.
  attr(out, "bearings") <- bearings

  # Return the initialized projection matrix.
  return(out)
}


#' Initialize Direction of Arrival (DOA) Solutions for MUSIC Algorithm
#'
#' This function initializes the data structure for storing Direction of Arrival (DOA) solutions
#' calculated by the MUSIC algorithm.
#'
#' @details
#' The function returns a list containing two sub-lists, one for \code{"single"} solutions and another for \code{"dual"} solutions:
#' - \code{"single"}: Contains placeholders for single DOA solutions:
#'   - \code{bearing}: The bearing angle (\code{NA_real_} by default).
#'   - \code{a}: The complex steering vector (\code{NA_complex_} by default).
#'   - \code{P}: The power spectrum value (\code{NA_complex_} by default).
#' - \code{"dual"}: Contains placeholders for dual DOA solutions:
#'   - \code{bearing}: The bearing angle (\code{NA_real_} by default).
#'   - \code{a}: The complex steering vector (\code{NA_complex_} by default).
#'   - \code{P}: A 2 x 2 complex matrix initialized to \code{NA_complex_}.
#'
#' @return A list with initialized placeholders for \code{"single"} and \code{"dual"} DOA solutions.
#'
#' @seealso
#' \code{\link{seasonder_MUSICInitCov}} for initializing covariance matrices.
#' \code{\link{seasonder_MUSICInitProjections}} for initializing projection matrices.
#'
seasonder_MUSICInitDOASolutions <- function() {

  # Initialize a list with placeholders for single and dual DOA solutions.
  out <- list(
    single = list(
      bearing = NA_real_,  # Placeholder for bearing angle
      a = NA_complex_,     # Placeholder for complex steering vector
      P = NA_complex_,      # Placeholder for power spectrum value
      PPFG = 9,
      PWFG = 9
    ),
    dual = list(
      bearing = NA_real_,  # Placeholder for bearing angle
      a = NA_complex_,     # Placeholder for complex steering vector
      P = matrix(rep(NA_complex_, 4), nrow = 2),  # Placeholder for 2x2 complex matrix
      PPFG = c(9,9),
      PWFG = c(9,9)
    )
  )

  # Return the initialized list
  return(out)
}


#' Initialize Eigenvalue Decomposition Structure for MUSIC Algorithm
#'
#' This function initializes the data structure for storing the eigenvalue decomposition
#' results used in the MUSIC algorithm.
#'
#' @details
#' The function returns a list with the following components:
#' - \code{values}: A vector of length 3, initialized with \code{NA_complex_}, to hold the eigenvalues.
#' - \code{vectors}: A 3 x 3 matrix, initialized with \code{NA_complex_}, to hold the eigenvectors.
#'
#' This structure is designed to support three-channel antenna configurations typical in SeaSondeR applications.
#'
#' @return A list with two elements:
#' - \code{values}: Eigenvalues as a complex vector.
#' - \code{vectors}: Eigenvectors as a complex matrix.
#'
#' @seealso
#' \code{\link{seasonder_MUSICInitCov}} for initializing covariance matrices.
#'
seasonder_MUSICInitEigenDecomp <- function() {

  # Initialize a list with placeholders for eigenvalues and eigenvectors.
  out <- list(
    values = rep(NA_complex_, 3),  # Placeholder for eigenvalues
    vectors = matrix(rep(NA_complex_, 9), nrow = 3)  # Placeholder for eigenvectors
  )

  # Return the initialized list
  return(out)
}


#' Initialize Interpolated Data for MUSIC Algorithm
#'
#' This function initializes the data structure for storing interpolated cross-spectral data
#' to be used in the MUSIC algorithm.
#'
#' @param seasonder_cs_object A SeaSondeR cross-spectral object containing metadata about the number
#'   of Doppler cells and range cells.
#'
#' @details
#' The function retrieves the number of Doppler cells and range cells from the provided cross-spectral object
#' and uses this information to initialize the interpolated data structure. The resulting structure is
#' compatible with the dimensions of the cross-spectral data used in SeaSondeR.
#'
#' The data structure is initialized using \code{\link{seasonder_initCSDataStructure}}, ensuring it contains
#' placeholders for components such as \code{SSA1}, \code{SSA2}, \code{SSA3}, \code{CS12}, \code{CS13},
#' \code{CS23}, and \code{QC}.
#'
#' @return A list containing the initialized interpolated data structure with placeholders for cross-spectral components.
#'
#' @seealso
#' \code{\link{seasonder_initCSDataStructure}} for details on the cross-spectral data structure.
#'
seasonder_MUSICInitInterpolatedData <- function(seasonder_cs_object) {

  # Get the number of Doppler cells for MUSIC interpolation
  nDoppler <- seasonder_getSeaSondeRCS_MUSIC_nDopplerCells(seasonder_cs_object)

  # Get the number of range cells
  nRanges <- seasonder_getnRangeCells(seasonder_cs_object)

  # Initialize the cross-spectral data structure for interpolation
  interpolated_data <- seasonder_initCSDataStructure(nRanges = nRanges, nDoppler = nDoppler)

  # Return the initialized interpolated data structure
  return(interpolated_data)
}



#' Initialize NULL Data Structure for SeaSondeR MUSIC Analysis
#'
#' This function initializes a NULL data structure for storing results of the MUSIC analysis in
#' SeaSondeR. The structure is designed as a tibble with pre-defined columns for range cells,
#' Doppler bins, and various MUSIC-related parameters.
#'
#' @details
#' The initialized tibble contains the following columns:
#' - \code{range_cell}: Numeric vector representing range cell indices.
#' - \code{doppler_bin}: Numeric vector for Doppler bin indices.
#' - \code{range}: Numeric vector for range values.
#' - \code{freq}: Numeric vector for frequencies.
#' - \code{radial_v}: Numeric vector for radial velocities.
#' - \code{cov}: A list to store covariance matrices.
#' - \code{eigen}: A list to store eigenvalue decompositions.
#' - \code{projections}: A list to store projection matrices.
#' - \code{DOA_solutions}: A list to store Direction of Arrival (DOA) solutions.
#' - \code{eigen_values_ratio}: Numeric vector for the ratio of eigenvalues.
#' - \code{P1_check}: Logical vector indicating if the P1 criterion is satisfied.
#' - \code{retained_solution}: Character vector for the type of retained solution (\code{"single"} or \code{"dual"}).
#' - \code{DOA}: A list to store final DOA results.
#' - \code{lonlat}: A list containing a data frame with longitude (\code{lon}) and latitude (\code{lat}) values.
#'
#' @return A tibble with pre-defined columns and empty values, ready to be populated with MUSIC analysis results.
#'
#' @seealso
#' \code{\link{seasonder_MUSICInitCov}} for initializing covariance matrices.
#' \code{\link{seasonder_MUSICInitEigenDecomp}} for initializing eigenvalue decompositions.
#' \code{\link{seasonder_MUSICInitProjections}} for initializing projection matrices.
#' \code{\link{seasonder_MUSICInitDOASolutions}} for initializing DOA solutions.
#'
seasonder_NULLSeaSondeRCS_MUSIC <- function() {

  # Initialize an empty data frame with predefined columns for MUSIC analysis results
  out <- data.frame(
    range_cell = numeric(0),              # Range cell indices
    doppler_bin = numeric(0),             # Doppler bin indices
    range = numeric(0),                   # Range values
    freq = numeric(0),                    # Frequency values
    radial_v = numeric(0),                # Radial velocity values
    cov = list(),                         # List of covariance matrices
    eigen = list(),                       # List of eigenvalue decompositions
    projections = list(),                 # List of projection matrices
    DOA_solutions = list(),               # List of DOA solutions
    eigen_values_ratio = numeric(0),      # Ratio of eigenvalues
    P1_check = logical(0),                # Logical check for P1 criterion
    retained_solution = character(0),     # Type of retained solution
    DOA = list(),                         # List of DOA results
    lonlat = list(data.frame(             # Data frame for longitude and latitude
      lon = numeric(0),
      lat = numeric(0)
    ))
  )

  # Convert the data frame to a tibble
  out <- tibble::as_tibble(out)

  # Return the initialized tibble
  return(out)
}


#' Initialize SeaSondeR MUSIC Data Structure
#'
#' This function initializes a data structure for storing MUSIC analysis results
#' for a given SeaSondeR cross-spectral object.
#'
#' @param seasonder_cs_object A SeaSondeR cross-spectral object containing metadata about the radar system.
#' @param range_cells An optional vector specifying the range cells to include. Defaults to all range cells in the object.
#' @param doppler_bins An optional vector specifying the Doppler bins to include. Defaults to all Doppler bins in the object.
#'
#' @details
#' The function creates a tibble with pre-computed range, frequency, and radial velocity values
#' for the specified range cells and Doppler bins. It also initializes placeholders for MUSIC-related
#' parameters such as covariance matrices, eigen decompositions, projections, DOA solutions, and more.
#'
#' Columns in the resulting tibble include:
#' - \code{range_cell}: Range cell indices.
#' - \code{doppler_bin}: Doppler bin indices.
#' - \code{range}: Computed range values for the specified range cells.
#' - \code{freq}: Computed frequency values for the specified Doppler bins.
#' - \code{radial_v}: Computed radial velocities for the specified Doppler bins.
#' - \code{cov}: Initialized covariance matrices (see \code{\link{seasonder_MUSICInitCov}}).
#' - \code{eigen}: Initialized eigen decompositions (see \code{\link{seasonder_MUSICInitEigenDecomp}}).
#' - \code{projections}: Initialized projection matrices (see \code{\link{seasonder_MUSICInitProjections}}).
#' - \code{DOA_solutions}: Initialized DOA solutions (see \code{\link{seasonder_MUSICInitDOASolutions}}).
#' - \code{eigen_values_ratio}: Placeholder for the ratio of eigenvalues.
#' - \code{P1_check}: Logical placeholder for the P1 criterion (default is \code{TRUE}).
#' - \code{retained_solution}: Placeholder for the type of retained solution (\code{"dual"} by default).
#' - \code{DOA}: Placeholder for final DOA results.
#' - \code{lonlat}: Placeholder for longitude and latitude data as a data frame.
#'
#' @return A tibble with initialized MUSIC analysis data for the specified range cells and Doppler bins.
#'
#' @seealso
#' \code{\link{seasonder_NULLSeaSondeRCS_MUSIC}} for a NULL initialized structure.
#' \code{\link{seasonder_MUSICInitCov}}, \code{\link{seasonder_MUSICInitEigenDecomp}},
#' \code{\link{seasonder_MUSICInitProjections}}, \code{\link{seasonder_MUSICInitDOASolutions}} for initializing individual components.
#'
seasonder_initSeaSondeRCS_MUSIC <- function(seasonder_cs_object, range_cells = NULL, doppler_bins = NULL) {

  # Initialize globals for seasonder_initSeaSondeRCS_MUSIC
  range_cell <- doppler_bin <- rlang::zap()

  # Default to all range cells and Doppler bins if not specified
  if (is.null(range_cells) || is.null(doppler_bins)) {

    if (is.null(range_cells)) {
      range_cells <- 1:seasonder_getnRangeCells(seasonder_obj = seasonder_cs_object)
    }

    if (is.null(doppler_bins)) {
      doppler_bins <- 1:seasonder_getSeaSondeRCS_MUSIC_nDopplerCells(seasonder_cs_object = seasonder_cs_object)
    }

    # Create a data frame of all combinations of range cells and Doppler bins
    out <- expand.grid(range_cell = range_cells, doppler_bin = doppler_bins)

  } else {
    # Create a data frame from the specified range cells and Doppler bins
    out <- data.frame(range_cell = range_cells, doppler_bin = doppler_bins)
  }

  # Convert to a tibble
  out <- tibble::as_tibble(out)

  # Add computed columns and initialize MUSIC-related placeholders
  out <- out %>%
    dplyr::mutate(
      range = seasonder_getCellsDistKm(seasonder_cs_object)[range_cell],
      freq = seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency(seasonder_cs_object)[doppler_bin],
      radial_v = seasonder_getSeaSondeRCS_MUSIC_BinsRadialVelocity(seasonder_cs_object)[doppler_bin],
      cov = list(seasonder_MUSICInitCov()),
      MA1S = NA_real_,
      MA2S = NA_real_,
      MA3S = NA_real_,
      eigen = list(seasonder_MUSICInitEigenDecomp()),
      projections = list(seasonder_MUSICInitProjections()),
      DOA_solutions = list(seasonder_MUSICInitDOASolutions()),
      eigen_values_ratio = NA_real_,
      P1_check = TRUE,
      retained_solution = "dual",
      DOA = list(c(NA_real_, NA_real_)),
      lonlat = list(data.frame(lon = NA_real_, lat = NA_real_)),
      VFLG = 0
    )

  # Return the initialized tibble
  return(out)
}


#' Initialize MUSIC Data for SeaSondeR
#'
#' This function initializes the MUSIC data structure for a SeaSondeR cross-spectral object, including
#' optional interpolation, parameter setup, and pre-computed placeholders for MUSIC analysis.
#'
#' @param seasonder_cs_object A SeaSondeR cross-spectral object containing metadata about the radar system.
#' @param range_cells An optional vector specifying the range cells to include. Defaults to all range cells in the object.
#' @param doppler_bins An optional vector specifying the Doppler bins to include. Defaults to all Doppler bins in the object.
#' @param NULL_MUSIC Logical. If \code{TRUE}, initializes the MUSIC structure with a NULL placeholder
#'   (see \code{\link{seasonder_NULLSeaSondeRCS_MUSIC}}). Defaults to \code{FALSE}.
#'
#' @details
#' The function performs the following steps:
#' 1. Ensures the SeaSondeR object has valid interpolation and parameter settings for MUSIC analysis.
#' 2. Initializes the MUSIC data structure. If \code{NULL_MUSIC} is \code{FALSE}, the structure is
#'    populated with range cell and Doppler bin combinations.
#' 3. Computes proportion of dual solutions for MUSIC using \code{\link{seasonder_MUSICComputePropDualSols}}.
#' 4. Initializes interpolated data for cross-spectral analysis using \code{\link{seasonder_MUSICInitInterpolatedData}}.
#'
#' The final object is ready for further MUSIC analysis steps, such as computing Direction of Arrival (DOA).
#'
#' @return The updated SeaSondeR cross-spectral object with initialized MUSIC-related attributes.
#'
#' @seealso
#' \code{\link{seasonder_NULLSeaSondeRCS_MUSIC}} for initializing a NULL structure.
#' \code{\link{seasonder_initSeaSondeRCS_MUSIC}} for range and Doppler-based initialization.
#' \code{\link{seasonder_MUSICInitInterpolatedData}} for interpolated data initialization.
#'
#' @examples
#' # Minimal example for initializing MUSIC data (all range cells and Doppler bins)
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' music_obj <- seasonder_initMUSICData(cs_obj)
#'
#' # Example: specific range cells and Doppler bins
#' music_obj2 <- seasonder_initMUSICData(
#'   cs_obj,
#'   range_cells = c(1, 2),
#'   doppler_bins = c(1, 2, 5, 10)
#' )
#' @export
seasonder_initMUSICData <- function(seasonder_cs_object, range_cells = NULL, doppler_bins = NULL, NULL_MUSIC = FALSE) {

  # Copy the input SeaSondeR object
  out <- seasonder_cs_object

  # Set interpolation and parameter settings for MUSIC
  out %<>% seasonder_setSeaSondeRCS_MUSIC_doppler_interpolation(
    seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation(out)
  )
  out %<>% seasonder_setSeaSondeRCS_MUSIC_parameters(
    seasonder_getSeaSondeRCS_MUSIC_parameters(out)
  )

  out %<>% seasonder_setSeaSondeRCS_MUSIC_options(
    seasonder_getSeaSondeRCS_MUSIC_options(out)
  )

  # Initialize MUSIC data structure
  MUSIC <- seasonder_NULLSeaSondeRCS_MUSIC()

  if (!NULL_MUSIC) {
    MUSIC <- seasonder_initSeaSondeRCS_MUSIC(
      out, range_cells = range_cells, doppler_bins = doppler_bins
    )
  }
  out %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Compute properties of dual solutions for MUSIC
  out %<>% seasonder_MUSICComputePropDualSols()

  # Initialize interpolated data
  out %<>% seasonder_setSeaSondeRCS_MUSIC_interpolated_data(
    seasonder_MUSICInitInterpolatedData(out)
  )

  # Return the updated object
  return(out)
}


#### Validation ####

#' Validate Doppler Interpolation Factor for SeaSondeRCS Objects
#'
#' This function validates the \code{doppler_interpolation} factor for a \code{SeaSondeRCS} object, ensuring it is within the allowed range and does not result in exceeding the maximum number of Doppler bins after interpolation.
#'
#' @param value An integer specifying the Doppler interpolation factor. Must be one of 1, 2, 3, or 4.
#' @param seasonder_cs_object A \code{SeaSondeRCS} object containing metadata for Doppler bin calculations.
#'
#' @details
#' Doppler interpolation is a process that increases the number of Doppler bins by the specified factor before radial processing.
#' The function performs the following validations:
#' - Ensures the \code{doppler_interpolation} factor is one of 1, 2, 3, or 4.
#' - Computes the total number of Doppler bins after applying the specified interpolation factor. If this number exceeds 2048, the function aborts with a descriptive error message.
#'
#' The maximum Doppler bins (2048) constraint is derived from CODAR's SeaSonde R8 Radial Config Setup, which specifies that the product of the interpolation factor and the original number of Doppler bins should not exceed this limit.
#'
#' @return
#' The validated \code{doppler_interpolation} factor as an integer.
#'
#' @section Warnings:
#' - Using Doppler interpolation factors of 3x or 4x is not recommended.
#' - Exceeding 2048 Doppler bins after interpolation will result in an error.
#'
#' @seealso
#' \code{\link{seasonder_getnDopplerCells}} for retrieving the number of Doppler bins,
#' \code{\link{seasonder_logAndAbort}} for error handling and logging.
#'
#' @importFrom glue glue
#'
SeaSondeRCS_MUSIC_validate_doppler_interpolation <- function(value, seasonder_cs_object){
  value <- as.integer(value)

  value %in% c(1L,2L,3L,4L) || seasonder_logAndAbort(glue::glue("doppler_interpolation should be one of 1, 2, 3 or 4, but is {value}"), calling_function = "SeaSondeRCS_MUSIC_validate_doppler_interpolation")

  ndoppler <- value * seasonder_getnDopplerCells(seasonder_cs_object)

  ndoppler <= 2048 || seasonder_logAndAbort(glue::glue("The number of interpolated doppler cells should not exceed 2048, and current doppler_interpolation settion of {value} would result in {ndoppler} interpolated doppler cells. Please check."), calling_function = "SeaSondeRCS_MUSIC_validate_doppler_interpolation")

  return(value)

}

#### Processing_steps ####


SeaSondeRCS_MUSIC_start_step_text <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC algorithm started.")
}


SeaSondeRCS_MUSIC_end_step_text <- function(seasonder_cs_object) {
  # Use glue to format the message with the current system time and the provided file path

  proportion <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_MUSIC_dual_solutions_proportion()
  glue::glue("{Sys.time()}: MUSIC algorithm ended with {sprintf('%0.1f',proportion*100)}% of dual solutions.")
}

SeaSondeRCS_MUSIC_compute_cov_start_step_text <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC covariance matrix computation started")
}

SeaSondeRCS_MUSIC_covariance_decomposition_start_step_text <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC covariance matrix decomposition started")
}

SeaSondeRCS_compute_DOA_functions_start_step_text <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC DOA functions computation started")
}


SeaSondeRCS_dual_solutions_testing_start_step_text  <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC dual solutions testing started")
}

SeaSondeRCS_peak_extraction_start_step_text  <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC peak extraction started")
}


SeaSondeRCS_doa_selection_start_step_text  <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC DOA selection started")
}


SeaSondeRCS_MUSIC_compute_cov_end_step_text <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC covariance matrix computation ended")
}

SeaSondeRCS_MUSIC_covariance_decomposition_end_step_text <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC covariance matrix decomposition ended")
}

SeaSondeRCS_compute_DOA_functions_end_step_text <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC DOA functions computation ended")
}


SeaSondeRCS_dual_solutions_testing_end_step_text  <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC dual solutions testing ended")
}

SeaSondeRCS_peak_extraction_end_step_text  <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC peak extraction ended")
}


SeaSondeRCS_doa_selection_end_step_text  <- function() {
  # Use glue to format the message with the current system time and the provided file path
  glue::glue("{Sys.time()}: MUSIC DOA selection ended")
}

#### Setters ####

#' Set MUSIC Options for a SeaSondeRCS Object
#'
#' This function updates the MUSIC options stored in a SeaSondeRCS object's MUSIC data attribute.
#' It merges the provided options with the default MUSIC options, ensuring that any missing fields
#' are filled with the defaults.
#'
#' @param seasonder_cs_object A SeaSondeRCS object that contains the MUSIC data as an attribute.
#' @param MUSIC_options A named list of MUSIC options. Defaults to the output of \code{seasonder_defaultMUSICOptions()}.
#'
#' @return The updated SeaSondeRCS object with the MUSIC options stored in its MUSIC data attribute.
#'
#' @details
#' The function uses \code{modifyList} to merge the default MUSIC options with any user-specified options.
#' This ensures that the resulting options list contains all required fields.
#'
#' @examples
#' # Example: update MUSIC options on a minimal CS object
#' header <- list(nRangeCells = 1, nDopplerCells = 1)
#' data <- list(
#'   SSA1 = matrix(0,1,1), SSA2 = matrix(0,1,1), SSA3 = matrix(0,1,1),
#'   CS12 = matrix(complex(real=0,imaginary=0),1,1),
#'   CS13 = matrix(complex(real=0,imaginary=0),1,1),
#'   CS23 = matrix(complex(real=0,imaginary=0),1,1), QC = matrix(0,1,1)
#' )
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(list(header = header, data = data), 
#' seasonder_apm_object = apm_obj)
#' cs_obj <- seasonder_setMUSICOptions(cs_obj, list(doppler_interpolation = 3))
#' opts <- seasonder_getMUSICOptions(cs_obj)
#' print(opts)
#' @export
seasonder_setMUSICOptions <- seasonder_setSeaSondeRCS_MUSIC_options <- function(seasonder_cs_object, MUSIC_options = seasonder_defaultMUSICOptions()) {

  MUSIC_options <- utils::modifyList(seasonder_defaultMUSICOptions(), MUSIC_options)

  attr(seasonder_cs_object, "MUSIC_data")$MUSIC_options <- MUSIC_options


  return(seasonder_cs_object)


}


#' Set a Specific MUSIC Option for a SeaSondeRCS Object
#'
#' This function updates a single MUSIC option in the MUSIC data of a SeaSondeRCS object.
#' It verifies that the provided option name is valid (i.e. exists in the default options),
#' then updates that field with the new value.
#'
#' @param seasonder_cs_object A SeaSondeRCS object that contains the MUSIC data as an attribute.
#' @param option_name A character string specifying the name of the MUSIC option to update.
#' @param option_value The new value to assign to the specified MUSIC option.
#'
#' @return The updated SeaSondeRCS object with the specified MUSIC option modified.
#'
#' @details
#' The function first checks if \code{option_name} is one of the valid options as provided by
#' \code{seasonder_defaultMUSICOptions()}. If not, it calls \code{seasonder_logAndAbort} to log an error.
#' Then, the current MUSIC options are retrieved, updated with the new value, and stored back into the object.
#'
#' @examples
#' # Example: set a specific MUSIC option on a minimal CS object
#' header <- list(nRangeCells = 1, nDopplerCells = 1)
#' data <- list(
#'   SSA1 = matrix(0,1,1), SSA2 = matrix(0,1,1), SSA3 = matrix(0,1,1),
#'   CS12 = matrix(complex(real=0,imaginary=0),1,1),
#'   CS13 = matrix(complex(real=0,imaginary=0),1,1),
#'   CS23 = matrix(complex(real=0,imaginary=0),1,1), QC = matrix(0,1,1)
#' )
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(list(header = header, data = data), 
#' seasonder_apm_object = apm_obj)
#' cs_obj <- seasonder_setMUSICOption(cs_obj, "smoothNoiseLevel", TRUE)
#' opts <- seasonder_getMUSICOptions(cs_obj)
#' print(opts$smoothNoiseLevel)
#' @export
seasonder_setMUSICOption <- seasonder_setSeaSondeRCS_MUSIC_option <- function(seasonder_cs_object, option_name, option_value) {
  # Get the valid option names from the default MUSIC options
  valid_options <- names(seasonder_defaultMUSICOptions())
  
  # Check if the provided option_name is valid
  if (!option_name %in% valid_options) {
    seasonder_logAndAbort(
      sprintf("Invalid MUSIC option '%s'. Valid options are: %s", option_name, paste(valid_options, collapse = ", ")),
      calling_function = "seasonder_setSeaSondeRCS_MUSIC_option"
    )
  }
  
  # Retrieve current MUSIC options from the object
  current_options <- seasonder_getSeaSondeRCS_MUSIC_options(seasonder_cs_object)
  
  # Set the desired option to the new value
  current_options[[option_name]] <- option_value
  
  # Update the SeaSondeRCS object with the modified MUSIC options
  seasonder_cs_object <- seasonder_setSeaSondeRCS_MUSIC_options(seasonder_cs_object, current_options)
  
  return(seasonder_cs_object)
}


#' Set MUSIC Parameters for a SeaSondeRCS Object
#'
#' This function updates the MUSIC algorithm parameters stored in a SeaSondeRCS object's
#' MUSIC data attribute. The parameters are updated in the MUSIC options under the
#' \code{MUSIC_parameters} field. Currently, no explicit validation of the provided parameters is performed.
#'
#' @param seasonder_cs_object A \code{SeaSondeRCS} object containing cross-spectral data and metadata.
#' @param MUSIC_parameters A numeric vector of parameters for the MUSIC algorithm. Defaults to the result of \code{seasonder_defaultMUSIC_parameters()}.
#'
#' @return The updated \code{SeaSondeRCS} object with the new MUSIC parameters stored in its MUSIC options.
#'
#' @details
#' The function assigns the provided \code{MUSIC_parameters} vector to the
#' \code{MUSIC_parameters} field within the \code{MUSIC_options} list, which is stored
#' in the object's \code{MUSIC_data} attribute. These parameters are used in various
#' steps of the MUSIC processing workflow.
#'
seasonder_setSeaSondeRCS_MUSIC_parameters <- function(seasonder_cs_object, MUSIC_parameters = seasonder_defaultMUSIC_parameters()) {

  # TODO: validate MUSIC parameters

  attr(seasonder_cs_object, "MUSIC_data")$MUSIC_options$MUSIC_parameters <- MUSIC_parameters


  return(seasonder_cs_object)


}


#' Set MUSIC Data in a SeaSondeRCS Object
#'
#' This function assigns MUSIC analysis results to a SeaSondeRCS object. The MUSIC data is stored
#' within the object's \code{MUSIC_data} attribute under the field \code{MUSIC}. Currently, no explicit
#' validation is performed on the provided MUSIC data.
#'
#' @param seasonder_cs_object A \code{SeaSondeRCS} object containing cross-spectral data and metadata.
#' @param MUSIC A data structure containing the MUSIC algorithm results. This is typically a list or tibble
#'        produced during the MUSIC processing workflow.
#'
#' @return The updated \code{SeaSondeRCS} object with the specified MUSIC data stored in its attributes.
#'
#' @details
#' This low-level setter function updates the SeaSondeRCS object's MUSIC data by assigning the provided
#' MUSIC results to the \code{MUSIC} field within the \code{MUSIC_data} attribute. It is intended for use
#' during the MUSIC processing workflow.
#'
seasonder_setSeaSondeRCS_MUSIC <- function(seasonder_cs_object, MUSIC) {

  # TODO: validate MUSIC

  attr(seasonder_cs_object, "MUSIC_data")$MUSIC <- MUSIC


  return(seasonder_cs_object)


}

#' Set Dual Solutions Proportion for MUSIC Analysis
#'
#' This function assigns the dual solutions proportion to the MUSIC data of a SeaSondeRCS object.
#' The dual solutions proportion represents the fraction of solutions identified as dual in the MUSIC
#' processing workflow. Currently, no explicit validation of the provided value is performed.
#'
#' @param seasonder_cs_object A \code{SeaSondeRCS} object containing MUSIC analysis results.
#' @param dual_solutions_proportion A numeric value representing the proportion of dual solutions.
#'
#' @return The updated \code{SeaSondeRCS} object with the dual solutions proportion stored in its MUSIC data.
#'
#' @details
#' The function updates the \code{dual_solutions_proportion} field within the \code{MUSIC_data} attribute.
#' This value is later used to assess the prevalence of dual bearing solutions in the MUSIC results.
#'
seasonder_setSeaSondeRCS_MUSIC_dual_solutions_proportion <- function(seasonder_cs_object, dual_solutions_proportion) {

  # TODO: validate dual_solutions_proportion

  attr(seasonder_cs_object, "MUSIC_data")$dual_solutions_proportion <- dual_solutions_proportion


  return(seasonder_cs_object)


}


#' Set the Doppler Interpolation Factor in a SeaSondeRCS Object
#'
#' This function validates and assigns the Doppler interpolation factor in the SeaSondeRCS object, updating the corresponding option in the \code{MUSIC_data} field.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing radar data and metadata.
#' @param doppler_interpolation An integer specifying the Doppler interpolation factor. Must be 1, 2, 3, or 4.
#'
#' @details
#' The function performs the following operations:
#' \enumerate{
#'   \item Validates the value of \code{doppler_interpolation} using the function \code{SeaSondeRCS_MUSIC_validate_doppler_interpolation}.
#'   \item Updates the attribute \code{MUSIC_options$doppler_interpolation} of the SeaSondeRCS object with the validated value.
#' }
#'
#' @return The SeaSondeRCS object with the updated Doppler interpolation option.
#'
#' @seealso
#' \code{\link{SeaSondeRCS_MUSIC_validate_doppler_interpolation}} for Doppler interpolation factor validation.
#'
#' @examples
#' # Create a valid SeaSondeRCS object for examples
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_object <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' # Set the Doppler interpolation factor to 2 (internal alias)
#' cs_object <- seasonder_setSeaSondeRCS_MUSIC_doppler_interpolation(cs_object, 2)
#' @export
#' @aliases seasonder_setMUSICDopplerInterpolation
seasonder_setSeaSondeRCS_MUSIC_doppler_interpolation <- seasonder_setMUSICDopplerInterpolation <- function(seasonder_cs_object, doppler_interpolation){
  
  doppler_interpolation <- SeaSondeRCS_MUSIC_validate_doppler_interpolation(doppler_interpolation, seasonder_cs_object)
  
  attr(seasonder_cs_object, "MUSIC_data")$MUSIC_options$doppler_interpolation <- doppler_interpolation
  
  return(seasonder_cs_object)
}


#' Set Interpolated MUSIC Data in a SeaSondeRCS Object
#'
#' This function assigns the interpolated cross-spectral data to the MUSIC data attribute of a SeaSondeRCS object.
#' It stores the provided interpolated data into the \code{interpolated_data} field of the MUSIC data.
#' If no data is provided, it defaults to the output of \code{seasonder_MUSICInitInterpolatedData()}.
#'
#' @param seasonder_cs_object A \code{SeaSondeRCS} object containing cross-spectral and MUSIC data.
#' @param interpolated_data A data structure (typically a list or tibble) representing the interpolated cross-spectral data.
#'        If \code{NULL}, the function uses \code{seasonder_MUSICInitInterpolatedData()} to initialize the structure.
#'
#' @return The updated \code{SeaSondeRCS} object with the \code{interpolated_data} field set in its MUSIC data attribute.
#'
#' @details
#' The function assigns the provided interpolated data (or initializes a new data structure) to the
#' \code{interpolated_data} field within the MUSIC data attribute. This structure is intended for use in further
#' MUSIC processing steps, where interpolated cross-spectral data is required for refining the estimation of Doppler bins.
#'
seasonder_setSeaSondeRCS_MUSIC_interpolated_data <- seasonder_setMUSICInterpolatedData <- function(seasonder_cs_object, interpolated_data){

  attr(seasonder_cs_object, "MUSIC_data")$interpolated_data <- interpolated_data %||% seasonder_MUSICInitInterpolatedData(seasonder_cs_object)



  return(seasonder_cs_object)

}


seasonder_setSeaSondeRCS_MUSIC_interpolated_doppler_cells_index <- function(seasonder_cs_object, interpolated_doppler_cells_index){

  # TODO: Valiate interpolated_doppler_cells_index (should be integer in the range of 1: (nDopplerCells))


  attr(seasonder_cs_object, "MUSIC_data")$interpolated_doppler_cells_index <- interpolated_doppler_cells_index

  return(seasonder_cs_object)

}
seasonder_setMUSICInterpolatedDopplerCellsIndex <- seasonder_setSeaSondeRCS_MUSIC_interpolated_doppler_cells_index

#### Getters ####


#' Retrieve MUSIC Parameters from a SeaSondeRCS Object
#'
#' This function extracts the MUSIC algorithm parameters from a SeaSondeRCS object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC data as an attribute.
#'
#' @return A numeric vector of MUSIC parameters.
#'
#' @details
#' The function checks for the presence of MUSIC parameters in the object's MUSIC_data attribute.
#' If not found, it defaults to the values returned by \code{seasonder_defaultMUSIC_parameters()}.
#'
#' @examples
#' # Minimal example for seasonder_getSeaSondeRCS_MUSIC_parameters
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' params <- seasonder_getSeaSondeRCS_MUSIC_parameters(cs_obj)
#' print(params)
#' @export
#' @aliases seasonder_getMUSICParameters
seasonder_getSeaSondeRCS_MUSIC_parameters <- seasonder_getMUSICParameters <- function(seasonder_cs_object) {


  out <- attr(seasonder_cs_object, "MUSIC_data", exact = TRUE)$MUSIC_options$MUSIC_parameters %||% seasonder_defaultMUSIC_parameters()


  return(out)


}


#' Retrieve MUSIC Options from a SeaSondeRCS Object
#'
#' This function extracts the MUSIC options from a SeaSondeRCS object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC data as an attribute.
#'
#' @return A list of MUSIC options.
#'
#' @details
#' The function retrieves the MUSIC options from the object's MUSIC_data attribute.
#' In the absence of user-defined options, it returns the default options provided by \code{seasonder_defaultMUSICOptions()}.
#'
#' @examples
#' # Minimal example for seasonder_getMUSICOptions
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' opts <- seasonder_getMUSICOptions(cs_obj)
#' print(opts)
#' @export
seasonder_getMUSICOptions <- seasonder_getSeaSondeRCS_MUSIC_options <- function(seasonder_cs_object) {


  out <- attr(seasonder_cs_object, "MUSIC_data", exact = TRUE)$MUSIC_options %||% seasonder_defaultMUSICOptions()


  return(out)


}

#' Retrieve MUSIC Data from a SeaSondeRCS Object
#'
#' This function extracts the MUSIC data structure from a SeaSondeRCS object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC data as an attribute.
#'
#' @return The MUSIC data structure, typically a data frame or tibble with MUSIC results.
#'
#' @details
#' If the MUSIC data does not exist in the object, the function initializes it via \code{seasonder_initSeaSondeRCS_MUSIC()}.
#'
#' @examples
#' # Minimal example for seasonder_getMUSIC
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' music_data <- seasonder_getSeaSondeRCS_MUSIC(cs_obj)
#' print(music_data)
#' @export
#' @aliases seasonder_getSeaSondeRCS_MUSIC
seasonder_getSeaSondeRCS_MUSIC <- seasonder_getMUSIC <-  function(seasonder_cs_object) {


  out <- attr(seasonder_cs_object, "MUSIC_data", exact = TRUE)$MUSIC %||% seasonder_initSeaSondeRCS_MUSIC(seasonder_cs_object)


  return(out)


}


#' Retrieve Proportion of Dual Solutions from MUSIC Data
#'
#' This function extracts the proportion of dual solutions from the MUSIC data in a SeaSondeRCS object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC data.
#'
#' @return A numeric value representing the dual solutions proportion, or NA if not set.
#'
#' @details
#' The function checks the MUSIC_data attribute for a dual_solutions_proportion value.
#' If not available, it defaults to NA_real_.
#'
#' @examples
#' # Minimal example for seasonder_getMUSICDualSolutionsProportion
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' dual_prop <- seasonder_getMUSICDualSolutionsProportion(cs_obj)
#' print(dual_prop)
#' @export
seasonder_getMUSICDualSolutionsProportion <- seasonder_getSeaSondeRCS_MUSIC_dual_solutions_proportion <- function(seasonder_cs_object) {


  out <- attr(seasonder_cs_object, "MUSIC_data", exact = TRUE)$dual_solutions_proportion %||% NA_real_


  return(out)


}


#' Retrieve the Doppler Interpolation Factor from MUSIC Options
#'
#' This function obtains the Doppler interpolation factor used in the MUSIC algorithm from a SeaSondeRCS object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC data and options.
#'
#' @return An integer representing the Doppler interpolation factor.
#'
#' @details
#' The function accesses the MUSIC_data attribute under MUSIC_options and retrieves the doppler_interpolation parameter.
#' If absent, it defaults to 1L.
#'
#' @examples
#' # Assuming cs_object is a valid SeaSondeRCS object.
#' # Minimal example for seasonder_getMUSICDopplerInterpolation
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' interp_factor <- seasonder_getMUSICDopplerInterpolation(cs_obj)
#' print(interp_factor)
#' @export
#' @aliases seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation
seasonder_getMUSICDopplerInterpolation <- seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation <-   function(seasonder_cs_object){


  out <- attr(seasonder_cs_object, "MUSIC_data", exact = TRUE)$MUSIC_options$doppler_interpolation %||% 1L

  return(out)

}



#' Retrieve Interpolated MUSIC Data from a SeaSondeRCS Object
#'
#' This function extracts the interpolated MUSIC cross-spectra data from a SeaSondeRCS object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing interpolated MUSIC data as an attribute.
#'
#' @return A list representing the interpolated cross-spectra data.
#'
#' @details
#' The function first checks if the interpolated data is set in the MUSIC_data attribute.
#' If absent, it initializes the data with \code{seasonder_MUSICInitInterpolatedData()}.
#'
#' @examples
#' # Minimal example for seasonder_getMUSICInterpolatedData
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' interp_data <- seasonder_getMUSICInterpolatedData(cs_obj)
#' str(interp_data)
#' @export
seasonder_getMUSICInterpolatedData <- seasonder_getSeaSondeRCS_MUSIC_interpolated_data <- function(seasonder_cs_object){



  out <- attr(seasonder_cs_object, "MUSIC_data", exact = TRUE)$interpolated_data
  
  if(is.null(out)){
    out <- seasonder_MUSICInitInterpolatedData(seasonder_cs_object)
  }
  
  return(out)
}


#' Retrieve Interpolated Doppler Cells Index from a SeaSondeRCS Object
#'
#' This function extracts the index of interpolated Doppler cells, stored in the MUSIC_data attribute of a SeaSondeRCS object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC data.
#'
#' @return A vector of indices corresponding to the interpolated Doppler cells.
#'
#' @details
#' The interpolated doppler cells index is part of the MUSIC_data and is used to identify
#' which Doppler bins were introduced during the interpolation process.
#'
#' @examples
#' # Minimal example for seasonder_getMUSICInterpolatedDopplerCellsIndex
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' doppler_index <- seasonder_getMUSICInterpolatedDopplerCellsIndex(cs_obj)
#' print(doppler_index)
#' @export
seasonder_getMUSICInterpolatedDopplerCellsIndex <- seasonder_getSeaSondeRCS_MUSIC_interpolated_doppler_cells_index <- function(seasonder_cs_object){

  out <-  attr(seasonder_cs_object, "MUSIC_data", exact = TRUE)$interpolated_doppler_cells_index

  return(out)

}


#' Retrieve the MUSIC Configuration from a SeaSondeRCS Object
#'
#' This function returns the key configuration parameters for the MUSIC algorithm from a SeaSondeRCS object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC data and options.
#'
#' @return A list containing:
#' \itemize{
#'   \item \code{doppler_interpolation}: The Doppler interpolation factor.
#'   \item \code{MUSIC_parameters}: The numeric vector of MUSIC parameters.
#' }
#'
#' @details
#' The configuration is aggregated from the MUSIC_data attribute of the object for easy access.
#'
#' @examples
#' # Minimal example for seasonder_getMUSICConfig
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' config <- seasonder_getMUSICConfig(cs_obj)
#' print(config)
#' @export
seasonder_getMUSICConfig <- seasonder_getSeaSondeRCS_MUSICConfig <- function(seasonder_cs_object){

  out <- list(doppler_interpolation = seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation(seasonder_cs_object),
              MUSIC_parameters = seasonder_getSeaSondeRCS_MUSIC_parameters(seasonder_cs_object))


  return(out)
}

#### Derived quantities ####

#' Compute the Proportion of Dual Solutions in MUSIC Data
#'
#' This function calculates the proportion of "dual" solutions in the MUSIC data
#' associated with a given `SeaSondeRCS` object. It updates the object with the computed
#' proportion as a new attribute.
#'
#' @param seasonder_cs_object A `SeaSondeRCS` object containing MUSIC data and other related attributes.
#'
#' @return A `SeaSondeRCS` object with the calculated proportion of "dual" solutions stored as
#' an attribute. This attribute can be accessed using a relevant getter function.
#'
#' @details
#' The function performs the following steps:
#' 1. Extracts the MUSIC data from the provided `SeaSondeRCS` object.
#' 2. Computes the proportion of entries in the `retained_solution` column of the MUSIC data
#'    that are labeled as "dual".
#' 3. Updates the `SeaSondeRCS` object by adding the computed proportion as an attribute
#'    using `seasonder_setSeaSondeRCS_MUSIC_dual_solutions_proportion`.
#'
#' @seealso
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC}} to retrieve the MUSIC data.
#' \code{\link{seasonder_setSeaSondeRCS_MUSIC_dual_solutions_proportion}} to set the computed proportion.
#'
#' @importFrom magrittr %<>%
#'
seasonder_MUSICComputePropDualSols <- function(seasonder_cs_object) {

  # Retrieve the MUSIC data associated with the SeaSondeRCS object
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # Calculate the proportion of "dual" solutions in the MUSIC data
  proportion <- sum(as.integer(MUSIC$retained_solution == "dual")) / nrow(MUSIC)

  # Update the SeaSondeRCS object by adding the computed proportion of "dual" solutions
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC_dual_solutions_proportion(proportion)

  # Return the updated SeaSondeRCS object
  return(seasonder_cs_object)
}



#' Compute Power Matrix
#'
#' This function calculates the power matrix based on the provided steering vector, eigenvalues, and eigenvectors. The computation differs depending on the number of columns in the steering vector matrix.
#'
#' @param eig A list containing the eigenvalues and eigenvectors of a covariance matrix. The list should include:
#'   - \code{values}: A numeric vector of eigenvalues.
#'   - \code{vectors}: A matrix where each column is an eigenvector.
#' @param a A complex matrix representing the steering vector(s). Each column corresponds to a direction of arrival.
#'
#' @return A complex matrix representing the power matrix, calculated based on the provided eigenvalues, eigenvectors, and steering vectors. If the number of columns in \code{a} is zero, the function returns \code{NULL}.
#'
#' @details
#' The function computes the power matrix using the following steps:
#' - If \code{a} has two columns:
#'   1. Select the first two eigenvalues and their corresponding eigenvectors.
#'   2. Compute the matrix \eqn{G = a^* \cdot \text{eigVector}}, where \eqn{a^*} is the conjugate transpose of \code{a}.
#'   3. Calculate the inverse of \code{G} and its conjugate transpose.
#'   4. Compute the power matrix \eqn{P = G_{\text{inv}}^* \cdot \text{diag(eigValues)} \cdot G_{\text{inv}}}.
#' - If \code{a} has one column:
#'   1. Select the first eigenvalue and its corresponding eigenvector.
#'   2. Follow similar steps as above with single-column operations.
#'
#' If \code{a} has no columns, the function returns \code{NULL}.
#'
#' @section Mathematical Formula:
#' For a steering vector matrix \eqn{a}, eigenvectors \eqn{\text{eigVector}}, and eigenvalues \eqn{\text{eigValues}}, the power matrix is calculated as:
#' \deqn{P = G_{\text{inv}}^* \cdot \text{diag(eigValues)} \cdot G_{\text{inv}}}
#' where:
#' \eqn{G = a^* \cdot \text{eigVector}}
#' and \eqn{G_{\text{inv}}} is the inverse of \eqn{G}.
#'
#' @section References:
#' - Paolo, T. de, Cook, T., & Terrill, E. (2007). Properties of HF RADAR Compact Antenna Arrays and Their Effect on the MUSIC Algorithm. \emph{OCEANS 2007}, 1–10. doi:10.1109/oceans.2007.4449265.
#'
seasonder_computePowerMatrix <- function(eig, a) {

  # Initialize the output matrix P as NULL, to store the computed power matrix
  P <- NULL
  # Check if the matrix `a` has at least one column
  if (ncol(a) > 0) {

    # Compute the conjugate transpose of matrix `a`
    a_star <- Conj(t(a))

    # If `a` has exactly two columns, compute the power matrix for two eigenvalues
    if (ncol(a) == 2) {

      # Select the first two eigenvectors and eigenvalues
      eigVector <- eig$vectors[, c(2, 1)]
      eigValues <- eig$values[c(2, 1)]

      # Compute the matrix G as the product of the conjugate transpose of `a` and the selected eigenvectors
      G <- a_star %*% eigVector

      # Compute the inverse of G
      G_inv <- solve(G)

      # Compute the conjugate transpose of G_inv
      G_inv_t <- Conj(t(G_inv))

      # Compute the power matrix P using G_inv_t, the diagonal matrix of eigenvalues, and G_inv
      P <- G_inv_t %*% diag(eigValues) %*% G_inv

      # If `a` has exactly one column, compute the power matrix for one eigenvalue
    } else if (ncol(a) == 1) {

      # Select the first eigenvector and eigenvalue
      eigVector <- eig$vectors[, c(1), drop = FALSE]
      eigValues <- eig$values[c(1)]

      # Compute the matrix G as the product of the conjugate transpose of `a` and the selected eigenvector
      G <- a_star %*% eigVector

      # Compute the inverse of G
      G_inv <- solve(G)

      # Compute the conjugate transpose of G_inv
      G_inv_t <- Conj(t(G_inv))

      # Compute the power matrix P using G_inv_t, the eigenvalue (in a matrix form), and G_inv
      P <- G_inv_t %*% matrix(eigValues) %*% G_inv
    }
  }

  # Return the computed power matrix P
  return(P)
}


#' Compute Signal Power Matrix for MUSIC Algorithm
#'
#' This function computes the signal power matrix for each direction of arrival (DOA) solution obtained
#' from the MUSIC algorithm. It updates the MUSIC data in the provided SeaSondeRCS object with the computed
#' power matrices.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC data, including eigenvalues, eigenvectors,
#' and DOA solutions.
#'
#' @return The updated SeaSondeRCS object with the MUSIC data containing the computed power matrices for
#' both dual and single solutions.
#'
#' @details
#' The function performs the following steps:
#' 1. Retrieves the MUSIC data from the SeaSondeRCS object.
#' 2. Defines an internal function to update the DOA solutions with computed power matrices:
#'    - For dual steering vectors (\code{DOA_sol$dual$a}), computes the power matrix using
#'      \code{seasonder_computePowerMatrix} and updates \code{DOA_sol$dual$P}.
#'    - For single steering vectors (\code{DOA_sol$single$a}), computes the power matrix using
#'      \code{seasonder_computePowerMatrix} and updates \code{DOA_sol$single$P}.
#' 3. Iterates through the MUSIC data, applying the update function to each set of eigenvalues and DOA solutions.
#' 4. Updates the SeaSondeRCS object with the modified MUSIC data.
#'
#' @seealso
#' \code{\link{seasonder_computePowerMatrix}}
#'
#' @importFrom dplyr mutate
#' @importFrom purrr map2
#' @importFrom rlang zap
#'
seasonder_MUSICComputeSignalPowerMatrix <- function(seasonder_cs_object) {

  # Initialize the DOA_solutions with a placeholder indicating no prior solutions
  DOA_solutions <- rlang::zap()

  # Retrieve the MUSIC data from the SeaSondeRCS object
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # Define an internal function to update the DOA solutions with computed power matrices
  update_DOA_solutions <- function(eig, DOA_sol) {

    # Copy the current DOA solution to modify and return
    out <- DOA_sol


    out$dual$P <- NULL
    out$single$P <- NULL

    # If dual steering vectors are available, compute the dual power matrix
    if (is.matrix(DOA_sol$dual$a) && ncol(DOA_sol$dual$a) > 0) {

      # Compute the dual power matrix
      P_dual <- seasonder_computePowerMatrix(eig, DOA_sol$dual$a)

      # Update the dual power matrix in the DOA solution if computation was successful
      if (!is.null(P_dual)) {
        out$dual$P <- P_dual
      }
    }

    # Compute the single power matrix for single steering vectors only if a matrix is provided
    if (is.matrix(DOA_sol$single$a) && ncol(DOA_sol$single$a) > 0) {
      P_single <- seasonder_computePowerMatrix(eig, DOA_sol$single$a)
      # Update the single power matrix in the DOA solution if computation was successful
      if (!is.null(P_single)) {
        out$single$P <- P_single
      }
    }

    return(out) # Return the updated DOA solution
  }

  # Update DOA solutions for each eigenvalue-eigenvector pair in the MUSIC data
  MUSIC %<>% dplyr::mutate(
    DOA_solutions = purrr::map2(eigen, DOA_solutions, update_DOA_solutions)
  )

  # Update the SeaSondeRCS object with the modified MUSIC data
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Return the updated SeaSondeRCS object
  return(seasonder_cs_object)
}


seasonder_getSeaSondeRCS_MUSIC_interpolated_dataMatrix <- function(seasonder_cs_object, matrix_name) {

  if(!matrix_name %in% c("SSA1","SSA2","SSA3","CS12","CS13","CS23","QC")){
    seasonder_logAndAbort(glue::glue("Unknown data matrix name '{matrix_name}'"),calling_function = "matrix_name", class = "seasonder_unknown_data_matrix_name", seasonder_matrix_name = matrix_name)
  }
  matrix <- seasonder_getSeaSondeRCS_MUSIC_interpolated_data(seasonder_cs_object = seasonder_cs_object)[[matrix_name]]

  return(matrix)

}

#' Retrieve the Interpolated Number of Doppler Cells for MUSIC
#'
#' This function calculates the interpolated number of Doppler cells for the MUSIC
#' data in a given `SeaSondeRCS` object. It applies a Doppler interpolation factor
#' to the original number of Doppler cells.
#'
#' @param seasonder_cs_object A `SeaSondeRCS` object containing metadata and configurations
#'        related to MUSIC data processing.
#'
#' @return An integer representing the number of Doppler cells adjusted by the
#'         Doppler interpolation factor.
#'
#' @details
#' The function performs the following steps:
#' 1. Retrieves the total number of Doppler cells using `seasonder_getnDopplerCells`.
#' 2. Retrieves the Doppler interpolation factor using `seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation`.
#' 3. Multiplies the number of Doppler cells by the interpolation factor to compute
#'    the interpolated number of Doppler cells.
#'
#' @seealso
#' \code{\link{seasonder_getnDopplerCells}} to obtain the base number of Doppler cells.
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation}} to retrieve the Doppler interpolation factor.
#'
seasonder_getSeaSondeRCS_MUSIC_nDopplerCells <- function(seasonder_cs_object) {

  # Retrieve the number of Doppler cells from the SeaSondeRCS object
  out <- seasonder_getnDopplerCells(seasonder_cs_object)

  # Retrieve the Doppler interpolation factor for MUSIC data
  doppler_interpolation <- seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation(seasonder_cs_object)

  # Adjust the number of Doppler cells by multiplying with the Doppler interpolation factor
  out <- out * doppler_interpolation

  # Return the adjusted number of Doppler cells
  return(out)
}


#' Retrieve the Adjusted Doppler Spectrum Resolution for MUSIC Analysis
#'
#' This function calculates the Doppler spectrum resolution adjusted by the Doppler
#' interpolation factor for a given `SeaSondeRCS` object. The adjustment ensures
#' that the spectrum resolution reflects the impact of interpolation applied in the
#' MUSIC analysis.
#'
#' @param seasonder_cs_object A `SeaSondeRCS` object containing the data and parameters
#'        for MUSIC analysis.
#'
#' @return A numeric value representing the adjusted Doppler spectrum resolution.
#'
#' @details
#' The function performs the following steps:
#' 1. Retrieves the base Doppler spectrum resolution using \code{\link{seasonder_getDopplerSpectrumResolution}}.
#' 2. Obtains the Doppler interpolation factor using \code{\link{seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation}}.
#' 3. Divides the base resolution by the interpolation factor to compute the adjusted resolution.
#'
#' This adjustment is critical for accurately interpreting MUSIC data in cases where
#' Doppler interpolation has been applied.
#'
#' @seealso
#' \code{\link{seasonder_getDopplerSpectrumResolution}} to retrieve the base Doppler spectrum resolution.
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation}} to retrieve the Doppler interpolation factor.
#'
seasonder_getSeaSondeRCS_MUSIC_DopplerSpectrumResolution <- function(seasonder_cs_object) {

  # Retrieve the Doppler spectrum resolution from the SeaSondeRCS object
  res <- seasonder_getDopplerSpectrumResolution(seasonder_cs_object)

  # Retrieve the Doppler interpolation factor specific to the MUSIC analysis
  doppler_interpolation <- seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation(seasonder_cs_object)

  # Adjust the spectrum resolution by dividing it by the Doppler interpolation factor
  out <- res / doppler_interpolation

  # Return the adjusted Doppler spectrum resolution
  return(out)
}


#' Calculate Doppler Bins Frequencies for MUSIC Analysis
#'
#' This function computes the Doppler bin frequencies for a given `SeaSondeRCS` object,
#' incorporating adjustments from the MUSIC analysis. The computation accounts for
#' Doppler interpolation and the spectrum resolution.
#'
#' @param seasonder_cs_object A `SeaSondeRCS` object containing the data and parameters
#'        for MUSIC analysis.
#' @param normalized Logical. If `TRUE`, the returned frequencies are normalized by the
#'        positive Bragg frequency. Default is `FALSE`, returning frequencies in Hz.
#'
#' @return A numeric vector representing the frequency values for each Doppler bin.
#'         If `normalized = TRUE`, these values are dimensionless, relative to the
#'         positive Bragg frequency. Otherwise, they are in Hz.
#'
#' @details
#' The function performs the following steps:
#' 1. Retrieves the central Doppler bin corresponding to 0 frequency using
#'    \code{\link{seasonder_getSeaSondeRCS_MUSIC_CenterDopplerBin}}.
#' 2. Retrieves the total number of Doppler cells (adjusted for interpolation)
#'    using \code{\link{seasonder_getSeaSondeRCS_MUSIC_nDopplerCells}}.
#' 3. Retrieves the Doppler spectrum resolution using
#'    \code{\link{seasonder_getSeaSondeRCS_MUSIC_DopplerSpectrumResolution}}.
#' 4. Computes the Doppler bin frequencies using
#'    \code{\link{seasonder_computeDopplerBinsFrequency}}.
#'
#' The resulting Doppler bins frequencies are crucial for analyzing the spectral
#' properties of the MUSIC output.
#'
#' @seealso
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_CenterDopplerBin}} to retrieve the central bin.
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_nDopplerCells}} for the number of Doppler cells.
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_DopplerSpectrumResolution}} for the adjusted spectrum resolution.
#' \code{\link{seasonder_computeDopplerBinsFrequency}} for the frequency calculation.
#'
seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency <- function(seasonder_cs_object, normalized = FALSE) {

  # Retrieve the central Doppler bin, which corresponds to the frequency 0
  center_bin <- seasonder_getSeaSondeRCS_MUSIC_CenterDopplerBin(seasonder_cs_object)

  # Retrieve the total number of Doppler cells, including interpolation adjustments
  nDoppler <- seasonder_getSeaSondeRCS_MUSIC_nDopplerCells(seasonder_cs_object)

  # Retrieve the Doppler spectrum resolution, adjusted for interpolation
  spectra_res <- seasonder_getSeaSondeRCS_MUSIC_DopplerSpectrumResolution(seasonder_cs_object)

  # Compute the Doppler bins frequency based on the total Doppler cells, central bin, and resolution
  out <- seasonder_computeDopplerBinsFrequency(seasonder_cs_object, nDoppler, center_bin, spectra_res, normalized = normalized)

  # Return the calculated frequencies for each Doppler bin
  return(out)
}




#' Retrieve the Center Doppler Bin for MUSIC Analysis
#'
#' This function calculates the center Doppler bin for a `SeaSondeRCS` object.
#' The center bin corresponds to the Doppler bin representing zero frequency,
#' and the computation accounts for adjustments from the MUSIC Doppler interpolation.
#'
#' @param seasonder_cs_object A `SeaSondeRCS` object containing data and parameters
#'        for MUSIC analysis.
#'
#' @return An integer representing the center Doppler bin.
#'
#' @details
#' The function performs the following steps:
#' 1. Retrieves the total number of Doppler cells, including adjustments for MUSIC
#'    interpolation, using \code{\link{seasonder_getSeaSondeRCS_MUSIC_nDopplerCells}}.
#' 2. Computes the center Doppler bin using \code{\link{seasonder_computeCenterDopplerBin}},
#'    which determines the bin corresponding to zero frequency.
#'
#' The center Doppler bin is a key parameter for organizing Doppler frequency data
#' around zero and is critical for spectral analysis.
#'
#' @seealso
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_nDopplerCells}} to retrieve the adjusted number of Doppler cells.
#' \code{\link{seasonder_computeCenterDopplerBin}} for the center bin calculation.
#'
seasonder_getSeaSondeRCS_MUSIC_CenterDopplerBin <- function(seasonder_cs_object) {

  # Retrieve the total number of Doppler cells, including adjustments from MUSIC interpolation
  nDoppler <- seasonder_getSeaSondeRCS_MUSIC_nDopplerCells(seasonder_cs_object)

  # Compute the center Doppler bin based on the total number of Doppler cells
  out <- seasonder_computeCenterDopplerBin(seasonder_cs_object, nDoppler)

  # Return the computed center Doppler bin
  return(out)
}


#' Retrieve Radial Velocities for MUSIC Doppler Bins
#'
#' This function calculates the radial velocities for MUSIC Doppler bins based on the given
#' SeaSonde cross-spectral object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object representing the cross-spectral data structure.
#'                         It contains necessary metadata and Doppler frequency information.
#'
#' @return A numeric vector containing the radial velocities corresponding to each MUSIC Doppler bin.
#'
#' @details
#' The function uses the following process:
#' - It retrieves the Doppler bin frequencies using \code{\link{seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency}}.
#' - It computes the radial velocities associated with the bins using \code{\link{seasonder_computeBinsRadialVelocity}}.
#'
#' The computed velocities are returned as a numeric vector, which can be used in subsequent
#' analyses or visualizations.
#'
#' @seealso
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency}},
#' \code{\link{seasonder_computeBinsRadialVelocity}}
#'
seasonder_getSeaSondeRCS_MUSIC_BinsRadialVelocity <- function(seasonder_cs_object) {
  # Retrieve Doppler bin frequencies for the given cross-spectral object.
  freq <- seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency(seasonder_cs_object)

  # Compute the radial velocities corresponding to these Doppler bins.
  out <- seasonder_computeBinsRadialVelocity(seasonder_cs_object, freq)

  # Return the computed radial velocities.
  return(out)
}


#' Map Doppler Frequencies to Doppler Bins
#'
#' This function maps specified Doppler frequency values to the corresponding Doppler bins
#' for a given SeaSonde cross-spectral object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object representing the cross-spectral data structure.
#'                         It contains metadata and configuration for Doppler frequency and bin mapping.
#' @param doppler_values A numeric vector of Doppler frequency values to be mapped to Doppler bins.
#'
#' @return A numeric vector of Doppler bins corresponding to the input Doppler frequency values.
#'
#' @details
#' The function performs the following steps:
#' - Retrieves the unnormalized Doppler bin frequencies using \code{\link{seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency}}.
#' - Retrieves the Doppler spectrum resolution using \code{\link{seasonder_getSeaSondeRCS_MUSIC_DopplerSpectrumResolution}}.
#' - Retrieves the total number of Doppler cells using \code{\link{seasonder_getSeaSondeRCS_MUSIC_nDopplerCells}}.
#' - Computes the Doppler bin indices corresponding to the input Doppler frequency values using
#'   \code{\link{seasonder_computeDopplerFreq2Bins}}.
#'
#' This mapping is essential for translating frequency-domain values into bin indices used
#' in further data processing or visualization.
#'
#' @seealso
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency}},
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_DopplerSpectrumResolution}},
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_nDopplerCells}},
#' \code{\link{seasonder_computeDopplerFreq2Bins}}
#'
seasonder_MUSIC_DopplerFreq2Bins <- function(seasonder_cs_object, doppler_values) {
  # Retrieve the Doppler bin frequencies for the given cross-spectral object without normalization.
  doppler_freqs <- seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency(seasonder_cs_object, normalized = FALSE)

  # Get the resolution of the Doppler spectrum.
  delta_freq <- seasonder_getSeaSondeRCS_MUSIC_DopplerSpectrumResolution(seasonder_cs_object)

  # Retrieve the total number of Doppler cells.
  nDoppler <- seasonder_getSeaSondeRCS_MUSIC_nDopplerCells(seasonder_cs_object)

  # Map the given Doppler frequency values to the corresponding Doppler bins.
  out <- seasonder_computeDopplerFreq2Bins(seasonder_cs_object, doppler_values, doppler_freqs, delta_freq, nDoppler)

  # Return the mapped Doppler bins.
  return(out)
}


#' Map Doppler Bins to Doppler Frequencies
#'
#' This function retrieves the Doppler frequencies corresponding to specified Doppler bins
#' for a given SeaSonde cross-spectral object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object representing the cross-spectral data structure.
#'                         It contains metadata and configuration for Doppler frequency and bin mapping.
#' @param bins A numeric or integer vector of bin indices for which Doppler frequencies are needed.
#'
#' @return A numeric vector of Doppler frequencies corresponding to the input bin indices.
#'
#' @details
#' The function retrieves the full set of unnormalized Doppler bin frequencies using
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency}} and returns the frequencies
#' corresponding to the provided bin indices. This is useful for translating bin-domain indices
#' into physical Doppler frequency values for analysis or visualization.
#'
#' @seealso
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency}}
#'
seasonder_MUSIC_Bins2DopplerFreq <- function(seasonder_cs_object, bins) {
  # Retrieve the Doppler bin frequencies for the given cross-spectral object without normalization.
  doppler_freqs <- seasonder_getSeaSondeRCS_MUSIC_DopplerBinsFrequency(seasonder_cs_object, normalized = FALSE)

  # Return the Doppler frequencies corresponding to the specified bins.
  return(doppler_freqs[bins])
}


#### Dual solution tests ####


seasonder_MUSICCheckTwoSolutions <- function(seasonder_cs_object){
# Initialize globals for seasonder_MUSICCheckTwoSolutions
  DOA_solutions <- rlang::zap()
  # Extract MUSIC data from the object
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # Compute the ratio of signal powers for dual-bearing solutions
  MUSIC %<>% dplyr::mutate(P0_check = purrr::map_dbl(DOA_solutions,\(DOA_sol){

   out <- length(DOA_sol$dual$bearing) == 2
    return(out)
  }), .before = "P1_check")

  # Mark solutions that fail the P0 test as "single"
  MUSIC$retained_solution[!MUSIC$P0_check] <- "single"

  # Update the MUSIC data in the object
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Return the modified object
  return(seasonder_cs_object)

}

#' Validate Eigenvalue Ratio Using MUSIC Algorithm
#'
#' This function implements the P1 test for solutions derived using the MUSIC algorithm.
#' The test checks the ratio between the largest and the second-largest eigenvalues, which
#' serves as an indicator of signal quality.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing the MUSIC solutions and related data.
#'
#' @details
#' The P1 test is based on the ratio of the largest eigenvalue (lambda1) to the second-largest eigenvalue (lambda2):
#' 
#'   Ratio = lambda1 / lambda2
#'
#' This ratio is compared to a threshold defined in the MUSIC parameters to determine whether the solution
#' is considered valid. Solutions failing this test are marked as "single."
#'
#' @return The updated SeaSondeRCS object with the following modifications:
#' - A new column `eigen_values_ratio` in the MUSIC data.
#' - A logical column `P1_check` indicating whether each solution passes the P1 test.
#' - Updated `retained_solution` values for solutions that fail the test.
#'
#' @seealso
#' `seasonder_getSeaSondeRCS_MUSIC`, `seasonder_setSeaSondeRCS_MUSIC`
#'
#' @importFrom dplyr pull
#' @importFrom purrr map list_c
#'
seasonder_MUSICCheckEigenValueRatio <- function(seasonder_cs_object){
  # Extract MUSIC data from the object
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # Extract eigenvalues and keep only the first two for each solution
  MUSIC_eigen_values <- MUSIC %>% dplyr::pull("eigen") %>% purrr::map(\(eig) eig$values[1:2])

  # Calculate the ratio of the largest to the second largest eigenvalue
  eigen_values_ratio <- MUSIC_eigen_values %>% purrr::map(\(values) values[1]/values[2]) %>% purrr::list_c()

  # Add the eigenvalues ratio to the MUSIC data
  MUSIC$eigen_values_ratio <- eigen_values_ratio

  # Retrieve the threshold parameter for the P1 test
  MUSIC_parameter <- seasonder_getSeaSondeRCS_MUSIC_parameters(seasonder_cs_object) %>% magrittr::extract(1)

  # Determine whether each solution passes the P1 test
  P1_check <- eigen_values_ratio < MUSIC_parameter

  # Add the P1 check results to the MUSIC data
  MUSIC$P1_check <- P1_check

  # Mark solutions that fail the P1 check as "single"
  MUSIC$retained_solution[!P1_check] <- "single"

  # Update the MUSIC data in the object
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Return the modified object
  return(seasonder_cs_object)
}


#' Validate Signal Power Ratios Using MUSIC Algorithm
#'
#' This function implements the P2 test for solutions derived using the MUSIC algorithm.
#' The test evaluates the ratio between the largest and smallest signal powers for dual-bearing solutions.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing the MUSIC solutions and related data.
#'
#' @details
#' The P2 test is based on the ratio of the largest signal power (\eqn{P_{max}}) to the smallest signal power (\eqn{P_{min}}):
#'
#' \deqn{Ratio = \frac{P_{max}}{P_{min}}}
#'
#' This ratio is compared to a threshold defined in the MUSIC parameters. Only solutions that meet the following criteria are retained:
#' - The solution has two bearings.
#' - The signal power ratio is below the threshold.
#'
#' Solutions failing this test are marked as "single."
#'
#' @return The updated SeaSondeRCS object with the following modifications:
#' - A new column \code{signal_power_ratio} in the MUSIC data.
#' - A logical column \code{P2_check} indicating whether each solution passes the P2 test.
#' - Updated \code{retained_solution} values for solutions that fail the test.
#'
#' @seealso
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC}}, \code{\link{seasonder_setSeaSondeRCS_MUSIC}}
#'
#' @importFrom dplyr mutate
#' @importFrom purrr map_dbl
#'
seasonder_MUSICCheckSignalPowers <- function(seasonder_cs_object){
  # Initialize globals for seasonder_MUSICCheckSignalPowers
  DOA_solutions <- signal_power_ratio <- rlang::zap()
  # Extract MUSIC data from the object
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # Compute the ratio of signal powers for dual-bearing solutions
  MUSIC %<>% dplyr::mutate(signal_power_ratio = purrr::map_dbl(DOA_solutions,\(DOA_sol){
    out <- NA_real_
    if(length(DOA_sol$dual$bearing) == 2){
      P_diag <- abs(diag(DOA_sol$dual$P))
      out <- max(P_diag)/min(P_diag)
    }
    return(out)
  }), .after = "P1_check")

  # Retrieve the threshold parameter for the P2 test
  MUSIC_parameter <- seasonder_getSeaSondeRCS_MUSIC_parameters(seasonder_cs_object) %>% magrittr::extract(2)

  # Determine whether each solution passes the P2 test
  MUSIC %<>% dplyr::mutate(P2_check = purrr::map_lgl(DOA_solutions, \(DOA_sol) length(DOA_sol$dual$bearing) == 2 ) & !is.na(signal_power_ratio) & signal_power_ratio < MUSIC_parameter, .after = "signal_power_ratio")

  # Mark solutions that fail the P2 test as "single"
  MUSIC$retained_solution[!MUSIC$P2_check] <- "single"

  # Update the MUSIC data in the object
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Return the modified object
  return(seasonder_cs_object)
}


#' Validate Signal Matrix Power Ratios Using MUSIC Algorithm
#'
#' This function implements the P3 test for solutions derived using the MUSIC algorithm.
#' The test evaluates the ratio between the diagonal (P_diag) and off-diagonal (P_off-diag)
#' elements of the signal covariance matrix. Specifically, the ratio is computed as:
#'
#'   Ratio = P_off_diag / P_diag
#'
#' where P_diag is the product of the absolute values of the diagonal elements and
#' P_off_diag is the square of the absolute value of the upper-left off-diagonal element.
#'
#' The computed ratio is compared with the threshold parameter (the third element in the MUSIC parameters).
#' For each dual-bearing solution (i.e. when exactly two bearings are present), if the ratio is less than
#' the reciprocal of the threshold, the solution passes the P3 test; otherwise, it is marked as "single".
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC data (including DOA solutions and power matrices).
#'
#' @return The updated SeaSondeRCS object in which:
#' \itemize{
#'   \item A new column \code{diag_off_diag_power_ratio} is added to the MUSIC data.
#'   \item A logical column \code{P3_check} indicates if each solution passes the P3 test.
#'   \item The \code{retained_solution} field of solutions that fail the test is updated to "single".
#' }
#'
#' @details
#' For each entry in the MUSIC data, the function:
#' \enumerate{
#'   \item Extracts the covariance matrix power from the dual DOA solution (\code{DOA_sol$dual$P}).
#'   \item Computes the ratio by taking the product of the absolute diagonal elements and the square of the absolute
#'         off-diagonal element.
#'   \item Retrieves the threshold parameter for the P3 test.
#'   \item Validates each solution by checking that:
#'       \itemize{
#'         \item The solution has exactly two bearings.
#'         \item The computed ratio is available (not NA) and less than 1 divided by the threshold.
#'       }
#'   \item Updates the \code{retained_solution} field to "single" for solutions that do not pass the test.
#' }
#'
#' @seealso
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC}} to retrieve MUSIC data,
#' \code{\link{seasonder_setSeaSondeRCS_MUSIC}} to update MUSIC data,
#' and \code{\link{seasonder_getSeaSondeRCS_MUSIC_parameters}} to retrieve MUSIC parameters.
#'
seasonder_MUSICCheckSignalMatrix <- function(seasonder_cs_object){
  # Initialize globals for seasonder_MUSICCheckSignalMatrix
DOA_solutions <- diag_off_diag_power_ratio <- rlang::zap()
  # Extract MUSIC data from the object
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)
  # Compute the ratio of diagonal to off-diagonal power for dual-bearing solutions
  MUSIC %<>% dplyr::mutate(diag_off_diag_power_ratio = purrr::map_dbl(DOA_solutions,\(DOA_sol){
    out <- NA_real_
    if(length(DOA_sol$dual$bearing) == 2){
      P <- DOA_sol$dual$P
      P_diag <- abs(diag(P)) %>% prod()
      P_off_diag <- abs(P[1,2])^2
      out <- P_off_diag/P_diag
    }
    return(out)
  }), .after = "P2_check")

  # Retrieve the threshold parameter for the P3 test
  MUSIC_parameter <- seasonder_getSeaSondeRCS_MUSIC_parameters(seasonder_cs_object) %>% magrittr::extract(3)

  # Determine whether each solution passes the P3 test
  MUSIC %<>% dplyr::mutate(P3_check = purrr::map_lgl(DOA_solutions, \(DOA_sol) length(DOA_sol$dual$bearing) == 2 ) & !is.na(diag_off_diag_power_ratio) & diag_off_diag_power_ratio < 1/MUSIC_parameter, .after = "diag_off_diag_power_ratio")

  # Mark solutions that fail the P3 test as "single"
  MUSIC$retained_solution[!MUSIC$P3_check] <- "single"

  # Update the MUSIC data in the object
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Return the modified object
  return(seasonder_cs_object)
}


seasonder_MUSICCheckBearingDistance <- function(seasonder_cs_object){
  DOA_solutions <- bearing_separation <- rlang::zap()
  # Extract MUSIC data from the object
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # Compute the bearing distance for dual-bearing solutions
  MUSIC %<>% dplyr::mutate(bearing_separation = purrr::map_dbl(DOA_solutions,\(DOA_sol){
    out <- NA_real_
    if(length(DOA_sol$dual$bearing) == 2){
      dist <- abs(diff(DOA_sol$dual$bearing))
      out <- min(dist, 360 - dist)
    }
    return(out)
  }), .after = "P3_check")

  # Retrieve the threshold parameter for the P4 test
  MUSIC_parameter <- seasonder_getSeaSondeRCS_MUSIC_parameters(seasonder_cs_object) %>% magrittr::extract(4)

  # Determine whether each solution passes the P4 test
  MUSIC %<>% dplyr::mutate(P4_check = purrr::map_lgl(DOA_solutions, \(DOA_sol) length(DOA_sol$dual$bearing) == 2 ) & !is.na(bearing_separation) & bearing_separation > MUSIC_parameter, .after = "bearing_separation")

  # Mark solutions that fail the P4 test as "single"
  MUSIC$retained_solution[!MUSIC$P4_check] <- "single"

  # Update the MUSIC data in the object
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Return the modified object
  return(seasonder_cs_object)
}



#' Test Dual-Bearing Solutions Using MUSIC Algorithm
#'
#' This function applies a sequence of tests (P1, P2, and P3) to validate dual-bearing solutions
#' derived using the MUSIC algorithm. The tests evaluate the quality of solutions based on
#' eigenvalue ratios, signal power ratios, and covariance matrix power ratios.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC solutions and related data.
#'
#' @details
#' The function applies the following sequence of tests:
#' 1. **P1: Eigenvalue Ratio Test**:
#'    - Evaluates the ratio between the largest and second-largest eigenvalues.
#' 2. **P2: Signal Power Ratio Test**:
#'    - Validates the ratio of signal powers for dual-bearing solutions.
#' 3. **P3: Signal Matrix Power Ratio Test**:
#'    - Checks the ratio of diagonal to off-diagonal powers in the covariance matrix.
#'
#' Each test updates the MUSIC solutions in the input object, marking solutions that fail the tests as "single."
#' The function also logs the start and end of the testing process as part of the object's processing steps.
#'
#' @return The updated SeaSondeRCS object with validated dual-bearing solutions and recorded processing steps.
#'
#' @seealso
#' \code{\link{seasonder_MUSICCheckEigenValueRatio}}, \code{\link{seasonder_MUSICCheckSignalPowers}},
#' \code{\link{seasonder_MUSICCheckSignalMatrix}}, \code{\link{seasonder_setSeaSondeRCS_ProcessingSteps}}
#'
seasonder_MUSICTestDualSolutions <- function(seasonder_cs_object) {



  # Log the start of the dual solutions testing process
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_dual_solutions_testing_start_step_text())

  seasonder_cs_object %<>% seasonder_MUSICCheckTwoSolutions()

  # Apply the P1 test to validate eigenvalue ratios
  seasonder_cs_object %<>% seasonder_MUSICCheckEigenValueRatio()

  # Apply the P2 test to validate signal power ratios for dual-bearing solutions
  seasonder_cs_object %<>% seasonder_MUSICCheckSignalPowers()

  # Apply the P3 test to validate diagonal to off-diagonal power ratios
  seasonder_cs_object %<>% seasonder_MUSICCheckSignalMatrix()

  # Apply the P4 test to validate the angular distance between solutions
  seasonder_cs_object %<>% seasonder_MUSICCheckBearingDistance()


  # Log the end of the dual solutions testing process
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_dual_solutions_testing_end_step_text())


  # Return the updated object with validated dual solutions
  return(seasonder_cs_object)
}


##### Doppler interpolation #####


#' Perform Doppler Interpolation for SeaSonde Cross-Spectra Data
#'
#' This function performs Doppler interpolation on the cross-spectra data of a SeaSondeRCS object, preparing the data for MUSIC processing.
#' Interpolation is achieved by inserting additional Doppler bins using linear interpolation, potentially increasing the number of detected vectors while possibly smoothing the radials. The function tries to mimic CODAR's AnalyzeSpectra tool interpolation, including the addition of a wraparound Doppler cell before interpolation.
#'
#' @param seasonder_cs_object A \code{SeaSondeRCS} object containing cross-spectra data and metadata for processing.
#'
#' @details
#' Doppler interpolation increases the number of Doppler bins by a factor of 2, 3, or 4 before radial processing.
#' This is accomplished by linearly interpolating between existing bins, increasing the number of radial vectors by approximately 15% for a 2x interpolation, and yielding smoother radials. The interpolation factor is configurable via the SeaSondeRCS object's \code{doppler_interpolation} attribute and it's setter \code{seasonder_setSeaSondeRCS_MUSIC_doppler_interpolation}. The number of Doppler bins after interpolation should not exceed 2048; exceeding this limit will result in an error.
#'
#' The interpolation process is as follows:
#' 1. A wraparound Doppler cell is added to the right of the data.
#' 2. For non-quality-control (QC) matrices, linear interpolation is applied to fill in the newly added Doppler bins.
#' 3. QC matrices are updated with a default value (-1) for interpolated bins.
#'
#' @note
#' - CODAR's SeaSonde R8 Radial Config Setup documentation advises against using 3x or 4x interpolation.
#' - The function ensures the number of Doppler bins after interpolation does not exceed 2048.
#' - Doppler interpolation is a preprocessing step typically performed by CODAR's AnalyzeSpectra tool before MUSIC processing.
#'
#' @return
#' A \code{SeaSondeRCS} object with updated interpolated cross-spectra data and metadata.
#'
#'
#' @seealso
#' \code{\link{seasonder_setSeaSondeRCS_MUSIC_interpolated_data}} for setting interpolated data,
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation}} for retrieving the interpolation factor,
#' \code{\link{seasonder_setSeaSondeRCS_MUSIC_doppler_interpolation}} for setting the interpolation factor,
#' \code{\link{seasonder_initCSDataStructure}} for initializing the interpolated data structure.
#'
#' @importFrom dplyr setdiff
#' @importFrom purrr map2 reduce
#' @importFrom zoo na.approx
#' @importFrom pracma Real Imag
#'
#' @examples
#' # Doppler interpolation
#' # Create a SeaSondeRCS object for interpolation example
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' # Perform Doppler interpolation
#' out <- seasonder_SeaSondeRCSMUSICInterpolateDoppler(cs_obj)
#' @export
seasonder_SeaSondeRCSMUSICInterpolateDoppler <- function(seasonder_cs_object){

  # Initialize the output object as a copy of the input object
  out <- seasonder_cs_object

  # Extract the existing cross-spectra data from the input object
  data <- seasonder_getSeaSondeRCS_data(seasonder_cs_object)

  # Set the initial interpolated data in the output object as a copy of the original data
  out %<>% seasonder_setSeaSondeRCS_MUSIC_interpolated_data(data)

  # Retrieve the Doppler interpolation factor, indicating how much to interpolate
  doppler_interpolation <- seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation(seasonder_cs_object)

  # Proceed only if the Doppler interpolation factor is greater than 1 (indicating interpolation is needed)
  if(doppler_interpolation > 1L){

    # Get the number of Doppler and range cells after interpolation
    nDoppler <- seasonder_getSeaSondeRCS_MUSIC_nDopplerCells(seasonder_cs_object)
    nRanges <- seasonder_getnRangeCells(seasonder_cs_object)

    # Initialize an empty data structure to store the interpolated results
    interpolated_data <- seasonder_initCSDataStructure(nRanges = nRanges, nDoppler = nDoppler)

    # Map the original indices to the interpolated indices based on the interpolation factor
    index_mapping <- data.frame(
      original = 1:ncol(data[[1]]),
      mapped = (0:(ncol(data[[1]]) - 1)) * doppler_interpolation + 1
    )

    # Determine the indices of the Doppler cells that will be interpolated
    interpolated_cells <- dplyr::setdiff(1:nDoppler, index_mapping$mapped)

    # Perform the interpolation for each matrix in the cross-spectra data structure
    interpolated_data %<>% purrr::map2(names(interpolated_data), \(matrix, name){

      # Validate that the matrix name exists in the original data
      if(!name %in% names(data)){
        seasonder_logAndAbort(
          glue::glue("{name} is not a data matrix name."),
          calling_function = "seasonder_SeaSondeRCSInterpolateDoppler"
        )
      }

      # Copy the original data into the interpolated matrix based on the mapped indices
      original_matrix <- data[[name]]
      matrix[, index_mapping$mapped] <- original_matrix[, index_mapping$original]

      # Special handling for the quality control (QC) matrix
      if(name == "QC"){
        # Set interpolated cells in the QC matrix to a default value (-1)
        matrix[, interpolated_cells] <- -1L
      } else {
        # For other matrices, perform linear interpolation for the interpolated cells
        matrix <- 1:nrow(matrix) %>% purrr::reduce(\(matrix_so_far, i){

          # Prepare the row data for interpolation, adding a wraparound element
          data <- c(matrix_so_far[i,, drop = TRUE], matrix_so_far[i, 1, drop = TRUE])

          if(!rlang::is_complex(data)){
            # Perform linear interpolation for real-valued data
            data[interpolated_cells] <- zoo::na.approx(abs(data))[interpolated_cells]
            data <- data[-length(data)]  # Remove the wraparound element
          } else {
            # Perform linear interpolation for complex-valued data
            data <- complex(
              real = zoo::na.approx(pracma::Real(data))[-length(data)],
              imaginary = zoo::na.approx(pracma::Imag(data))[-length(data)]
            )
          }

          # Update the row in the matrix
          matrix_so_far[i, ] <- data

          return(matrix_so_far)
        }, .init = matrix)
      }

      return(matrix)
    })

    # Set the indices of the interpolated Doppler cells in the output object
    out %<>% seasonder_setSeaSondeRCS_MUSIC_interpolated_doppler_cells_index(interpolated_cells)

    # Set the interpolated cross-spectra data in the output object
    out %<>% seasonder_setSeaSondeRCS_MUSIC_interpolated_data(interpolated_data)
  }

  # Return the updated object with interpolated data
  return(out)
}



#### MUSIC algorithm ####


#' Calculate the MUSIC Covariance Matrix for each Given Cell Range and Doppler Bin
#'
#' This function computes the Multiple Signal Classification (MUSIC) covariance matrix
#' for each cell range and Doppler bin from SeaSonde Cross Spectra (CS) data. The MUSIC
#' algorithm is used in direction finding and spectral estimation.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing the cross-spectra data.
#'
#' @return A SeaSondeRCS object updated with a computed 3x3 complex covariance matrix for each cell range and Doppler bin.
#'         The covariance matrix is stored in the MUSIC data field. Each matrix element \eqn{C_{ij}} is calculated
#'         based on auto-spectra (for diagonal elements) or cross-spectra (for off-diagonal elements).
#'         - Diagonal elements (\eqn{i = j}) are derived from auto-spectra `SSA{i}`.
#'         - Off-diagonal elements (\eqn{i \neq j}) are derived from cross-spectra `CSij`.
#'         - Auto-spectra values for the third antenna (\code{SSA3}) are taken as absolute values to comply
#'           with CODAR's recommendation to handle negative values indicating noise or interference.
#'
#' @details
#' The MUSIC algorithm estimates the direction of arrival (DOA) of signals, requiring the computation of a
#' covariance matrix from sensor data. This function constructs the covariance matrix by iterating through
#' the auto-spectra (`SSA{i}`) and cross-spectra (`CSij`) fields of the cross-spectra data.
#'
#' For diagonal elements (\eqn{i = j}), the matrix uses data from the auto-spectra field corresponding to
#' the antenna index (\code{SSA1}, \code{SSA2}, or \code{SSA3}). Negative values in \code{SSA3}, which
#' indicate noise or interference, are converted to their absolute values before use, as per the
#' Cross Spectra File Format Version 6 guidelines.
#'
#' Off-diagonal elements (\eqn{i \neq j}) are derived from cross-spectra fields, such as \code{CS12} or \code{CS23}.
#' If the row index is greater than the column index, the conjugate of the value is used.
#'
#' @seealso
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC}}, \code{\link{seasonder_setSeaSondeRCS_MUSIC}}
#'
#' @references
#' Cross Spectra File Format Version 6, CODAR. (2016).
#' Paolo, T. de, Cook, T. & Terrill, E. Properties of HF RADAR Compact Antenna Arrays and Their Effect on the MUSIC Algorithm. OCEANS 2007 1–10 (2007) doi:10.1109/oceans.2007.4449265.
#'
seasonder_MUSICComputeCov <- function(seasonder_cs_object) {

  # Initialize globals for seasonder_MUSICComputeCov
  range_cell <- doppler_bin <- rlang::zap()
  # Log the start of the MUSIC covariance matrix computation
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_ProcessingSteps(
    SeaSondeRCS_MUSIC_compute_cov_start_step_text()
  )

  # Retrieve the MUSIC data object from the input and ensure it's a data frame
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)
  MUSIC <- as.data.frame(MUSIC)

SSA1 <- seasonder_getSeaSondeRCS_MUSIC_interpolated_dataMatrix(seasonder_cs_object, "SSA1")
SSA2 <- seasonder_getSeaSondeRCS_MUSIC_interpolated_dataMatrix(seasonder_cs_object, "SSA2")
SSA3 <- seasonder_getSeaSondeRCS_MUSIC_interpolated_dataMatrix(seasonder_cs_object, "SSA3")
CS12 <- seasonder_getSeaSondeRCS_MUSIC_interpolated_dataMatrix(seasonder_cs_object, "CS12")
CS13 <- seasonder_getSeaSondeRCS_MUSIC_interpolated_dataMatrix(seasonder_cs_object, "CS13")
CS23 <- seasonder_getSeaSondeRCS_MUSIC_interpolated_dataMatrix(seasonder_cs_object, "CS23")
  # Update the MUSIC data by computing the covariance matrix for each range cell and Doppler bin
  MUSIC %<>% dplyr::mutate(
    cov = purrr::map2(range_cell, doppler_bin, \(r, d, ss1 = SSA1, cs12 = CS12, ss2 = SSA2, cs13 = CS13, ss3 = SSA3, cs23 = CS23) {
      # Initialize a 3x3 complex matrix for the covariance computation
      out <- seasonder_MUSICInitCov()

      # Iterate over the rows and columns of the covariance matrix
      for (i in 1:3) {
        for (j in 1:3) {
          if (i == j) {
            # Diagonal elements: Retrieve the auto-spectra from the corresponding SSA matrix
            value <- get(sprintf("ss%d",i))[r, d]

            # For the third antenna, take the absolute value of the auto-spectra
            if (i == 3) {
              value <- abs(value)
            }
          } else {
            # Off-diagonal elements: Retrieve the cross-spectra from the corresponding CS matrix
            
            # Conjugate the value if the row index is greater than the column index
            if (i > j) {
              value <- get(sprintf("cs%d%d",j,i))[r, d]
              value <- Conj(value)
            }else{
              value <- get(sprintf("cs%d%d",i,j))[r, d]
            }
          }

          # Assign the computed value to the covariance matrix
          out[i, j] <- value
        }
      }

      # Return the computed covariance matrix
      return(out)
    })
  )

  # Update the MUSIC data object with the computed covariance matrices
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Log the end of the MUSIC covariance matrix computation
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_ProcessingSteps(
    SeaSondeRCS_MUSIC_compute_cov_end_step_text()
  )

  # Return the updated SeaSondeRCS object
  return(seasonder_cs_object)
}



seasonder_computeSignalSNR <- function(seasonder_cs_object){

  # Initialize globals for seasonder_computeSignalSNR
  MA1S <- MA2S <- MA3S <- rlang::zap()
  # Retrieve the MUSIC data object from the input
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)
NL1 <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_NoiseLevel(dB = TRUE, antenna = 1)
NL2 <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_NoiseLevel(dB = TRUE, antenna = 2)
NL3 <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_NoiseLevel(dB = TRUE, antenna = 3)

new_columns <-  purrr::map2(MUSIC$cov, MUSIC$range_cell, \(C, rc,n1 = NL1, n2 = NL2, n3 = NL3){
  C_diag <- diag(C)
  C_diag_db <- seasonder_SelfSpectra2dB(seasonder_cs_object, C_diag)
 MA1S <- C_diag_db[1]- n1[rc]
  MA2S <- C_diag_db[2]- n2[rc]
  MA3S <- C_diag_db[3]- n3[rc]
  list(MA1S = MA1S,
             MA2S = MA2S,
             MA3S = MA3S)
} )  %>% purrr::transpose() %>% purrr::map(unlist) %>% as.data.frame()



  MUSIC <- MUSIC %>% dplyr::select(-c(MA1S,MA2S,MA3S)) %>% dplyr::bind_cols(new_columns)
  # Update the MUSIC data object with the computed covariance matrices
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  return(seasonder_cs_object)

}

seasonder_SNRCheck <- function(seasonder_cs_object, discard_low_SNR  = TRUE){
  # Initialize globals for seasonder_SNRCheck
  MA1S <- MA2S <- MA3S <- VFLG <- rlang::zap()
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)
  noisefact <- seasonder_getFOR_noisefact(seasonder_cs_object)

  SNR_threshold <- 10*log10(noisefact)

  passes_SNR_check <- function(MA1S,MA2S,MA3S, SNR_threshold)  MA3S >= SNR_threshold & (MA1S >= SNR_threshold | MA2S >= SNR_threshold)

  if(discard_low_SNR){

    MUSIC %<>% dplyr::filter(passes_SNR_check(MA1S,MA2S,MA3S, SNR_threshold))
  }else{

    MUSIC %<>% dplyr::mutate(VFLG = VFLG + 64L * as.integer(!passes_SNR_check(MA1S,MA2S,MA3S, SNR_threshold)))

  }

  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  return(seasonder_cs_object)

}

seasonder_eigen_decomp_C <- function(C){
  # Initialize an empty list to store the eigen decomposition results
  out <- seasonder_MUSICInitEigenDecomp()

  # Perform the eigen decomposition of the symmetric complex covariance matrix
  eigen_decomp <- eigen(C, symmetric = TRUE)

  # Extract eigenvalues and store them in the output list
  out$values <- eigen_decomp$values

  vectors <- eigen_decomp$vectors

  for(i in 1:3){

    theta <- Arg(vectors[3,i])
    vectors[,i] <- vectors[,i] * exp(-1i * theta)
  }

  # Extract eigenvectors and store them in the output list
  out$vectors <- vectors

  # Return the list containing eigenvalues and eigenvectors
  return(out)
}




#' Eigen Decomposition of the MUSIC Covariance Matrix
#'
#' Performs the eigen decomposition of a MUSIC covariance matrix to obtain the eigenvalues
#' and eigenvectors. This decomposition is a critical step in the MUSIC algorithm for spectral
#' estimation and direction finding, as it enables the identification of the signal and noise subspaces.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing the covariance matrices derived from cross-spectra data.
#'
#' @return An updated SeaSondeRCS object where each Doppler cell includes the eigenvalues and eigenvectors of its covariance matrix.
#'         The eigenvalues are sorted in descending order, and the eigenvectors are aligned accordingly. The updates include:
#'         \itemize{
#'           \item \code{eigen$values}: A numeric vector containing the sorted eigenvalues for each Doppler cell.
#'           \item \code{eigen$vectors}: A 3x3 matrix of the corresponding eigenvectors for each Doppler cell, aligned with the eigenvalues.
#'         }
#'
#' @details
#' The covariance matrix represents one Doppler cell of the averaged cross-spectra of three received signals. This matrix captures
#' the summation of signals from all bearings (plus noise) received by the antennas. To estimate the direction of arrival (DOA),
#' the covariance matrix is subjected to eigenvalue decomposition (diagonalization) to estimate the signal and noise subspaces.
#'
#' In practical HF radar systems, there are two primary sources of noise:
#' 1. \emph{System (thermal) noise}: Generated by the receiving equipment and assumed to be uncorrelated between antennas.
#' 2. \emph{Spatial noise field}: Includes wind-wave noise and current noise, modeled as Gaussian, which introduces correlation.
#'
#' The eigenvalue decomposition produces:
#' - Three eigenvalues, ordered from largest to smallest.
#' - Three corresponding eigenvectors forming a 3-dimensional orthonormal basis.
#'
#' Based on the largest eigenvalues:
#' - If there is one signal present, the first eigenvector defines a 1-dimensional signal subspace, and the remaining eigenvectors
#'   represent a 2-dimensional noise subspace.
#' - If two signals are present, the first two eigenvectors form a 2-dimensional signal subspace, while the remaining eigenvector
#'   represents a 1-dimensional noise subspace.
#'
#' The signal and noise subspaces are orthogonal. This decomposition facilitates identifying the signal's direction by finding the
#' antenna manifold that best fits the signal subspace.
#'
#' @references Paolo, T. de, Cook, T. & Terrill, E. Properties of HF RADAR Compact Antenna Arrays and Their Effect on the MUSIC Algorithm. OCEANS 2007 1–10 (2007) doi:10.1109/oceans.2007.4449265.
#'
#' @seealso
#' \code{\link{seasonder_MUSICComputeCov}} for computing the covariance matrix.
#'
seasonder_MUSICCovDecomposition <- function(seasonder_cs_object){

  # Initialize globals for seasonder_MUSICCovDecomposition
  cov <- rlang::zap()
  # Add a processing step entry indicating the start of the MUSIC covariance decomposition
  seasonder_cs_object  %<>% seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_MUSIC_covariance_decomposition_start_step_text())

  # Retrieve the MUSIC data structure from the SeaSondeRCS object
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # Apply the eigen decomposition function to each covariance matrix in the MUSIC data
  MUSIC %<>% dplyr::mutate(eigen = purrr::map(cov, seasonder_eigen_decomp_C))

  # Update the MUSIC data in the SeaSondeRCS object with the new eigen decomposition results
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Add a processing step entry indicating the end of the MUSIC covariance decomposition
  seasonder_cs_object  %<>% seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_MUSIC_covariance_decomposition_end_step_text())

  # Return the updated SeaSondeRCS object
  return(seasonder_cs_object)
}




#' Compute Antenna Pattern Projections for the MUSIC Algorithm
#'
#' This function computes the projection of the antenna pattern vector onto the noise subspace,
#' a critical step in the Multiple Signal Classification (MUSIC) algorithm. It is used to estimate
#' the direction of arrival (DOA) by identifying the bearing that minimizes this projection.
#'
#' @param En A matrix containing the eigenvectors of the noise subspace, derived from the
#'           covariance matrix of the signals.
#' @param a A complex-valued vector representing the antenna manifold response for a specific
#'          bearing. Each element corresponds to the response of an antenna element.
#' @param Conj_t_a The conjugate transpose of the antenna manifold vector \code{a}.
#'
#' @return A complex scalar representing the magnitude of the projection of the antenna
#'         manifold vector onto the noise subspace. This value indicates how close the
#'         antenna manifold vector is to being orthogonal to the noise subspace.
#'
#' @details
#' The MUSIC algorithm leverages the property that the antenna manifold vector is orthogonal
#' to the noise subspace eigenvectors in an ideal scenario. However, in practice, noise in the
#' covariance matrix perturbs the noise subspace, resulting in a small but non-zero projection.
#' This function calculates the magnitude of this projection using the formula:
#'
#' \deqn{P = a^H (En E_n^H) a}
#'
#' where:
#' - \eqn{a} is the antenna manifold vector.
#' - \eqn{En} is the noise subspace eigenvector matrix.
#' - \eqn{H} denotes the Hermitian (conjugate transpose) operator.
#'
#' The bearing that produces the smallest projection is considered the best estimate of the signal bearing,
#' as it corresponds to the direction where the signal is strongest relative to the noise.
#'
#' @section References:
#' - Paolo, T. de, Cook, T., & Terrill, E. (2007). Properties of HF RADAR Compact Antenna Arrays and Their Effect on the MUSIC Algorithm. \emph{OCEANS 2007}, 1–10. doi:10.1109/oceans.2007.4449265.
#'
seasonder_compute_antenna_pattern_proyections <- function(En, a, Conj_t_a){
  # Computes the projection of the antenna pattern onto the noise subspace.
  # This function is used to calculate the Euclidean distance in the MUSIC algorithm.
  # En: A matrix containing eigenvectors corresponding to the noise subspace.
  # a: A vector representing the antenna pattern response for a specific bearing.
  # Conj_t_a: The conjugate transpose of the antenna pattern vector.
  # The formula used is a conjugate transpose of 'a' times En, times En's conjugate transpose, times 'a'.
  Conj_t_a %*% (En %*% Conj(t(En))) %*% a
}


#' Compute DOA Functions Using the MUSIC Algorithm
#'
#' This function calculates the Direction of Arrival (DOA) functions based on the MUSIC algorithm
#' for a given SeaSonde cross-spectra (CS) object. It projects the antenna patterns onto the noise
#' subspace for each Doppler bin and computes single and dual signal solutions, following the MUSIC method.
#'
#' @param seasonder_cs_object An object representing the cross-spectra (CS) data from SeaSonde.
#'
#' @return The updated `seasonder_cs_object` with the MUSIC DOA functions computed and appended.
#'
#' @details
#' The function operates as follows:
#' 1. It sets a processing step indicating the start of DOA function computation.
#' 2. Retrieves the Antenna Pattern Measurement (APM) and bearings associated with the CS object.
#' 3. Iteratively computes projections of antenna pattern responses into the noise subspace for each
#'    Doppler bin using the MUSIC algorithm. This includes:
#'    - Initializing storage for projection results.
#'    - Calculating projections for single (m = 1) and dual (m = 2) signal solutions using
#'      the eigenvectors defining the noise subspace.
#'    - For each bearing, projecting the antenna manifold vector onto the noise subspace,
#'      as described by the formula:
#'      \deqn{DOA(\theta) = \frac{1}{A^*(\theta) E_n E_n^* A(\theta)}}
#'      where:
#'      - \eqn{E_n} is the eigenvector matrix of the noise subspace.
#'      - \eqn{A(\theta)} is the antenna pattern response vector at bearing \eqn{\theta}.
#'      - \eqn{A^*(\theta)} is its conjugate transpose.
#' 4. Appends the computed DOA functions to the MUSIC data of the CS object.
#' 5. Updates the processing step to indicate completion.
#'
#' @section References:
#' - Paolo, T. de, Cook, T., & Terrill, E. (2007). Properties of HF RADAR Compact Antenna Arrays and Their Effect on the MUSIC Algorithm. \emph{OCEANS 2007}, 1–10. doi:10.1109/oceans.2007.4449265.
#'
#' @seealso
#' \code{\link{seasonder_compute_antenna_pattern_proyections}} for computing projections.
#'
#'
#' @importFrom dplyr mutate
#' @importFrom purrr map
#' @importFrom magrittr %<>%
#'
seasonder_MUSICComputeDOAProjections <- function(seasonder_cs_object){
  # Sets a processing step message indicating the start of distance computation.
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_compute_DOA_functions_start_step_text())

  # Retrieves the Antenna Pattern Measurement (APM) object associated with the cross spectra (CS) object.
  seasonder_apm_obj <- seasonder_getSeaSondeRCS_APM(seasonder_cs_object)

  # Extrapolate APM
  seasonder_apm_obj %<>% seasonder_extrapolateAPM(n = 1)

  # Retrieves the bearings associated with the APM object.
  bearings <- seasonder_getSeaSondeRAPM_BEAR(seasonder_apm_obj)

  # Retrieves the MUSIC-related data associated with the CS object.
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  a_list <- purrr::map(seq_along(bearings), \(i){
    # Extracts the antenna pattern response for the current bearing.
    a <- seasonder_apm_obj[,i, drop =FALSE]
    names(a) <- NULL
    return(a)
  })

  Conj_t_a_list <- purrr::map(a_list, \(a){
    # Computes the conjugate transpose of the antenna pattern response.
    Conj(t(a))
  })

  # Updates the MUSIC object by calculating proyections of the antenna pattern into the noise sub-space for each Doppler bin.
  MUSIC %<>% dplyr::mutate(projections = purrr::map(eigen, \(eigen_analysis){
    # Initializes an empty matrix to store the projections for single and dual solutions.
    out <- seasonder_MUSICInitProjections(bearings = bearings)

    # Iterates over the two possible signal solutions (single and dual).
    for(i in 1:2){
      # Extracts the noise subspace eigenvectors corresponding to the current solution.
      En <- eigen_analysis$vectors[,(i+1):3, drop =FALSE]

      # Iterates over all bearings to calculate the projection.
      for(j in seq_along(bearings)){

        # Calculates the projection for the current solution and bearing.
        out[i,j] <- seasonder_compute_antenna_pattern_proyections(En,a_list[[j]],Conj_t_a_list[[j]])
      }
    }

    # Returns the calculated projections for the current Doppler bin.
    return(out)
  }))

  # Updates the CS object with the modified MUSIC data.
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Sets a processing step message indicating the end of projections computation.
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_compute_DOA_functions_end_step_text())

  # Returns the updated CS object.
  return(seasonder_cs_object)
}

#' Extract Direction of Arrival (DOA) Solutions Using the MUSIC Algorithm
#'
#' This function processes a set of MUSIC projection data to extract Direction of Arrival (DOA)
#' solutions for radar signals. It implements the approach described in Paolo and Terril (2007) for HF radar
#' analysis by first reversing the projection distances to enhance peak visibility, then detecting peaks for both
#' single and dual solution cases. Finally, it maps the detected peak locations back to bearing values.
#'
#' @param projections A numeric matrix of projection data where each column represents a set of MUSIC spectra
#'   for single and dual solutions. The matrix must have an attribute named \code{"bearings"} that contains
#'   the corresponding bearing angles (in degrees) for each column.
#' @param valid_bearings A numeric vector of valid bearing values (in degrees) that are acceptable.
#'   Detected bearing peaks falling outside this set will be disregarded.
#' @param seasonder_apm_obj A matrix or similar object representing the Antenna Pattern Matrix (APM).
#'   The columns of \code{seasonder_apm_obj} correspond to bearings and are used to extract antenna response
#'   information for the detected peaks.
#'
#' @return A list with two components corresponding to single and dual DOA solutions. Each component is a list
#'   containing:
#'   \itemize{
#'     \item \code{bearing}: The detected bearing(s) for the solution (in degrees).
#'     \item \code{a}: A subset of the APM data (columns) corresponding to the detected peak.
#'     \item \code{peak_resp}: The peak response value(s) at the detected peak(s), expressed in dB.
#'     \item \code{peak_width}: The width of the peak(s) calculated from the 3 dB limit, in degrees.
#'   }
#'
#' @details
#' The function proceeds as follows:
#' \enumerate{
#'   \item It retrieves the bearing angles from the attribute \code{"bearings"} of the \code{projections} matrix.
#'   \item It computes the inverse of the absolute projection values for both 'single' and 'dual' solution modes to
#'         enhance peaks.
#'   \item For single solutions, it detects the highest peak using \code{\link[pracma]{findpeaks}},
#'         then checks if the corresponding bearing is within the set of valid bearings.
#'   \item If a valid single peak is found, it calculates the response in dB and determines the peak width by finding
#'         the indices where the response exceeds the (peak response - 3 dB) threshold.
#'   \item For dual solutions, it similarly detects up to two peaks, filters them by valid bearings, and computes the
#'         response and peak width for each.
#'   \item Finally, the function populates and returns a DOA solutions structure containing both single and dual solution fields.
#' }
#'
#' @references
#' Paolo, S., & Terril, E. (2007). Detection and characterization of signals in HF radar cross-spectra using the MUSIC algorithm.
#' \emph{Journal of Atmospheric and Oceanic Technology}.
#'
#' @seealso
#' \code{\link{seasonder_MUSICExtractPeaks}}, \code{\link[pracma]{findpeaks}}
#'
seasonder_MUSICExtractDOASolutions <- function(projections, valid_bearings, seasonder_apm_obj){

  # Initialize an empty DOA solutions structure
  out <- seasonder_MUSICInitDOASolutions()

  # Extract the bearings attribute for projections
  bearings <- attr(projections,"bearings",exact = TRUE)

  # Reverse the single and dual solution projection matrices
  rev_single_solution_dist = 1/abs(projections['single',,drop = TRUE])
  rev_dual_solution_dist = 1/abs(projections['dual',,drop = TRUE])

  # Detect peaks in the single solution projections
  single_peaks_results <- pracma::findpeaks(rev_single_solution_dist, npeaks = 1, sortstr = TRUE)
  single_peak <- single_peaks_results[,2,drop = TRUE] # Peak location


  single_bearing <- bearings[single_peak]
  is_single_bearing_valid <- !is.null(single_peaks_results) && single_bearing %in% valid_bearings
  single_peak_resp <- NA
  single_peak_width <- NA
  if(is_single_bearing_valid){
    single_peak <- match(single_bearing, valid_bearings)


    # If valid peaks are found, calculate the corresponding responses


    single_peak_resp <- 10*log10(single_peaks_results[,1,drop = TRUE])
    limits_3db <- which(10*log10(rev_single_solution_dist) > (single_peak_resp - 3))

    limits_3db <- purrr::keep(split(limits_3db,cumsum(c(TRUE, diff(limits_3db) != 1))), function(x){
      single_peaks_results[,2,drop = TRUE] %in%x
    }) %>% unlist() %>% range()

    limits_3db[1] <- max(limits_3db[1]-1,1)
    limits_3db[2] <- min(limits_3db[2]+1,length(bearings))
      single_peak_width <- abs(diff(bearings[limits_3db]))

    # Populate the single DOA solution fields
    out$single$bearing <-  single_bearing



  }

  out$single$a <- seasonder_apm_obj[,single_peak, drop = FALSE]
  out$single$peak_resp <- single_peak_resp
  out$single$peak_width <- single_peak_width
  # Detect peaks in the dual solution projections
  dual_peaks_results <- pracma::findpeaks(rev_dual_solution_dist, npeaks = 2, sortstr = TRUE)

  # Populate the dual DOA solution fields
  dual_peaks_resp <- NA
  dual_peaks_width <- c(NA,NA)
  dual_peaks <- dual_peaks_results[,2,drop = TRUE]
  if (!is.null(dual_peaks_results)){



    dual_bearings <- bearings[dual_peaks]
    are_dual_bearings_valid <- dual_bearings %in% valid_bearings


    dual_peaks <- dual_peaks[are_dual_bearings_valid]
    dual_bearings <- bearings[dual_peaks]
    dual_peaks <- match(dual_bearings, valid_bearings)

    if (!is.null(dual_peaks_results) && any(are_dual_bearings_valid)) {
      dual_peaks_resp <- 10*log10(dual_peaks_results[,1,drop = TRUE])[are_dual_bearings_valid]

      for(i in 1:nrow(dual_peaks_results)){
      dual_row <- dual_peaks_results[i,, drop = FALSE]
      limits_3db <- which((10*log10(rev_dual_solution_dist)) > (dual_peaks_resp[i] - 3))

      limits_3db <- purrr::keep(split(limits_3db,cumsum(c(TRUE, diff(limits_3db) != 1))), function(x){
         dual_row[,2] %in%x
       }) %>% unlist() %>% range()

      limits_3db[1] <- max(limits_3db[1]-1,1)
      limits_3db[2] <- min(limits_3db[2]+1,length(bearings))
        dual_peaks_width[i] <- abs(diff(bearings[limits_3db]))
    }

      out$dual$bearing <- dual_bearings


    }
  }
  out$dual$a <- seasonder_apm_obj[,dual_peaks, drop = FALSE]
  out$dual$peak_resp <- dual_peaks_resp
  out$dual$peak_width <- dual_peaks_width

  return(out)
}

#' Validate Retained Solution in MUSIC Algorithm Peak Extraction
#'
#' This function verifies and adjusts the retained solution type ("single" or "dual") based on the
#' Direction of Arrival (DOA) solutions extracted using the MUSIC algorithm.
#'
#' @param ret_sol A character string specifying the initial solution type to retain. Valid values are \code{"single"} or \code{"dual"}.
#' @param DOA_sol A list containing extracted DOA solutions, as returned by \code{\link{seasonder_MUSICExtractDOASolutions}}.
#'
#' @return A character string indicating the validated solution type:
#' \itemize{
#'   \item \code{"single"}: If only one single solution bearing is valid.
#'   \item \code{"dual"}: If valid dual solution bearings are detected.
#'   \item \code{"none"}: If no valid bearings are found.
#' }
#'
#' @details
#' The function performs the following checks:
#' \enumerate{
#'   \item If the retained solution is "dual" but no valid dual solution bearings exist, it defaults to "single" if valid.
#'   \item If the retained solution is "single" but no valid single solution bearings exist, it defaults to "none".
#' }
#'
#' This validation ensures the output solutions are consistent with the detected peaks, addressing potential discrepancies
#' in the initial assumptions about the solution type.
#'
#' @seealso \code{\link{seasonder_MUSICExtractPeaks}}, \code{\link{seasonder_MUSICExtractDOASolutions}}
#'
seasonder_MUSICExtractPeaksCheckRetainedSolution <- function(ret_sol, DOA_sol){

  out <- ret_sol

  # Validate dual solutions
  if (ret_sol == "dual") {

    if (length(DOA_sol$dual$bearing) < 2) {
      if (length(DOA_sol$single$bearing ) == 1) {
        out <- "single"
      }else{
        out <- "none"
      }
    }

    # Validate single solutions
  }else if (ret_sol == "single") {
    if (length(DOA_sol$single$bearing ) != 1) {
      out <- "none"
    }
  }

  return(out)
}

#' Extract and Validate DOA Peaks Using MUSIC Algorithm
#'
#' This function processes a \code{SeaSondeRCS} object to extract Direction of Arrival (DOA) solutions using the MUSIC algorithm
#' and validates the retained solutions based on the extracted peaks.
#'
#' @param seasonder_cs_object An object of class \code{SeaSondeRCS} containing cross-spectra data processed with the MUSIC algorithm.
#'
#' @return An updated \code{SeaSondeRCS} object with the following fields modified:
#' \itemize{
#'   \item \code{MUSIC}: Contains the extracted DOA solutions.
#'   \item \code{ProcessingSteps}: Includes a log of the peak extraction process.
#' }
#'
#' @details
#' The function performs the following operations:
#' \enumerate{
#'   \item Initializes the peak extraction process and logs the start.
#'   \item Extracts DOA solutions for each set of projections using \code{\link{seasonder_MUSICExtractDOASolutions}}.
#'   \item Validates and adjusts the retained solution types using \code{\link{seasonder_MUSICExtractPeaksCheckRetainedSolution}}.
#'   \item Updates the \code{SeaSondeRCS} object with the extracted and validated solutions.
#'   \item Logs the completion of the peak extraction process.
#' }
#'
#' The MUSIC algorithm's implementation follows the theoretical framework outlined by Paolo and Terril (2007), emphasizing the identification of signal directions
#' in HF radar cross-spectra.
#'
#' @references
#' Paolo, S., & Terril, E. (2007). Detection and characterization of signals in HF radar cross-spectra using the MUSIC algorithm.
#' \emph{Journal of Atmospheric and Oceanic Technology}.
#'
#' @seealso \code{\link{seasonder_MUSICExtractDOASolutions}}, \code{\link{seasonder_MUSICExtractPeaksCheckRetainedSolution}}
#'
seasonder_MUSICExtractPeaks <- function(seasonder_cs_object){

  projections <- retained_solution <- DOA_solutions  <- rlang::zap()

  # Add a log entry to indicate the start of the peak extraction process
  seasonder_cs_object  %<>% seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_peak_extraction_start_step_text())

  # Retrieve the Antenna Pattern Matrix (APM) from the object
  seasonder_apm_obj <- seasonder_getSeaSondeRCS_APM(seasonder_cs_object)

  valid_bearings <- seasonder_getSeaSondeRAPM_BEAR(seasonder_apm_obj)

  # Retrieve the MUSIC data structure
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # Iterate over the projections matrix in the MUSIC object to extract DOA solutions
  MUSIC %<>% dplyr::mutate(DOA_solutions = purrr::map(projections, \(projection) seasonder_MUSICExtractDOASolutions(projection, valid_bearings, seasonder_apm_obj)))

  # Update the retained solution field based on DOA solutions
  MUSIC %<>% dplyr::mutate(retained_solution = purrr::map2_chr(retained_solution, DOA_solutions, seasonder_MUSICExtractPeaksCheckRetainedSolution))

  # Update the MUSIC field in the `SeaSondeRCS` object
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Log the end of the peak extraction process
  seasonder_cs_object  %<>% seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_peak_extraction_end_step_text())

  # Return the updated object
  return(seasonder_cs_object)
}



seasonder_MUSIC_remove_no_solutions <- function(seasonder_cs_object){

  # Initialize globals for seasonder_MUSIC_remove_no_solutions
  retained_solution <- rlang::zap()

  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)
  # Remove entries with no retained solution
  MUSIC %>% dplyr::filter(retained_solution != "none")

  seasonder_cs_object %<>%
    seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  return(seasonder_cs_object)

}


#' Select Direction of Arrival (DOA) from MUSIC Algorithm Results
#'
#' This function processes the results of the MUSIC algorithm, selects the relevant Direction of Arrival (DOA)
#' based on the specified retained solution, and updates the corresponding `SeaSondeRCS` object with the selected
#' DOA and updated processing steps.
#'
#' @param seasonder_cs_object A `SeaSondeRCS` object containing the results of the MUSIC algorithm and associated metadata.
#'
#' @return An updated `SeaSondeRCS` object with the selected DOA stored in the MUSIC results and updated processing steps.
#'
#' @details
#' The function performs the following steps:
#' 1. Updates the processing steps to indicate the start of the DOA selection process.
#' 2. Retrieves the MUSIC algorithm results from the `SeaSondeRCS` object.
#' 3. Maps the retained solution index to the corresponding DOA solution for each entry in the MUSIC results.
#' 4. Stores the updated MUSIC results, including the selected DOA, back into the `SeaSondeRCS` object.
#' 5. Updates the processing steps to indicate the end of the DOA selection process.
#'
#' @section Processing Steps:
#' The function appends the following processing steps to the `ProcessingSteps` attribute of the `SeaSondeRCS` object:
#' - Start of DOA selection.
#' - End of DOA selection.
#'
#' @seealso
#' \code{\link{seasonder_setSeaSondeRCS_ProcessingSteps}} to manage processing steps.
#' \code{\link{seasonder_getSeaSondeRCS_MUSIC}} to retrieve MUSIC results.
#' \code{\link{seasonder_setSeaSondeRCS_MUSIC}} to update MUSIC results.
#'
#' @importFrom purrr map2
#' @importFrom dplyr mutate
#' @importFrom magrittr %<>%
#'
seasonder_MUSICSelectDOA <- function(seasonder_cs_object) {

  # Placeholder variables to initialize empty values for DOA solutions and retained solutions
  DOA_solutions <- retained_solution <- rlang::zap()

  # Append a processing step to indicate the start of the DOA selection process
  # This logs and tracks the beginning of this operation in the processing history
  seasonder_cs_object %<>%
    seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_doa_selection_start_step_text())

  # Retrieve the MUSIC algorithm results from the SeaSondeRCS object
  # These results contain potential DOA solutions and metadata
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # Define a helper function to extract the relevant DOA solution
  extract_DOA_sol <- function(DOA_sol, retained_sol) {
    out <- NULL
    # If a retained solution is specified, select the corresponding DOA solution
    if (retained_sol != "none") {
      out <- DOA_sol[[retained_sol]]
    }
    return(out)
  }

  # For each row in the MUSIC results, map the retained solution to the corresponding DOA solution
  # This updates the MUSIC results to include only the selected DOA
  MUSIC %<>%
    dplyr::mutate(DOA = purrr::map2(DOA_solutions, retained_solution, extract_DOA_sol))

  # Update the MUSIC results in the SeaSondeRCS object with the selected DOA
  seasonder_cs_object %<>%
    seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Append a processing step to indicate the end of the DOA selection process
  # This logs and tracks the completion of this operation in the processing history
  seasonder_cs_object %<>%
    seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_doa_selection_end_step_text())

  # Return the updated SeaSondeRCS object
  # This object now includes the selected DOA in the MUSIC results
  return(seasonder_cs_object)
}


seasonder_checkPPMIN <- function(seasonder_cs_object) {

  # Initialize globals for seasonder_checkPPMIN
  DOA <- retained_solution <- rlang::zap()
  PPMIN <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_MUSIC_options() %>% purrr::pluck("PPMIN",.default = NULL)
 if(!is.null(PPMIN)){
   MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)
  check_PPMIN <- function(DOA_sol, retained_sol, PPMIN) {
    DOA_sol$PPFG <- as.integer(DOA_sol$peak_resp < PPMIN) * 3 + 1
  return(DOA_sol)
}


  MUSIC %<>%
    dplyr::mutate(DOA = purrr::map2(DOA, retained_solution, \(x, y) check_PPMIN(x, y, PPMIN)))

  seasonder_cs_object %<>%
    seasonder_setSeaSondeRCS_MUSIC(MUSIC)

 }
  return(seasonder_cs_object)
}

seasonder_checkPWMAX <- function(seasonder_cs_object) {
  # Initialize globals for seasonder_checkPWMAX
  DOA <- retained_solution <- rlang::zap()  

  PWMAX <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_MUSIC_options() %>% purrr::pluck("PWMAX",.default = NULL)
  if(!is.null(PWMAX)){
    MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)
    check_PWMAX <- function(DOA_sol, retained_sol, PWMAX) {
      DOA_sol$PWFG <- as.integer(DOA_sol$peak_width > PWMAX) * 3 + 1
      return(DOA_sol)
    }


    MUSIC %<>%
      dplyr::mutate(DOA = purrr::map2(DOA, retained_solution, \(x, y) check_PWMAX(x, y, PWMAX)))

    seasonder_cs_object %<>%
      seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  }
  return(seasonder_cs_object)
}

#' Execute the MUSIC Algorithm on a SeaSondeRCS Object
#'
#' This function performs the MUSIC (MUltiple SIgnal Classification) algorithm on a given
#' SeaSondeRCS object, executing a series of processing steps to extract direction-of-arrival (DOA) information
#' and other related metrics from the radar cross-spectrum data.
#'
#' @param seasonder_cs_object A \code{SeaSondeRCS} object that contains the radar cross-spectrum data
#'        and metadata. This object is modified in place to include the results of the MUSIC algorithm.
#'
#' @return A \code{SeaSondeRCS} object with updated MUSIC-related attributes. Specifically:
#' \itemize{
#'   \item Processing steps annotated with the MUSIC start and end points.
#'   \item Updated attributes and fields for covariance matrix computations, DOA estimations, and other MUSIC-related metrics.
#' }
#'
#' @details
#' The MUSIC algorithm is executed in a series of sequential steps:
#' \enumerate{
#'   \item Log the start of the MUSIC algorithm.
#'   \item Update the processing steps of the \code{SeaSondeRCS} object to include the MUSIC start text.
#'   \item Perform the following computations:
#'     \itemize{
#'       \item Compute the covariance matrix from the cross-spectrum data.
#'       \item Perform eigen decomposition on the covariance matrix.
#'       \item Compute the DOA functions using MUSIC-specific methods.
#'       \item Extract peaks from the DOA functions, corresponding to possible signal directions.
#'       \item Calculate the signal power matrix.
#'       \item Test for dual solutions and compute their proportions.
#'       \item Select the final set of DOAs from the computed data.
#'       \item Convert the selected DOAs to geographical coordinates (latitude and longitude).
#'     }
#'   \item Log the completion of the MUSIC algorithm.
#' }
#'
#' @seealso
#' \code{\link{seasonder_MUSICComputeCov}}: Compute the covariance matrix.
#' \code{\link{seasonder_MUSICCovDecomposition}}: Perform eigen decomposition of the covariance matrix.
#' \code{\link{seasonder_MUSICComputeDOAProjections}}: Compute the direction-of-arrival functions.
#' \code{\link{seasonder_MUSICExtractPeaks}}: Extract peaks from the DOA functions.
#' \code{\link{seasonder_MUSICComputeSignalPowerMatrix}}: Calculate the signal power matrix.
#' \code{\link{seasonder_MUSICTestDualSolutions}}: Test and analyze dual solutions in the DOA.
#' \code{\link{seasonder_MUSICComputePropDualSols}}: Compute proportions for dual solutions.
#' \code{\link{seasonder_MUSICSelectDOA}}: Select final DOA estimations.
#' \code{\link{seasonder_MUSIC_LonLat}}: Convert DOA estimations to geographical coordinates.
#'
#' @examples
#' \donttest{
#' # Prepare a SeaSondeRCS object with MUSIC data
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' cs_obj <- seasonder_initMUSICData(cs_obj)
#' cs_obj <- seasonder_initMUSICData(
#'  cs_obj,
#'  range_cells = c(rep(5,11), rep(4,11)),
#'  doppler_bins = c(c(669:679),c(674:684))
#' )
#' cs_obj <- seasonder_SeaSondeRCSMUSICInterpolateDoppler(cs_obj)
#' # Run the MUSIC algorithm
#' cs_obj <- seasonder_runMUSIC(cs_obj)
#' # Check the updated processing steps
#' print(seasonder_getSeaSondeRCS_ProcessingSteps(cs_obj))
#' }
#' @export
seasonder_runMUSIC <- function(seasonder_cs_object){

  # Log the start of the MUSIC algorithm.
  seasonder_logAndMessage("seasonder_runMUSIC: MUSIC algorithm started.", "info")
options <- seasonder_getSeaSondeRCS_MUSIC_options(seasonder_cs_object)
  
discard_low_SNR <- options$discard_low_SNR
discard_no_solution <- options$discard_no_solution

doppler_interpolation <- options$doppler_interpolation
  # Create a copy of the input object to store the results of the processing.
  out <- seasonder_cs_object

  out %<>% seasonder_setSeaSondeRCS_MUSIC_options(options)

  MUSIC_options <- out %>% seasonder_getSeaSondeRCS_MUSIC_options()

  # Update the processing steps to indicate the start of the MUSIC algorithm.
  out %<>% seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_MUSIC_start_step_text())

  # Set the Doppler interpolation level for the MUSIC algorithm
  out %<>% seasonder_setSeaSondeRCS_MUSIC_doppler_interpolation(doppler_interpolation)

  # Apply Doppler interpolation
  out %<>% seasonder_SeaSondeRCSMUSICInterpolateDoppler()


  # Compute the covariance matrix from the cross-spectrum data.
  out %<>% seasonder_MUSICComputeCov()


out %<>% seasonder_computeNoiseLevel(antenna = 1,smoothed= MUSIC_options$smoothNoiseLevel)
  out %<>% seasonder_computeNoiseLevel(antenna = 2,smoothed= MUSIC_options$smoothNoiseLevel)
  out %<>% seasonder_computeNoiseLevel(antenna = 3,smoothed= MUSIC_options$smoothNoiseLevel)

  out %<>% seasonder_computeSignalSNR()



  out %<>% seasonder_SNRCheck(discard_low_SNR = discard_low_SNR)

  # Perform eigen decomposition of the covariance matrix.
  out %<>% seasonder_MUSICCovDecomposition()

  # Compute the direction-of-arrival (DOA) functions using the MUSIC algorithm.
  out %<>% seasonder_MUSICComputeDOAProjections()

  # Extract peaks from the DOA functions, representing potential signal directions.
  out %<>% seasonder_MUSICExtractPeaks()


if(discard_no_solution){
  # Remove entries with no valid solutions from the MUSIC results.
  out %<>% seasonder_MUSIC_remove_no_solutions()
}

  # Compute the signal power matrix for the extracted peaks.
  out %<>% seasonder_MUSICComputeSignalPowerMatrix()

  # Test for the presence of dual solutions in the DOA estimations.
  out %<>% seasonder_MUSICTestDualSolutions()

  # Compute the proportion of dual solutions across the data.
  out %<>% seasonder_MUSICComputePropDualSols()

  # Select the final DOAs based on the computed data.
  out %<>% seasonder_MUSICSelectDOA()



  # Convert the selected DOAs into geographical coordinates (longitude and latitude).
  out %<>% seasonder_MUSIC_LonLat()


  if(!is.null(MUSIC_options$PPMIN)){
    out %<>% seasonder_checkPPMIN()
  }

  if(!is.null(MUSIC_options$PWMAX)){
    out %<>% seasonder_checkPWMAX()
  }

  # Update the processing steps to indicate the end of the MUSIC algorithm.
  out %<>% seasonder_setSeaSondeRCS_ProcessingSteps(SeaSondeRCS_MUSIC_end_step_text(out))

  # Log the completion of the MUSIC algorithm.
  seasonder_logAndMessage("seasonder_runMUSIC: MUSIC algorithm finished.", "info")

  # Return the updated object with the MUSIC results.
  return(out)
}


#' Run MUSIC Algorithm on FOR Data
#'
#' This function integrates the MUSIC (Multiple Signal Classification) algorithm into a SeaSondeRCS object that has First Order Regions (FOR) initialized.
#' It first applies Doppler interpolation to the cross-spectra data, then extracts the FOR boundaries for each range cell by transforming the negative 
#' and positive FOR Doppler bins into frequency values and subsequently mapping these frequencies back to Doppler bins. Finally, the function initializes 
#' the MUSIC data structure and invokes the full MUSIC algorithm to update the SeaSondeRCS object.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing cross-spectra data and FOR information.
#'   (Note: Although this parameter is specified as an argument in the documentation, the actual Doppler interpolation factor is retrieved from 
#'   the MUSIC options stored in the object.)
#'
#' @return A SeaSondeRCS object with its MUSIC data updated after applying Doppler interpolation, FOR extraction, and the complete MUSIC processing.
#'
#' @details
#' This function performs the following sequence of operations:
#' \enumerate{
#'   \item It retrieves the Doppler interpolation factor from the MUSIC options of the input object.
#'   \item It obtains the FOR data from the object using \code{seasonder_getSeaSondeRCS_FOR}.
#'   \item For each range cell in the FOR data:
#'     \enumerate{
#'       \item It processes the negative FOR bins by:
#'         \enumerate{
#'           \item Determining the frequency range corresponding to the negative bins via \code{seasonder_Bins2DopplerFreq}.
#'           \item Mapping these frequencies to new Doppler bin indices with \code{seasonder_MUSIC_DopplerFreq2Bins} and adjusting 
#'               the indices based on the interpolation factor.
#'         }
#'       \item It processes the positive FOR bins in an analogous manner.
#'       \item If valid Doppler bin indices are obtained, a data frame is created recording the range cell and Doppler bin information.
#'     }
#'   \item The function compiles the extracted FOR information from all range cells into a single data frame.
#'   \item It initializes the MUSIC data structure for the specified range cells and Doppler bins using \code{seasonder_initMUSICData}.
#'   \item Finally, it calls \code{seasonder_runMUSIC} to execute the MUSIC algorithm on the updated object.
#' }
#'
#' @importFrom purrr map compact
#' @importFrom dplyr bind_rows
#'
#' @examples
#' # Prepare a SeaSondeRCS object with MUSIC data (including FOR segments)
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' FOR <- seasonder_getSeaSondeRCS_FOR(cs_obj)
#' cs_obj <- seasonder_setSeaSondeRCS_FOR(cs_obj,FOR[4:5])
#' # Run MUSIC algorithm in FOR context
#' result <- seasonder_runMUSICInFOR(cs_obj)
#' # View processing steps
#' print(seasonder_getSeaSondeRCS_ProcessingSteps(result))
#' @export
#' @aliases seasonder_runMUSIC_in_FOR
seasonder_runMUSICInFOR <- seasonder_runMUSIC_in_FOR <- function(seasonder_cs_object){

  # Initialize the output as a copy of the input SeaSondeRCS object
  out <- seasonder_cs_object


doppler_interpolation <- seasonder_getSeaSondeRCS_MUSIC_options(out)$doppler_interpolation
  
  # Retrieve the First Order Regions (FOR) from the SeaSondeRCS object
  FOR <- seasonder_getSeaSondeRCS_FOR(seasonder_cs_object)

  # Iterate over each range cell in the FOR data
  FOR <- seq_len(length(FOR)) %>% purrr::map(\(range_cell) {
    o <- NULL
    doppler_bins <- integer(0)
    # Process negative FOR bins
    neg_bins <- FOR[[range_cell]]$negative_FOR
    if (length(neg_bins) > 0) {
      neg_range <- range(neg_bins)

      neg_range_freq <- seasonder_Bins2DopplerFreq(seasonder_cs_object, neg_range)
      new_neg_range_bins <- seasonder_MUSIC_DopplerFreq2Bins(out, neg_range_freq)
      new_neg_range_bins[1] <- new_neg_range_bins[1]-1*(doppler_interpolation-1)
      new_neg_range_bins[2] <- new_neg_range_bins[2]-1*(doppler_interpolation-1)
      if(new_neg_range_bins[1] != new_neg_range_bins[2]){
       doppler_bins <- c(doppler_bins, new_neg_range_bins[1]:new_neg_range_bins[2])
      }
    }

    # Process positive FOR bins
    pos_bins <- FOR[[range_cell]]$positive_FOR
    if (length(pos_bins) > 0) {
      pos_range <- range(pos_bins)
      pos_range_freq <- seasonder_Bins2DopplerFreq(seasonder_cs_object, pos_range)
      new_pos_range_bins <- seasonder_MUSIC_DopplerFreq2Bins(out, pos_range_freq)
      new_pos_range_bins[1] <- new_pos_range_bins[1]-1*(doppler_interpolation-1)
      new_pos_range_bins[2] <- new_pos_range_bins[2]-1*(doppler_interpolation-1)
      if(new_pos_range_bins[1] != new_pos_range_bins[2]){
        doppler_bins <- c(doppler_bins, new_pos_range_bins[1]:new_pos_range_bins[2])
      }

    }

    # Create a data frame if Doppler bins are found
    if (length(doppler_bins) > 0) {
      o <- data.frame(range_cell = range_cell, doppler_bin = doppler_bins)
    }

    return(o)
  }) %>% purrr::compact() %>% dplyr::bind_rows()
if(nrow(FOR) == 0){
  seasonder_logAndAbort("No valid FOR data found. Please run seasonder_computeFORs first.", calling_function = "seasonder_runMUSICInFOR")
  }
  # Initialize MUSIC data for the specified range cells and Doppler bins
  out %<>% seasonder_initMUSICData(range_cells = FOR$range_cell, doppler_bins = FOR$doppler_bin, NULL_MUSIC = nrow(FOR) == 0)

  # Run the MUSIC algorithm on the updated SeaSondeRCS object
  out %<>% seasonder_runMUSIC()

  # Return the updated SeaSondeRCS object
  return(out)


}


#### Utils ####

#' Convert MUSIC Bearings to Geographic Bearings
#'
#' This function converts MUSIC bearings (relative to the antenna) into geographic bearings using the antenna's bearing information from a `SeaSondeRAPM` object.
#'
#' @param bearings A list of numeric vectors containing MUSIC bearings in degrees. Each vector corresponds to a set of bearings relative to the antenna.
#' @param seasonder_apm_object A `SeaSondeRAPM` object containing the antenna's metadata, including the antenna's bearing.
#'
#' @details
#' The geographic bearing is calculated by:
#' 1. Multiplying the MUSIC bearings by -1 to invert their direction.
#' 2. Adjusting the angles to the range [0, 360) using modulo 360.
#' 3. Adding the antenna bearing to each value and wrapping the result to the range [0, 360) again using modulo 360.
#'
#' The formula for each bearing is: \eqn{geo_bearing = ((-1 * music_bearing \\%\\% 360) + antenna_bearing) \\%\\% 360}.
#'
#' @return A list of numeric vectors containing the geographic bearings in degrees.
#'
#' @seealso
#' - \code{\link{seasonder_getSeaSondeRAPM_AntennaBearing}}
#' - \code{\link[magrittr]{\%>\%}}
#' - \code{\link[purrr]{map}}
#'
#' @importFrom purrr map
#'
seasonder_MUSICBearing2GeographicalBearing <- function(bearings, seasonder_apm_object) {

  # Retrieve the antenna's bearing from the SeaSondeRAPM object
  antennaBearing <- seasonder_apm_object %>%
    seasonder_getSeaSondeRAPM_AntennaBearing()

  # Convert MUSIC bearings to geographic bearings
  # 1. Multiply bearings by -1 to invert their direction
  # 2. Normalize to [0, 360) using modulo 360
  # 3. Add the antenna bearing and normalize to [0, 360) again
  bearings %<>% purrr::map(\(angles) ((-1 * angles %% 360) + antennaBearing) %% 360)

  # Return the converted geographic bearings
  return(bearings)
}


#' Compute Geographic Coordinates from Origin, Distance, and Bearing
#'
#' This function calculates the geographic coordinates (latitude and longitude) for a given distance and bearing from a specified origin.
#'
#' @param origin_lon A numeric value representing the longitude of the origin point in decimal degrees.
#' @param origin_lat A numeric value representing the latitude of the origin point in decimal degrees.
#' @param dist A numeric value representing the distance from the origin in kilometers.
#' @param bearing A numeric vector of bearings (in degrees) indicating the direction from the origin.
#'
#' @details
#' The function uses the geodetic formulas provided by the `geosphere` package to compute the destination point based on:
#' - Origin longitude and latitude
#' - Distance in meters (converted from kilometers)
#' - Bearing in degrees
#'
#' The calculation employs the `geosphere::destPoint` function, which handles the spherical geometry of the Earth.
#'
#' @return A data frame with two columns:
#' \itemize{
#'   \item `lon`: The longitude of the computed geographic coordinates.
#'   \item `lat`: The latitude of the computed geographic coordinates.
#' }
#'
#' @seealso
#' \code{\link[geosphere]{destPoint}}
#'
#' @importFrom geosphere destPoint
#' @importFrom magrittr %>%
#'
#' @examples
#' # Example with a point at 100 km to the north of the origin
#' result <- seasonder_computeLonLatFromOriginDistBearing(-123.3656, 48.4284, 100, 0)
#' print(result)
#' @export
seasonder_computeLonLatFromOriginDistBearing <- function(origin_lon, origin_lat, dist, bearing) {

  # Calculate the geographic destination point
  # - Converts distance from kilometers to meters
  # - Uses geosphere::destPoint for spherical geometry calculations
  out <- geosphere::destPoint(c(origin_lon, origin_lat), bearing, dist * 1000) %>%
    # Convert the result to a data frame for easier handling
    as.data.frame()
out$lon[is.nan(out$lon)] <- NA_real_
out$lat[is.nan(out$lat)] <- NA_real_

  

  # Return the computed geographic coordinates
  return(out)

}




#' Map MUSIC Bearings to Geographic Coordinates
#'
#' This function calculates geographic coordinates (latitude and longitude) for each MUSIC detection based on the range and direction of arrival (DOA) bearings from a `SeaSondeRCS` object.
#'
#' @param seasonder_cs_object A `SeaSondeRCS` object containing MUSIC detection data.
#'
#' @details
#' This function performs the following operations:
#' 1. Retrieves MUSIC data and original geographic coordinates (latitude and longitude) from the `seasonder_cs_object`. If these coordinates are not available, the origin is derived from the associated Antenna Pattern (APM) data.
#' 2. Converts DOA bearings from MUSIC detections into geographic bearings using the APM object.
#' 3. Computes latitude and longitude for each MUSIC detection based on the range and geographic bearings using \code{\link{seasonder_computeLonLatFromOriginDistBearing}}
#' 4. Updates the `seasonder_cs_object` with the newly computed coordinates.
#'
#' @return A `SeaSondeRCS` object with updated MUSIC data, including geographic coordinates for each detection.
#'
#' @seealso
#' - \code{\link{seasonder_getSeaSondeRCS_MUSIC}}
#' - \code{\link{seasonder_getSeaSondeRCS_APM}}
#' - \code{\link{seasonder_MUSICBearing2GeographicalBearing}}
#' - \code{\link{seasonder_computeLonLatFromOriginDistBearing}}
#'

#' @importFrom purrr map map2
#'
 #'
 #' @export
 #' @examples
 #' # Create a SeaSondeRCS object for MUSIC example
 #' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
 #' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
 #' apm_obj  <- seasonder_readSeaSondeRAPMFile(apm_file)
 #' cs_obj   <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' cs_obj <- seasonder_initMUSICData(
#'  cs_obj,
#'  range_cells = c(rep(5,11), rep(4,11)),
#'  doppler_bins = c(c(669:679),c(674:684))
#' )
#' cs_obj <- seasonder_runMUSIC(cs_obj)
 #' updated_obj <- seasonder_MUSICLonLat(cs_obj)
 #' print(updated_obj)
 #' @aliases seasonder_MUSIC_LonLat
seasonder_MUSICLonLat <- seasonder_MUSIC_LonLat <- function(seasonder_cs_object) {

  # Retrieve MUSIC data from the SeaSondeRCS object
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # Retrieve original longitude and latitude; fallback to APM origin if unavailable
  longitude <- seasonder_getfLongitude(seasonder_cs_object)
  latitude <- seasonder_getfLatitude(seasonder_cs_object)
  if (is.null(longitude) || is.null(latitude) || abs(longitude) < 1e-10 || abs(latitude) < 1e-10) {
    # Extract the origin from the APM data
    origin <- seasonder_cs_object %>%
      seasonder_getSeaSondeRCS_APM() %>%
      seasonder_getSeaSondeRAPM_SiteOrigin()
    latitude <- origin[1]
    longitude <- origin[2]
  }

  # Extract range and DOA bearings from MUSIC data
  range <- MUSIC$range
  bearings <- MUSIC$DOA %>% purrr::map("bearing")

  # Retrieve the APM object from the SeaSondeRCS object
  seasonder_apm_object <- seasonder_cs_object %>%
    seasonder_getSeaSondeRCS_APM()

  # Convert MUSIC bearings to geographic bearings using the APM object
  bearings %<>% seasonder_MUSICBearing2GeographicalBearing(seasonder_apm_object)

  # Compute geographic coordinates for each MUSIC detection
  MUSIC$lonlat <- purrr::map2(range, bearings, \(dist, bear) {
    if (length(bear) > 0) {
      # Calculate geographic coordinates for valid detections
      seasonder_computeLonLatFromOriginDistBearing(longitude, latitude, dist = dist, bearing = bear)
    } else {
      # Return NA for invalid detections
      data.frame(lon = NA_real_, lat = NA_real_)
    }
  })

  # Update the SeaSondeRCS object with updated MUSIC data
  seasonder_cs_object %<>% seasonder_setSeaSondeRCS_MUSIC(MUSIC)

  # Return the updated SeaSondeRCS object
  return(seasonder_cs_object)
}



#' Export MUSIC Table from SeaSondeRCS Object
#'
#' This function generates a table containing detailed MUSIC detection data from a `SeaSondeRCS` object. The output table includes geographic coordinates, signal parameters, and other metadata for each MUSIC detection.
#'
#' @param seasonder_cs_object A `SeaSondeRCS` object containing MUSIC detection data and related metadata.
#'
#' @details
#' This function performs the following operations:
#' 1. Retrieves the timestamp (`nDateTime`) from the header of the `SeaSondeRCS` object. Defaults to `as.POSIXct(0)` if unavailable.
#' 2. Initializes an empty data frame with predefined columns.
#' 3. Retrieves MUSIC detection data, processes the Direction of Arrival (DOA) and geographic coordinates (`lonlat`), and unnests these fields.
#' 4. Converts MUSIC bearings to geographic bearings using the associated Antenna Pattern Matrix (APM) object.
#' 5. Computes additional metrics such as signal power in dB, signal-to-noise ratio (SNR), and DOA peak response in dB.
#' 6. Appends the timestamp to the table and reorders columns for clarity.
#'
#' @return A data frame with the following columns:
#' - `datetime`: Timestamp of the data.
#' - `longitude`: Geographic longitude of the detection.
#' - `latitude`: Geographic latitude of the detection.
#' - `range_cell`: Range cell number.
#' - `range`: Range in kilometers.
#' - `doppler_bin`: Doppler bin number.
#' - `doppler_freq`: Doppler frequency.
#' - `radial_velocity`: Radial velocity in m/s.
#' - `signal_power`: Signal power.
#' - `bearing`: Geographic bearing in degrees.
#' - `bearing_raw`: Original MUSIC bearing in degrees.
#' - `noise_level`: Noise level in dB.
#' - `signal_power_db`: Signal power in dB.
#' - `SNR`: Signal-to-noise ratio in dB.
#' - `DOA_peak_resp_db`: DOA peak response in dB.
#'
#' @seealso
#' - \code{\link{seasonder_getSeaSondeRCS_MUSIC}}
#' - \code{\link{seasonder_MUSICBearing2GeographicalBearing}}
#' - \code{\link{seasonder_getSeaSondeRAPM_AntennaBearing}}
#'
#' @importFrom dplyr select mutate
#' @importFrom tidyr unnest
#' @importFrom purrr map
#' @importFrom pracma Real
#'
#' @export
#'
#' @examples
#' \donttest{
#'   # Load sample CSS and APM files
#'   cs_file  <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#'   apm_file <- system.file("css_data/MeasPattern.txt",       package = "SeaSondeR")
#'   apm_obj  <- seasonder_readSeaSondeRAPMFile(apm_file)
#'   # Create SeaSondeRCS object with APM
#'   cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' FOR <- seasonder_getSeaSondeRCS_FOR(cs_obj)
#' cs_obj <- seasonder_setSeaSondeRCS_FOR(cs_obj,FOR[4:5])
#' 
#'   # Run MUSIC algorithm (in FOR context) if MUSIC data is available:
#'   cs_obj <- seasonder_runMUSICInFOR(cs_obj)
#'   # Export MUSIC table
#'   music_table <- seasonder_exportMUSICTable(cs_obj)
#'   print(music_table)
#' }
seasonder_exportMUSICTable <- function(seasonder_cs_object) {
  # Initialize globals for seasonder_exportMUSICTable
  range_cell <- doppler_bin <- freq <- radial_v <- DOA <- lonlat <- 
  signal_power <- signal_power_db <- noise_level <- lon <- lat <- 
  doppler_freq <- radial_velocity <- bearing <- bearing_raw <- SNR <- 
  DOA_peak_resp_db <- rlang::zap()
  # Retrieve timestamp from the header; default to POSIXct(0) if unavailable
  datetime <- seasonder_getSeaSondeRCS_headerField(seasonder_cs_object, "nDateTime") %||% as.POSIXct(0)

  # Initialize an empty data frame with predefined columns
  out <- data.frame(
    longitude = numeric(0),
    latitude = numeric(0),
    range_cell = integer(0),
    range = numeric(0),
    doppler_bin = integer(0),
    doppler_freq = numeric(0),
    radial_velocity = numeric(0),
    signal_power = numeric(0),
    bearing = numeric(0),
    noise_level = numeric(0),
    signal_power_db = numeric(0),
    SNR = numeric(0),
    DOA_peak_resp_db = numeric(0)
  )

  # Retrieve MUSIC data from the SeaSondeRCS object
  MUSIC <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  # If MUSIC data is non-empty, process it
  if (nrow(MUSIC) > 0) {

    # Select and rename relevant columns from MUSIC data
    out <- MUSIC %>% dplyr::select(
      range_cell, doppler_bin, range,
      doppler_freq = freq, radial_velocity = radial_v, DOA, lonlat
    )

    # Process and unnest DOA and lonlat columns
    out %<>% dplyr::mutate(
      DOA = purrr::map(DOA, \(DOA_sol) {
        # Ensure DOA_sol is a list with expected elements
        if (is.list(DOA_sol) && !is.null(DOA_sol$bearing) && length(DOA_sol$bearing) > 0) {
          data.frame(
            bearing = DOA_sol$bearing,
            signal_power = abs(diag(DOA_sol$P)),
            DOA_peak_resp_db = DOA_sol$peak_resp
          )
        } else {
          data.frame(
            bearing = NA_real_,
            signal_power = NA_real_,
            DOA_peak_resp_db = NA_real_
          )
        }
      })
    ) %>%
      tidyr::unnest(c(DOA, lonlat))

    # Retrieve APM object for bearing conversion
    seasonder_apm_object <- seasonder_cs_object %>%
      seasonder_getSeaSondeRCS_APM()

    # Preserve raw MUSIC bearing values
    out$bearing_raw <- out$bearing

    # Convert MUSIC bearings to geographic bearings
    out$bearing %<>% seasonder_MUSICBearing2GeographicalBearing(seasonder_apm_object) %>% unlist()

    # Retrieve noise level and compute additional metrics
    FOR_config <- seasonder_getSeaSondeRCS_FORConfig(seasonder_cs_object)
    out$noise_level <- FOR_config$NoiseLevel[out$range_cell] %>% magrittr::set_names(NULL)

    out %<>% dplyr::mutate(
      signal_power_db = self_spectra_to_dB(signal_power, seasonder_getReceiverGain_dB(seasonder_cs_object)),
      SNR = signal_power_db - noise_level
    )

    # Reorder and rename columns for final output
    out %<>% dplyr::select(
      longitude = lon, latitude = lat, range_cell, range,
      doppler_bin, doppler_freq, radial_velocity, signal_power,
      bearing, bearing_raw, noise_level, signal_power_db, SNR,
      DOA_peak_resp_db
    )
  }

  # Add timestamp as the first column
  out %<>% dplyr::mutate(datetime = datetime, .before = 1)

  # Return the processed table
  return(out)
}


#' Export MUSIC Table to CSV
#'
#' This function exports the MUSIC detection table from a `SeaSondeRCS` object to a CSV file.
#'
#' @param seasonder_cs_object A `SeaSondeRCS` object containing MUSIC detection data.
#' @param filepath A character string specifying the path to the output CSV file.
#'
#' @details
#' This function performs the following steps:
#' 1. Generates a MUSIC table using \code{seasonder_exportMUSICTable}.
#' 2. Converts the resulting table to a data frame.
#' 3. Writes the data frame to the specified CSV file using \code{data.table::fwrite}.
#'
#' @return The function returns \code{NULL} invisibly. The output is saved to the specified file.
#'
#' @seealso
#' - \code{\link{seasonder_exportMUSICTable}}
#' - \code{\link[data.table]{fwrite}}
#'
#' @importFrom data.table fwrite
#'
#' @examples
#'   # Prepare a SeaSondeRCS object for examples, including APM
#'   apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#'   apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#'   specs_path <- seasonder_defaultSpecsFilePath("CS")
#'   cs_obj <- seasonder_createSeaSondeRCS(
#'     system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR"),
#'     specs_path = specs_path,
#'     seasonder_apm_object = apm_obj
#'   )
#' cs_obj <- seasonder_initMUSICData(
#'  cs_obj,
#'  range_cells = c(rep(5,11), rep(4,11)),
#'  doppler_bins = c(c(669:679),c(674:684))
#' )
#' cs_obj <- seasonder_runMUSIC(cs_obj)
#'   # Export MUSIC table to a temporary CSV file
#'   tmpfile <- tempfile(fileext = ".csv")
#'   seasonder_exportCSVMUSICTable(cs_obj, tmpfile)
#'   print(tmpfile)
#' @export
seasonder_exportCSVMUSICTable <- function(seasonder_cs_object, filepath) {

  # Generate the MUSIC table from the SeaSondeRCS object
  table <- seasonder_exportMUSICTable(seasonder_cs_object) %>%
    # Convert the resulting table to a data frame
    as.data.frame()

  # Write the data frame to the specified CSV file
  data.table::fwrite(table, file = filepath)

  # Return NULL invisibly to indicate successful execution
  invisible(NULL)
}

seasonder_computeMDRJ <- function(music){

  out <- rep(0,nrow(music))

out <- out + as.integer(!music$P1_check) * 1

out <- out + as.integer(!music$P2_check) * 2 * as.integer(music$P0_check)

out <- out + as.integer(!music$P3_check) * 4 * as.integer(music$P0_check)

out <- out + as.integer(!music$P4_check) * 8 * as.integer(music$P0_check)

out <- out + as.integer(!music$P0_check) * 16


  return(out)

}

#' Export Range Information from a SeaSondeRCS Object
#'
#' This function computes and exports range-related information based on the MUSIC data stored in a SeaSondeRCS object.
#' The output table includes range cell identifiers, range values, noise levels (for each antenna), first-order region
#' (FOR) boundaries, and counts of detections classified as single or dual solutions.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC data and associated metadata.
#'
#' @return A data frame with the following columns:
#' \describe{
#'   \item{SPRC}{Range cell identifier.}
#'   \item{RNGC}{Range (in appropriate units).}
#'   \item{NF01}{Noise level (in dB) for antenna 1.}
#'   \item{NF02}{Noise level (in dB) for antenna 2.}
#'   \item{NF03}{Noise level (in dB) for antenna 3.}
#'   \item{ALM1}{Lower FOR boundary (after Doppler interpolation).}
#'   \item{ALM2}{Upper FOR boundary (after Doppler interpolation) for the first boundary set.}
#'   \item{ALM3}{Lower FOR boundary (after Doppler interpolation) for the second boundary set.}
#'   \item{ALM4}{Upper FOR boundary (after Doppler interpolation) for the second boundary set.}
#'   \item{NVSC}{Count of detections classified as "single".}
#'   \item{NVDC}{Count of detections classified as "dual".}
#'   \item{NVAC}{Total adjusted count (NVSC plus twice NVDC).}
#' }
#'
#' @details
#' The function performs the following operations:
#' \enumerate{
#'   \item Extracts key fields from the MUSIC data: range cell, range, and the retained solution type.
#'   \item Aggregates counts of detections classified as single versus dual solutions.
#'   \item Retrieves noise levels (in dB) for each antenna.
#'   \item Obtains FOR boundaries using a dedicated export function and adjusts them based on the Doppler
#'         interpolation factor.
#'   \item Merges the aggregated detection counts, noise levels, and FOR boundaries by range cell.
#'   \item Selects and reorders the output columns.
#' }
#'
#' @examples
#' \donttest{
#' # Prepare a SeaSondeRCS object with MUSIC data
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' FOR <- seasonder_getSeaSondeRCS_FOR(cs_obj)
#' cs_obj <- seasonder_setSeaSondeRCS_FOR(cs_obj,FOR[4:5])
#' # Run MUSIC algorithm to populate MUSIC data
#' cs_obj <- seasonder_runMUSICInFOR(cs_obj)
#' range_info <- seasonder_exportRangeInfo(cs_obj)
#' head(range_info)
#' }
#' @export
seasonder_exportRangeInfo <- function(seasonder_cs_object){
  # Initialize globals for seasonder_exportRangeInfo
  range_cell <- retained_solution <- MSEL <- is_single <- is_dual <- NVSC <- NVDC <- SPRC <- RNGC <- rlang::zap()
  cols <- c("SPRC", "RNGC", "NF01", "NF02", "NF03", "ALM1", "ALM2", "ALM3", "ALM4", "NVSC", "NVDC", "NVAC")

  rm <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object) %>% dplyr::select(SPRC = range_cell, RNGC = range, MSEL = retained_solution)

  rc <- rm %>% dplyr::mutate(is_single = MSEL == "single", is_dual = MSEL == "dual") %>% dplyr::summarise(NVSC = sum(is_single), NVDC = sum(is_dual), NVAC = NVSC + NVDC * 2,.by = c(SPRC, RNGC))

  NF01 <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_NoiseLevel(dB = TRUE, antenna = 1)
  NF02 <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_NoiseLevel(dB = TRUE, antenna = 2)
  NF03 <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_NoiseLevel(dB = TRUE, antenna = 3)
  Noise <- data.frame(NF01 = NF01, NF02 = NF02, NF03 = NF03) %>% dplyr::mutate(SPRC = dplyr::row_number())

  FOR <-   seasonder_SeaSondeRCSExportFORBoundaries(seasonder_cs_object) %>% magrittr::set_colnames(c("SPRC","ALM1","ALM2","ALM3","ALM4"))

# FOR <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_FOR()
# FOR %<>% purrr::map_dfr(\(range_d) c(range(range_d$negative_FOR),range(range_d$positive_FOR)) %>% as.list() %>% magrittr::set_names(c("ALM1", "ALM2", "ALM3", "ALM4")) %>% as.data.frame()) %>% dplyr::mutate(SPRC = dplyr::row_number(), .before = 1)
interp <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation()
FOR %<>% dplyr::mutate(dplyr::across(dplyr::all_of(c("ALM1","ALM2","ALM3","ALM4")), \(x) (x-1)* interp ))


out <- rc %>% dplyr::left_join(FOR, by = "SPRC") %>% dplyr::left_join(Noise, by = "SPRC")


out %<>% dplyr::select(dplyr::all_of(cols))

attr(out, "radial_metrics") <- rm

return(out)
}

#' Export Radial Metrics from a SeaSondeRCS Object
#'
#' This function extracts and formats radial metrics from a SeaSondeRCS object for export.  
#' It processes the MUSIC table, computes various spectral metrics, applies antenna pattern corrections,  
#' and combines the results into a final data frame formatted according to predefined column specifications.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC detection data and related metadata.
#' @param AngSeg An optional list of angular segments to be applied to the vector flag field (VFLG).  
#'        Each element should be a numeric vector of length 3 defining a segment. Default is an empty list.
#'
#' @return A data frame with 34 columns containing the computed radial metrics. The columns include geographic  
#'         coordinates, velocity components, range, bearing information, signal power metrics, noise thresholds,  
#'         and computed spectral parameters.
#'
#' @details
#' The function proceeds as follows:
#' \enumerate{
#'   \item Retrieves the MUSIC table using \code{seasonder_getSeaSondeRCS_MUSIC} and the associated APM object.
#'   \item Defines a template row with 34 predefined columns, initializing most numeric values to NA, except for  
#'         specific defaults such as \code{MSA1}, \code{MDA1}, and \code{MDA2} (set to 1440L).
#'   \item Copies basic numeric fields and computes additional fields from the MUSIC table, such as the radial  
#'         velocity (scaled by 100), range, range cell, doppler cell (shifted by -1), eigenvalue ratio, signal power  
#'         ratio, and offset power ratio.
#'   \item Computes the metric \code{MDRJ} by applying the function \code{seasonder_computeMDRJ} on the MUSIC row.
#'   \item Extracts eigen decomposition results from each MUSIC row to populate the eigenvalue fields (\code{MEI1},  
#'         \code{MEI2}, \code{MEI3}).
#'   \item Processes the DOA solutions stored in each MUSIC row:  
#'       - For solutions retained as "single", geographic bearing corrections are applied to populate \code{MSA1}.  
#'       - For dual-bearing solutions, the first two elements of the DOA bearings populate \code{MDA1} and  
#'         \code{MDA2}, respectively.
#'   \item Computes additional spectral metrics such as the self-spectra conversion to dB (fields \code{MA1S},  
#'         \code{MA2S}, and \code{MA3S}) after subtracting the noise level (obtained for each antenna).
#'   \item Based on the retained solution type (either "single" or "dual"), assigns location data (if available),  
#'         sets selection flags, and computes additional output metrics (e.g., \code{PPFG} and \code{PWFG}).
#'   \item Finally, all rows are combined into a data frame. If angular segments are provided, additional modifications  
#'         to the vector flag (\code{VFLG}) are applied.
#' }
#'
#' @examples
#' \donttest{
#' # Prepare a SeaSondeRCS object with MUSIC data
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' FOR <- seasonder_getSeaSondeRCS_FOR(cs_obj)
#' cs_obj <- seasonder_setSeaSondeRCS_FOR(cs_obj,FOR[4:5])
#' # Run MUSIC algorithm to populate MUSIC data
#' cs_obj <- seasonder_runMUSICInFOR(cs_obj)
#' radial_metrics <- seasonder_exportRadialMetrics(cs_obj, AngSeg = list(c(5, 30, 60)))
#' head(radial_metrics)
#' }
#' @export
seasonder_exportRadialMetrics <- function(seasonder_cs_object, AngSeg = list()) {
  # Initialize globals for seasonder_exportRadialMetrics
  DOA_solutions <- radial_v <- doppler_bin <- diag_off_diag_power_ratio <- length_dual <- 
  MPKR <- MDP1 <- MDP2 <- MDR1 <- MDR2 <- MDW1 <- MDW2 <- cov <- 
  SPRC <- MAXS <- Noise_3 <- Noise_2 <- Noise_1 <- .id <- lonlat <- DOA <- 
  MSA1 <- MDA1 <- MDA2 <- retained_solution <- BEAR <- VFLG <- PPFG <- PWFG <- 
  VELO <- HEAD <- length_single <- MSP1 <- MSW1 <- MSR1 <- rlang::zap()
  # Obtain the MUSIC table using the function seasonder_getSeaSondeRCS_MUSIC from the global environment. This allows the function to be overridden using local_redefs.
  music <- seasonder_getSeaSondeRCS_MUSIC(seasonder_cs_object)

  seasonder_apm_obj <- seasonder_getSeaSondeRCS_APM(seasonder_cs_object)

  # Define the expected 34 columns
  cols <- c("LOND", "LATD", "VELU", "VELV", "VFLG", "RNGE", "BEAR", "VELO", "HEAD",
            "SPRC", "SPDC", "MSEL", "MSA1", "MDA1", "MDA2", "MEGR", "MPKR", "MOFR",
            "MP13", "MP23",
            "MSP1", "MDP1", "MDP2", "MSW1", "MDW1", "MDW2", "MSR1", "MDR1", "MDR2",
            "MA1S", "MA2S", "MA3S", "MEI1", "MEI2", "MEI3", "MDRJ","PPFG","PWFG")

  # List to collect output rows
  

  receiver_gain <- seasonder_getReceiverGain_dB(seasonder_cs_object)

out <- data.frame(matrix(NA, nrow = 1, ncol = length(cols))) %>% magrittr::set_colnames(cols)
if(nrow(music) > 0){
  out <- music

  out <- out %>% 
  dplyr::rename(
    c(
      RNGE = "range",
      SPRC = "range_cell",
      MEGR = "eigen_values_ratio",
      MPKR = "signal_power_ratio"
      )
    ) %>% 
  dplyr::mutate(
    .id = dplyr::row_number(),
    length_single = purrr::map_int(DOA_solutions, \(x) length(which(!is.na(x$single$bearing)))),
    length_dual = purrr::map_int(DOA_solutions, \(x) length(which(!is.na(x$dual$bearing)))),
    VELO = radial_v * 100,
    SPDC = doppler_bin -1,
    MOFR = tidyr::replace_na(diag_off_diag_power_ratio, 0),
    MEI1 = purrr::map_dbl(eigen, \(x) x$values[1]),
    MEI2 = purrr::map_dbl(eigen, \(x) x$values[2]),
    MEI3 = purrr::map_dbl(eigen, \(x) x$values[3]),
    MSA1 = purrr::map_dbl(DOA_solutions, \(x) x$single$bearing %||% NA),
    MDA1 = purrr::map_dbl(DOA_solutions, \(x) x$dual$bearing[1] %||% NA),
    MDA2 = purrr::map_dbl(DOA_solutions, \(x) x$dual$bearing[2] %||% NA),
    MSA1 = dplyr::case_when(!is.na(MSA1) ~ unlist(seasonder_MUSICBearing2GeographicalBearing(MSA1,seasonder_apm_obj)), TRUE ~ 1440),
    MDA1 = dplyr::case_when(!is.na(MDA1) ~ unlist(seasonder_MUSICBearing2GeographicalBearing(MDA1,seasonder_apm_obj)), TRUE ~ 1440),
    MDA2 = dplyr::case_when(!is.na(MDA2) ~ unlist(seasonder_MUSICBearing2GeographicalBearing(MDA2,seasonder_apm_obj)), TRUE ~ 1440),
    MSR1 = purrr::map_dbl(DOA_solutions, \(x) 10^(x$single$peak_resp/10)),
    MDR1 = purrr::map_dbl(DOA_solutions, \(x) 10^(x$dual$peak_resp[1]/10)),
    MDR2 = purrr::map_dbl(DOA_solutions, \(x) 10^(x$dual$peak_resp[2]/10)),
    MSW1 = purrr::map_dbl(DOA_solutions, \(x) x$single$peak_width),
    MDW1 = purrr::map_dbl(DOA_solutions, \(x) x$dual$peak_width[1]),
    MDW2 = purrr::map_dbl(DOA_solutions, \(x) x$dual$peak_width[2]),
    MSP1 = purrr::map2_dbl(DOA_solutions,length_single, \(x,y) ifelse(y==1, 10*log10(abs(x$single$P)),NA)),
    MDP1 = purrr::map2_dbl(DOA_solutions,length_dual, \(x,y) ifelse(y>0,10*log10(abs(x$dual$P[1,1])),NA)),
    MDP2 = purrr::map2_dbl(DOA_solutions,length_dual, \(x,y) ifelse(y>1,10*log10(abs(x$dual$P[2,2])),NA)),
    MPKR = tidyr::replace_na(MPKR, 0),
    MSP1 = tidyr::replace_na(MSP1, 0),
    MDP1 = tidyr::replace_na(MDP1, 0),
    MDP2 = tidyr::replace_na(MDP2, 0),
    MSR1 = tidyr::replace_na(MSR1, 0),
    MDR1 = tidyr::replace_na(MDR1, 0),
    MDR2 = tidyr::replace_na(MDR2, 0),
    MSW1 = tidyr::replace_na(MSW1, 0),
    MDW1 = tidyr::replace_na(MDW1, 0),
    MDW2 = tidyr::replace_na(MDW2, 0),
    MP13 = purrr::map_dbl(cov, \(x) Arg(x[1,3])*180/pi),
    MP23 = purrr::map_dbl(cov, \(x) Arg(x[2,3])*180/pi),
    Noise_3 = seasonder_getSeaSondeRCS_NoiseLevel(seasonder_cs_object, dB = TRUE, antenna = 3)[SPRC],
    Noise_2 = seasonder_getSeaSondeRCS_NoiseLevel(seasonder_cs_object, dB = TRUE, antenna = 2)[SPRC],
    Noise_1 = seasonder_getSeaSondeRCS_NoiseLevel(seasonder_cs_object, dB = TRUE, antenna = 1)[SPRC],
    MAXS = purrr::map(cov, \(x) self_spectra_to_dB(c(x[1,1], x[2,2], x[3,3]), receiver_gain)),
    MAS3 = purrr::map2_dbl(MAXS, Noise_3, \(x,y, cs=seasonder_cs_object) x[3]- y),
    MAS2 = purrr::map2_dbl(MAXS, Noise_2, \(x,y, cs=seasonder_cs_object) x[2]- y),
    MAS1 = purrr::map2_dbl(MAXS, Noise_1, \(x,y, cs=seasonder_cs_object) x[1]- y)
  )

  out$MDRJ <- seasonder_computeMDRJ(out)

  sol_out <- out %>% 
  dplyr::mutate(data = purrr::pmap(list(.id,lonlat, DOA, MSA1, MDA1, MDA2, retained_solution), \(id,ll,d,ms,md1,md2,sol){
    
    # browser(expr = sol == "single")
    if(sol %in% c("single","dual")){
      o <- list(
        .id = id,
        LOND = NA_real_,
          LATD = NA_real_,
          MSEL = ifelse(sol == "single", 1L, 2L),
          BEAR = ifelse(sol == "single", ms, md1),
          PPFG = d$PPFG[1],
          PWFG = d$PWFG[1]
      )
      if(!is.null(ll) && nrow(ll)>0){

          o$LOND <- ll$lon[1]
          o$LATD <- ll$lat[1]
        
      }
      if (sol == "dual") {
        o2 <- list(
          .id = id,
          LOND = NA_real_,
          LATD = NA_real_,
          MSEL = 3L,
          BEAR = md2,
          PPFG = d$PPFG[2],
          PWFG = d$PWFG[2]
        )
        if(!is.null(ll) && nrow(ll)>1){

            o2$LOND <- ll$lon[2]
            o2$LATD <- ll$lat[2]
          
        }
        o <- c(list(o), list(o2))
       
      }
    }else {
       o <- list(.id = id, LOND = NA_real_, LATD = NA_real_,MSEL= NA_integer_,BEAR = NA_real_, PPFG = NA_real_, PWFG = NA_real_)
    }
    return(o)
  } ) 
  ) %>% dplyr::pull("data")
  
  sol_out <- sol_out %>% purrr::reduce(\(so_far,x){
    if(length(x) == 2){
      so_far <- append(so_far,x)
    }else{
      so_far <- append(so_far, list(x))
    }
    so_far
  },.init = list()) %>% purrr::transpose() %>% purrr::map(unlist) %>% as.data.frame()

  
  out <- out %>% 
  dplyr::left_join(sol_out, by = ".id") %>% dplyr::select(-.id)

  out <- out %>% dplyr::mutate(HEAD = (BEAR -180) %% 360)

  
  # If no rows, return empty data.frame with correct columns
  if (nrow(out) > 0) {
    # result <- do.call(rbind, lapply(out_rows, as.data.frame))
    out <- out %>% dplyr::mutate( VFLG = VFLG + 4096L * as.integer(! PPFG %in% c(1,9) | !PWFG %in% c(1,9)),
                               VELU = VELO * sin(HEAD * pi / 180),
                               VELV = VELO * cos(HEAD * pi / 180)#,
                             
                               )


    if(length(AngSeg) > 0){
      out <- purrr::reduce(AngSeg,\(result_so_far,seg){
        if(length(seg) == 3 && seg[1] %in% out$SPRC && seg[2] <= seg[3]){
          result_so_far <- result_so_far %>% dplyr::mutate(VFLG = VFLG + 128L * as.integer(SPRC ==seg[1] & dplyr::between(BEAR,seg[2], seg[3])))
        }
        return(result_so_far)
      },.init = out)
    }

  } else {
    out <- data.frame(matrix(ncol = length(cols), nrow = 0))
    colnames(out) <- cols
  }

  out <- out %>% dplyr::select(dplyr::all_of(cols))
}
  return(out)
}


#' Export LLUV Radial Metrics to a File
#'
#' This function extracts radial metrics from a SeaSondeRCS object and formats them for export using
#' defined mustache templates. The formatted output, which includes MUSIC parameters, antenna pattern corrections,
#' noise thresholds, and other spectral metrics, is written to a specified file. Additionally, the function
#' returns the computed radial metrics as a data frame.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing MUSIC detection data and related metadata.
#' @param LLUV_path A character string specifying the output file path for the LLUV radial metrics.
#' @param ... Additional arguments passed to \code{seasonder_exportRadialMetrics}.
#'
#' @details
#' The function performs the following steps:
#' \enumerate{
#'   \item Retrieves the radial metrics from the SeaSondeRCS object using \code{seasonder_exportRadialMetrics}.
#'   \item Obtains MUSIC parameters and antenna pattern attributes from the object.
#'   \item Formats numeric values using predefined formats for each column.
#'   \item Renders a data template (from "LLUV_RDM1_data.mustache") with the formatted radial metrics.
#'   \item Generates a deterministic UUID from the rendered data.
#'   \item Renders an overall LLUV template (from "LLUV_RDM1.mustache") that incorporates the radial parameters,
#'       formatted data, header information, and the generated UUID.
#'   \item Writes the rendered LLUV content to the file specified by \code{LLUV_path}.
#' }
#'
#' @return Invisibly returns a data frame containing the radial metrics used in the export.
#'
#' @examples
#' # Prepare a SeaSondeRCS object with MUSIC data
#' apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#' apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#' cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#' cs_obj <- seasonder_createSeaSondeRCS(cs_file, seasonder_apm_object = apm_obj)
#' FOR <- seasonder_getSeaSondeRCS_FOR(cs_obj)
#' cs_obj <- seasonder_setSeaSondeRCS_FOR(cs_obj,FOR[4:5])
#' # Optionally, run MUSIC in FOR context to populate MUSIC data
#' cs_obj <- seasonder_runMUSICInFOR(cs_obj)
#' radial_metrics <- seasonder_exportLLUVRadialMetrics(cs_obj, tempfile(fileext = ".ruv"))
#' head(radial_metrics)
#' @export
seasonder_exportLLUVRadialMetrics <- function(seasonder_cs_object, LLUV_path,...) {
  
  apm_object <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_APM()

  radial_metrics <- seasonder_exportRadialMetrics(seasonder_cs_object,...)

  MUSIC_params <- seasonder_cs_object %>% seasonder_getSeaSondeRCS_MUSIC_parameters() %>% magrittr::extract(1:3)
  # Removed unused variable 'header'
  APM_attributes <- attributes(apm_object)

  sprintf_vector <- function(x, format, sep = " ") {
  vec_format <- rep(format, length(x)) %>% paste0(collapse = sep)
  do.call(sprintf, c(list(vec_format), as.list(x)))
  }

  get_col_format <- function(col_name) {

  col_formats <- list(
    # Coordinates (in degrees)
    list(cols = c("LOND"), format = "%12.7f"),
    list(cols = c("LATD"), format = "%12.7f"),
    # Velocity components (cm/s)
    list(cols = c("VELU", "VELV"), format = "%9.3f"),
    list(cols = c("VELO"), format = "%11.3f"),
    # Vector code (integer)
    list(cols = c("VFLG"), format = "%11d"),
    # Distances (km)
    list(cols = c("RNGE"), format = "%10.4f"),
    # Heading angle (e.g., Bearing in degrees)
    list(cols = c("BEAR"), format = "%8.1f"),
    # Direction (e.g., Head)
    list(cols = c("HEAD"), format = "%10.1f"),
    # Associated cells (e.g., RngCell, DopCell and selection flag)
    list(cols = c("SPRC"), format = "%10d"),
    list(cols = c("SPDC"), format = "%9d"),
    list(cols = c("MSEL"), format = "%6d  "),
    # Measurements associated with MusicSngl/MusicDual (numeric values with one decimal)
    list(cols = c("MSA1", "MDA1", "MDA2"), format = "%9.1f "),
    # Eigen ratio (Eigen Ratio)
    list(cols = c("MEGR"), format = "%14.4f"),
    # Power ratio (Power Ratio)
    list(cols = c("MPKR"), format = "%13.5f"),
    # Offset ratio (Off Ratio)
    list(cols = c("MOFR"), format = "%13.6f"),
    # Phases A13 and A23 (angles)
    list(cols = c("MP13", "MP23"), format = "%8.1f "),
    # Columns associated with Pwr, Pk Width, Peak Resp, S/N, etc. (displayed with one decimal)
    list(cols = c("MSP1", "MDP1", "MDP2"), format = "%10.1f"),
    list(cols = c("MSW1", "MDW1", "MDW2"), format = "%9.1f "),
    list(cols = c("MSR1", "MDR1", "MDR2"), format = "%10.1f"),
    list(cols = c("MA1S", "MA2S", "MA3S"), format = "%8.1f  "),
    # Very small values in scientific notation
    list(cols = c("MEI1", "MEI2", "MEI3"), format = "%14.5e"),
    # Count columns (peaks, rejections, etc.)
    list(cols = c("MDRJ", "PPFG", "PWFG"), format = "%6d  ")
  )

  fmt <- purrr::keep(col_formats, \(fmt) col_name %in% fmt$cols)
  out <- NULL
  if(length(fmt) > 0){
    out <- fmt %>% magrittr::extract2(1) %>% purrr::pluck("format")
  }

  return(out)
  }

  radial_metrics_fmt <- radial_metrics %>% dplyr::mutate(dplyr::across(dplyr::everything(), \(x) sprintf(get_col_format(dplyr::cur_column()), x)))

  radial_metrics_fmt <- as.list(radial_metrics_fmt) %>% purrr::transpose()
  
  data_string <- radial_metrics_fmt %>% purrr::map_chr(\(x) paste0(x[c("LOND","LATD","VELU","VELV","VFLG","RNGE","BEAR","VELO","HEAD","SPRC","SPDC","MSEL","MSA1","MDA1","MDA2","MEGR","MPKR","MOFR","MP13","MP23","MSP1","MDP1","MDP2","MSW1","MDW1","MDW2","MSR1","MDR1","MDR2","MA1S","MA2S","MA3S","MEI1","MEI2","MEI3","MDRJ","PPFG","PWFG")], collapse = "")) %>% paste0(collapse = "\n")
  # Render the data template from radial_metrics_fmt
  # data_string <- whisker::whisker.render(template_data, data= list(data=radial_metrics_fmt))

  # Calculate UUID_data deterministically using data_string as seed
  UUID_data <- toupper(uuid::UUIDfromName("80113bfb-8db7-018b-011d-fe529c9ec978",data_string))

  # NEW: Compute range info table using tableStart = "2"
  rangeInfo <- seasonder_exportCTFRangeInfo_string(seasonder_cs_object, tableStart = "2")

  data <- list(
  RadialMusicParameters = sprintf_vector(MUSIC_params, "%0.3f", " "),
  ncols = ncol(radial_metrics),
  nrows = nrow(radial_metrics),
  PatternPhaseCorrections = sprintf_vector(APM_attributes$PhaseCorrections, "%0.2f", " "),
  PatternAmplitudeCorrections = sprintf_vector(APM_attributes$AmplitudeFactors, "%0.4f", " "),
  RadialBraggNoiseThreshold = sprintf("%0.3f", seasonder_getFOR_noisefact(seasonder_cs_object)),
  RadialBraggPeakNull = sprintf("%0.3f", seasonder_getFOR_fdown(seasonder_cs_object)),
  RadialBraggPeakDropOff = sprintf("%0.3f", seasonder_getFOR_flim(seasonder_cs_object)),
  data = data_string,
  TimeStamp = format(as.POSIXct(seasonder_getSeaSondeRCS_headerField(seasonder_cs_object, "nDateTime"), origin = "1970-01-01"), "%Y %m %d  %H %M %S"),
  TransmitCenterFreqMHz = sprintf("%0.6f", seasonder_getCenterFreqMHz(seasonder_cs_object)),
  TransmitBandwidthKHz = sprintf("%0.6f",
                   -1^(seasonder_getSeaSondeRCS_headerField(seasonder_cs_object, "bSweepUp") == 0) *
                   seasonder_getSeaSondeRCS_headerField(seasonder_cs_object, "fBandwidthKHz")),
  TransmitSweepRateHz = sprintf("%0.6f", seasonder_getSeaSondeRCS_headerField(seasonder_cs_object, "fRepFreqHz")),
  RangeResolutionKMeters = sprintf("%0.6f", seasonder_getSeaSondeRCS_headerField(seasonder_cs_object, "fRangeCellDistKm")),
  Site = seasonder_getSeaSondeRCS_headerField(seasonder_cs_object, "nSiteCodeName"),
  TimeZone = seasonder_getSeaSondeRCS_headerField(seasonder_cs_object, "szTimeZone"),
  fHoursFromUTC = sprintf("%+0.3f", seasonder_getSeaSondeRCS_headerField(seasonder_cs_object, "fHoursFromUTC")),
  TimeCoverage = sprintf("%0.3f", seasonder_getSeaSondeRCS_headerField(seasonder_cs_object, "nCoverMinutes")),
  Origin = paste(APM_attributes$SiteOrigin, collapse = " "),
  UUID = UUID_data,
  PatternUUID = APM_attributes$FileID,
  AntennaBearing = sprintf("%0.1f", APM_attributes$AntennaBearing),
  PatternType = APM_attributes$Type %||% tools::file_path_sans_ext(APM_attributes$FileName),
  PatternDate = format(as.POSIXct(APM_attributes$PatternDate), "%Y %m %d  %H %M %S"),
  PatternResolution = sprintf("%0.1f", APM_attributes$PatternResolution),
  PatternSmoothing = sprintf("%d", APM_attributes$Smoothing),
  RangeStart = 1,
  RangeEnd = seasonder_getnRangeCells(seasonder_cs_object),
  RangeCells = seasonder_getnRangeCells(seasonder_cs_object),
  DopplerInterpolation = sprintf("%d", seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation(seasonder_cs_object)),
  DopplerCells = sprintf(
    "%d", 
    seasonder_getnDopplerCells(seasonder_cs_object) * 
    seasonder_getSeaSondeRCS_MUSIC_doppler_interpolation(seasonder_cs_object)
  ),
  BraggSmoothingPoints = sprintf("%d", seasonder_getFOR_nsm(seasonder_cs_object)),
  CurrentVelocityLimit = sprintf("%0.1f", seasonder_getFOR_currmax(seasonder_cs_object)),
  range_info_table = rangeInfo$out_str,
  ProcessedTimeStamp = format(Sys.time(), "%Y %m %d  %H %M %S"),
  SeaSondeRVersion = utils::packageVersion("SeaSondeR")
  )

  template <- system.file("templates", "LLUV_RDM1.mustache", package = "SeaSondeR") %>%
  readLines() %>% paste0(collapse = "\n")

  LLUV <- whisker::whisker.render(template, data = data)

  LLUV %>% writeLines(LLUV_path)

  return(radial_metrics)
}



seasonder_exportCTFRangeInfo_string <- function(seasonder_cs_object, tableStart = "") {
  range_info <- seasonder_exportRangeInfo(seasonder_cs_object)
  n <- nrow(range_info)
  tableRows <- n
  template_path <- system.file("templates", "music_range_template.mustache", package = "SeaSondeR")
  template <- paste(readLines(template_path, warn = FALSE), collapse = "\n")
  data_list <- list(
    TableRows = tableRows,
    TableStart = tableStart,
    rows = lapply(1:n, function(i) {
      row_data <- range_info[i, ]
      row_str <- sprintf("%% %4s %9s %8s %8s %8s %8s %8s %8s %8s %8s %8s %8s",
        formatC(row_data$SPRC, width = 4, flag = " "),
        formatC(row_data$RNGC, format = "f", digits = 4, width = 9),
        formatC(row_data$NF01, format = "f", digits = 1, width = 8),
        formatC(row_data$NF02, format = "f", digits = 1, width = 8),
        formatC(row_data$NF03, format = "f", digits = 1, width = 8),
        formatC(row_data$ALM1, format = "d", width = 8),
        formatC(row_data$ALM2, format = "d", width = 8),
        formatC(row_data$ALM3, format = "d", width = 8),
        formatC(row_data$ALM4, format = "d", width = 8),
        formatC(row_data$NVSC, format = "d", width = 8),
        formatC(row_data$NVDC, format = "d", width = 8),
        formatC(row_data$NVAC, format = "d", width = 8)
      )
      list(row = row_str)
    })
  )
  out_str <- whisker::whisker.render(template, data = data_list)
  return(list(out_str = out_str, range_info = range_info))
}

#' Export CTF Range Information to a File
#'
#' This function writes the formatted CTF range information, generated from a SeaSondeRCS object, to a specified file.
#'
#' @param seasonder_cs_object A SeaSondeRCS object containing the relevant MUSIC processing data.
#' @param file A character string specifying the output file path where the range information will be written.
#' @param tableStart A character string to prepend to the table output. Defaults to an empty string.
#'
#' @details
#' The function internally calls \code{seasonder_exportCTFRangeInfo_string} to obtain a formatted string of range information.
#' It then writes this output string to the specified file. Additionally, it returns the extracted range information
#' invisibly, allowing further processing if necessary.
#'
#' @return Invisibly returns a data frame containing the range information.
#'
#' @examples
#'   # Prepare a SeaSondeRCS object with valid data
#'   apm_file <- system.file("css_data/MeasPattern.txt", package = "SeaSondeR")
#'   apm_obj <- seasonder_readSeaSondeRAPMFile(apm_file)
#'   cs_file <- system.file("css_data/CSS_TORA_24_04_04_0700.cs", package = "SeaSondeR")
#'   cs_obj <- seasonder_createSeaSondeRCS(
#'     cs_file,
#'     seasonder_apm_object = apm_obj
#'   )
#'   # Export CTF range information to a temporary text file
#'   range_info <- seasonder_exportCTFRangeInfo(cs_obj, tempfile(fileext = ".txt"))
#'
#' @export
seasonder_exportCTFRangeInfo <- function(seasonder_cs_object, file, tableStart = "") {
  res <- seasonder_exportCTFRangeInfo_string(seasonder_cs_object, tableStart)
  writeLines(res$out_str, con = file)
  invisible(res$range_info)
}