#!/bin/bash
######################################################################
### Script to modify the orbital parameters of a file "startfi.nc" ###
### according to the date set in the file "run_PEM.def"            ###
######################################################################

# Name of the file to be modified
name_file="startfi.nc"

# Name of the file containing the orbital data
orb_data="obl_ecc_lsp.asc"


######################################################################
# Check if files necessary for the script exist
if [ ! -f "$orb_data" ]; then
    echo "Error: file \"$orb_data\" not found!"
    exit 1
fi
if [ ! -f "startfi.nc" ]; then
    echo "Error: file \"$name_file\" not found!"
    exit 1
fi
if [ ! -f "run_PEM.def" ]; then
    echo "Error: file \"run_PEM.def\" not found!"
    exit 1
fi

# Get the number of Earth years before present to start the PEM run from "run_PEM.def"
eyears_bp_ini=$(grep 'year_earth_bp_ini=' run_PEM.def | cut -d '=' -f 2 | tr -d '[:space:]')
echo "The starting date is $eyears_bp_ini Earth years according to \"run_PEM.def\"."

# Get the new values for the orbital parameters
yearlask=$(echo "$eyears_bp_ini" | bc -l)
found=false
read -r y1 obl1 ecc1 lsp1 < "$orb_data"
while read -r y2 obl2 ecc2 lsp2; do
    if [ "$first_line" = true ]; then
        first_line=false
        continue # Skip the first line as it's already read outside the loop
    fi
    if [ "$(echo "$y1 >= $yearlask && $yearlask > $y2" | bc)" -eq 1 ]; then
        found=true
        # Calculate linear interpolations for each orbital parameter
        new_obl=$(echo "($obl2 - $obl1)/($y2 - $y1)*($yearlask - $y1) + $obl1" | bc -l)
        new_ecc=$(echo "($ecc2 - $ecc1)/($y2 - $y1)*($yearlask - $y1) + $ecc1" | bc -l)
        if [ "$(echo "$lsp2 - $lsp1 > 300. || $lsp2 - $lsp1 < -300." | bc)" -eq 1 ]; then # If modulo is "crossed" through the interval
            if [ "$(echo "$lsp1 < $lsp2" | bc -l)" -eq 1 ]; then # Lsp should be decreasing
                lsp1=$(echo "$lsp1 + 360" | bc -l)
            else # Lsp should be increasing
                lsp2=$(echo "$lsp2 + 360" | bc -l)
            fi
        fi
        new_Lsp=$(echo "($lsp2 - $lsp1)/($y2 - $y1)*($yearlask - $y1) + $lsp1" | bc -l)
        echo "New obliquity    = $new_obl"
        echo "New eccentricity = $new_ecc"
        echo "New Lsp          = $new_Lsp"
        echo "Done!"
        break
    fi
    y1=$y2
    obl1=$obl2
    ecc1=$ecc2
    lsp1=$lsp2
done < "$orb_data"
if [ "$found" = false ]; then
    echo "Error: the current year could not be found in the file \"$orb_data\"!"
    exit 1
fi


# Compute the modified values in the file
if [ ! -f "$name_file" ]; then
    echo "Error: file \"$name_file\" not found!"
    exit 1
fi
year_day=669.
halfaxe=227.94
pi=$(echo "4.*a(1.)" | bc -l)
degrad=$(echo "180./$pi" | bc -l)

periheli=$(echo "$halfaxe*(1. - $new_ecc)" | bc -l)
aphelie=$(echo "$halfaxe*(1. + $new_ecc)" | bc -l)

tan=$(echo "s(0.5*$new_Lsp/$degrad)/c(0.5*$new_Lsp/$degrad)" | bc -l)
zx0=$(echo "-2.*a($tan*sqrt((1. - $new_ecc)/(1. + $new_ecc)))" | bc -l)
if [ $(echo "$zx0 <= 0." | bc -l) -eq 1 ]; then
    zx0=$(echo "$zx0 + 2.*$pi" | bc -l)
fi
peri_day=$(echo "$year_day*(1. - ($zx0 - $new_ecc*s($zx0))/(2.*$pi))" | bc -l)

ztheta=$(echo "($new_iniLs - $new_Lsp)/$degrad" | bc -l)
tan=$(echo "s(0.5*$ztheta)/c(0.5*$ztheta)" | bc -l)
zx0=$(echo "2.*a($tan*sqrt((1. - $new_ecc)/(1. + $new_ecc)))" | bc -l)
xref=$(echo "$zx0 - $new_ecc*s($zx0)" | bc -l)
new_inisol=$(echo "$peri_day + $xref*$year_day/(2.*$pi)" | bc -l)
if [ $(echo "$new_inisol < 0." | bc -l) -eq 1 ]; then
    new_inisol=$(echo "$new_inisol + $year_day" | bc -l)
fi
if [ $(echo "$new_inisol >= $year_day" | bc -l) -eq 1 ]; then
    new_inisol=$(echo "$new_inisol - $year_day" | bc -l)
fi

# Update the netCDF file
# controle(15) = periheli ! min. Sun-Mars distance (Mkm)  ̃206.66
# controle(16) = aphelie  ! max. SUn-Mars distance (Mkm)  ̃249.22
# controle(17) = peri_day ! date of perihelion (sols since N. spring)
# controle(18) = obliquit ! Obliquity of the planet (deg)  ̃23.98
# controle(3)  = day_ini + int(time) ! initial day
ncap2 -O -s "controle(17)=$new_obl" \
         -s "controle(14)=$periheli" \
         -s "controle(15)=$aphelie" \
         -s "controle(16)=$peri_day" \
         -s "controle(2)=$new_inisol" \
         $name_file $name_file

echo "In \"$name_file\":"
echo "New obliquit     = $new_obl"
echo "New eccentricity = $new_ecc -> new periheli = $periheli"
echo "                              -> new aphelie  = $aphelie"
echo "New Lsp          = $new_Lsp -> new peri_day = $peri_day"
echo "New initial Ls   = $new_iniLs -> New initial sol = $new_inisol"
echo "Done!"