| 1 | #!/bin/bash |
|---|
| 2 | ###################################################################### |
|---|
| 3 | ### Script to modify the orbital parameters of a file "startfi.nc" ### |
|---|
| 4 | ###################################################################### |
|---|
| 5 | set -e |
|---|
| 6 | trap 'echo -e "\033[31mError: an issue occurred in the script on line $LINENO! Please review the command and try again.\033[0m"' ERR |
|---|
| 7 | |
|---|
| 8 | ###################################################################### |
|---|
| 9 | # Modify here the parameters for the script |
|---|
| 10 | ########################################### |
|---|
| 11 | # Name of the file |
|---|
| 12 | name_file="startfi.nc" |
|---|
| 13 | |
|---|
| 14 | # New values for the orbital parameters |
|---|
| 15 | new_obl=25.18941689 |
|---|
| 16 | new_ecc=0.09340902 |
|---|
| 17 | new_Lsp=251.06881714 |
|---|
| 18 | new_iniLs=0. |
|---|
| 19 | ###################################################################### |
|---|
| 20 | |
|---|
| 21 | |
|---|
| 22 | # Calculate modified values |
|---|
| 23 | if [ ! -f "$name_file" ]; then |
|---|
| 24 | echo "Error: file \"$name_file\" not found!" |
|---|
| 25 | exit 1 |
|---|
| 26 | fi |
|---|
| 27 | year_day=669. |
|---|
| 28 | halfaxe=227.94 |
|---|
| 29 | pi=$(echo "4.*a(1.)" | bc -l) |
|---|
| 30 | degrad=$(echo "180./$pi" | bc -l) |
|---|
| 31 | |
|---|
| 32 | periheli=$(echo "$halfaxe*(1. - $new_ecc)" | bc -l) |
|---|
| 33 | aphelie=$(echo "$halfaxe*(1. + $new_ecc)" | bc -l) |
|---|
| 34 | |
|---|
| 35 | tan=$(echo "s(0.5*$new_Lsp/$degrad)/c(0.5*$new_Lsp/$degrad)" | bc -l) |
|---|
| 36 | zx0=$(echo "-2.*a($tan*sqrt((1. - $new_ecc)/(1. + $new_ecc)))" | bc -l) |
|---|
| 37 | if [ $(echo "$zx0 <= 0." | bc -l) -eq 1 ]; then |
|---|
| 38 | zx0=$(echo "$zx0 + 2.*$pi" | bc -l) |
|---|
| 39 | fi |
|---|
| 40 | peri_day=$(echo "$year_day*(1. - ($zx0 - $new_ecc*s($zx0))/(2.*$pi))" | bc -l) |
|---|
| 41 | |
|---|
| 42 | ztheta=$(echo "($new_iniLs - $new_Lsp)/$degrad" | bc -l) |
|---|
| 43 | tan=$(echo "s(0.5*$ztheta)/c(0.5*$ztheta)" | bc -l) |
|---|
| 44 | zx0=$(echo "2.*a($tan*sqrt((1. - $new_ecc)/(1. + $new_ecc)))" | bc -l) |
|---|
| 45 | xref=$(echo "$zx0 - $new_ecc*s($zx0)" | bc -l) |
|---|
| 46 | new_inisol=$(echo "$peri_day + $xref*$year_day/(2.*$pi)" | bc -l) |
|---|
| 47 | if [ $(echo "$new_inisol < 0." | bc -l) -eq 1 ]; then |
|---|
| 48 | new_inisol=$(echo "$new_inisol + $year_day" | bc -l) |
|---|
| 49 | fi |
|---|
| 50 | if [ $(echo "$new_inisol >= $year_day" | bc -l) -eq 1 ]; then |
|---|
| 51 | new_inisol=$(echo "$new_inisol - $year_day" | bc -l) |
|---|
| 52 | fi |
|---|
| 53 | |
|---|
| 54 | # Update the netCDF file |
|---|
| 55 | # controle(15) = periheli ! min. Sun-Mars distance (Mkm) ̃206.66 |
|---|
| 56 | # controle(16) = aphelie ! max. Sun-Mars distance (Mkm) ̃249.22 |
|---|
| 57 | # controle(17) = peri_day ! date of perihelion (sols since N. spring) |
|---|
| 58 | # controle(18) = obliquit ! Obliquity of the planet (deg) ̃23.98 |
|---|
| 59 | # controle(3) = day_ini + int(time) ! initial day |
|---|
| 60 | ncap2 -O -s "controle(17)=$new_obl" \ |
|---|
| 61 | -s "controle(14)=$periheli" \ |
|---|
| 62 | -s "controle(15)=$aphelie" \ |
|---|
| 63 | -s "controle(16)=$peri_day" \ |
|---|
| 64 | -s "controle(2)=$new_inisol" \ |
|---|
| 65 | $name_file $name_file |
|---|
| 66 | |
|---|
| 67 | echo "In \"$name_file\":" |
|---|
| 68 | echo "New obliquit = $new_obl" |
|---|
| 69 | echo "New eccentricity = $new_ecc -> new periheli = $periheli" |
|---|
| 70 | echo " -> new aphelie = $aphelie" |
|---|
| 71 | echo "New Lsp = $new_Lsp -> new peri_day = $peri_day" |
|---|
| 72 | echo "New initial Ls = $new_iniLs -> New initial sol = $new_inisol" |
|---|
| 73 | echo "Done!" |
|---|
