source: trunk/LMDZ.COMMON/libf/evolution/deftank/modify_startfi_orbit.sh

#!/bin/bash
######################################################################
### Script to modify the orbital parameters of a file "startfi.nc" ###
######################################################################

# Name of the file
name_file="startfi.nc"

# New values for the orbital parameters
new_obl=25.18941689
new_ecc=0.09340902
new_Lsp=251.06881714
new_iniLs=0.


######################################################################
# Calculate modified values
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!"