@pwessel
One last bit … sorry it’s for my sanity … Is there an easy way to get rid of the first folded “September” ?
I tried to play with -Ba1O+phase but it was not conclusive…
#!/bin/bash
# basemap, begin, coast, grdimage, grdmath, info, inset, legend, makecpt, math, plot, text
gmt math -foT --TIME_UNIT=d -T0/365/1 2 PI MUL T 1 SUB MUL 365 DIV STO@r POP \
0.006818 \
@r COS 0.39912 MUL SUB \
@r SIN 0.0757 MUL ADD \
2 @r MUL COS 0.006758 MUL SUB \
2 @r MUL SIN 0.000907 MUL ADD \
3 @r MUL COS 0.002697 MUL SUB \
3 @r MUL SIN 0.00148 MUL ADD \
180 PI DIV MUL = declination.txt
active_day="2012-12-21T"
active_day_rel=$(gmt math -Q $active_day $(gmt math -Q -fT $active_day --FORMAT_DATE_OUT=yyyy --TIME_UNIT=y --FORMAT_CLOCK_OUT=- =)-01-01T SUB --TIME_UNIT=d =)
active_line=$(gmt math -Q $active_day_rel 1 ADD =)
wrap_ref="1970-09-23T"
wrap_shift=$(gmt math -Q $wrap_ref 1970-12-31T SUB --TIME_UNIT=d 365 DIV =)
wrap_close=$(gmt math -Q $wrap_shift 1 ADD =)
wrap_day=$(gmt math -Q $active_day_rel 365 DIV STO@r $wrap_close GT 1 MUL STO@r2 POP @r @r2 SUB =)
delta=$(echo $(awk -v var=$active_line 'NR==var {print $2}' declination.txt))
delta_max=$(gmt info declination.txt -C -o3)
gmt begin test png
gmt plot declination.txt -R$wrap_shift/$wrap_close/-25/25 -wy -JX10c -BW -Byaf+u@. --MAP_FRAME_TYPE=graph
echo $wrap_day $delta | gmt plot -Sc0.1c -Gred -Wfaint,black
gmt basemap -R$wrap_ref/$(gmt math -Q -fT $wrap_ref 365 ADD --TIME_UNIT=d =)/-25/25 -BS -Bxa1Of --FORMAT_DATE_MAP=o --FORMAT_TIME_PRIMARY_MAP=a --MAP_FRAME_TYPE=graph -Y5c
gmt text -Dj0.2c -N -F+f+j -Y-5c <<- EOF
1971-06-21T 0 12p,Times-Roman BC Boreal summer
1970-12-31T 0 12p,Times-Roman BC Austral summer
$wrap_ref 30 12p,Symbol RT d
EOF
gmt inset begin -DjBR+w3.5c+o-3c/0.25c
gmt grdmath -Rg -I0.5 X SIND = illu.nc
gmt makecpt -Cgray -T-1/0 -I > illu.cpt
gmt coast -Rg -JG0/0/3c+w$(gmt math -Q $delta_max -1 MUL =) -Dc -Ggray -A5000 -Bg -Wfaint -X0.25 -Y0.25
gmt grdimage -JG0/0/3c+w$(gmt math -Q $delta $delta_max ADD -1 MUL =) illu.nc -C -t25
gmt plot -JG0/0/3.5c+w$(gmt math -Q $delta_max -1 MUL =) -Wthin,black -Am -X-0.25 -Y-0.25 <<- EOF
0 -89
0 89
EOF
gmt plot -Sv0.15c+ea+h0+a60+gred -Wthicker,red -N <<- EOF
-90 $delta $(gmt math -Q $delta $delta_max ADD -1 MUL =) 0.3c
EOF
gmt inset end
gmt legend -DjTR+o-3c/-1c -F <<- EOF
H 14p,Helvetica-Bold Solar declination
D 5p 1p
S 0.5c c 0.1c red faint 1c Observed day
S 0.5c v0.15c+h0+ea+a60 0.3c red thicker,red 1c Solar equatorial plan
EOF
gmt end show