#!/bin/bash
### dclrinexplot.sh
### run fortran program dclrinex and produces plots with GMT
### version 20210720

#set -x
### functions

_usage()
{
    cat << END

 Usage:   dclrinexplot.sh <sss1> <sss2> <yydoy[.xxxx]> <dur> "[<OPTION>]"

 where:

        sss1 = 4- or 9-character name for station 1 (travel)
        sss2 = 4- or 9-character name for station 2 (local)
       yydoy = first day of Rinex
       .xxxx = option to specify precise start time
         dur = Number of days to process
      OPTION = P1ISC1: Use C1 measurements instead of P1, for a station with no P1 data
               LZYYYY: Use files LZYYYYmj.day to report station 2 to the same reference as station 1
               GPSSOL to process only GPS (default = all four GNSS)
               BDSGEO to process BeiDou geostationary satellites (def=no)
               NOBDS2 to remove BeiDou 2 satellites (def=process)
               NOBDS3 to remove BeiDou 3 satellites (def=process)
               NOBDS to remove BeiDou satellites (def=process)
               NOGAL to remove Galileo satellites (def=process)
               NOGLO to remove GLONASS satellites (def=process)
               ALLCOD to process all codes including non CGGTTS V2E (def=only CGGTTS codes)
               BIGDIF to write all raw code differences in the output file .dif (def=smaller file)
               ELMIxx to remove elev<xx deg (def=05)
               If several options are used, they should be entered between " ".

 Output file = sss1sss2yydoy_dur.ext
 Warning: Remove files which have the same filename as this output file!

END
}

_tmp_cleanup()
{
    rm -f ./*[CcPpEeBbRrDd][1235678c].tmp
    rm -f ./*[CcPpEeBbRrDd][1235678c].err
}

_get_lin_axes() # $FNAME
{
    FNAME=${1}

    ymax=$(awk -v tmax=-99999999 '{ if ($2 > tmax) tmax=$2 } END { print tmax } ' "$FNAME")
    ymin=$(awk -v tmin=99999999  '{ if ($2 < tmin) tmin=$2 } END { print tmin } ' "$FNAME")
    ymax=$(echo "$ymax" | awk '{ print int($1)+1 }')
    ymin=$(echo "$ymin" | awk '{ print int($1)-1 }')

    yspan=$(echo "$ymax $ymin" | awk '{ print $1-$2 }')

    yint=$(echo "$yspan" | awk '{ print $1/30 }')

    yatic=100
    if [[ ${yspan} -lt 200 ]]; then
        yatic=10
    fi
    if [[ ${yspan} -lt 20 ]]; then
        yatic=1
    fi
    yftic=$(echo ${yatic} | awk '{ print $1/5 }')

    xmax=$(awk -v tmax=-99999999 ' { if ($1 > tmax) tmax=$1 } END { print tmax } ' "$FNAME")
    xmin=$(awk -v tmin=99999999  ' { if ($1 < tmin) tmin=$1 } END { print tmin } ' "$FNAME")
    xmax=$(echo "${xmax}" | awk '{ print int($1)+1 }')
    xmin=$(echo "${xmin}" | awk '{ print int($1) }')

    xspan=$(echo "${xmax} ${xmin}" | awk '{ print $1-$2 }')
    xftic=$(echo "${xspan}" | awk '{ print int($1/10)+1 }')
    xatic=$(echo "${xftic}" | awk '{ print $1*2 }')
    xint=$(echo "${xspan}" | awk '{ print $1/30 }')
}

_get_log_axes() # $FNAME
{
    FNAME=${1}

    xmax=$(awk -v tmax=0 '{ if ($1 > tmax) tmax=$1 } END { print tmax } ' "${FNAME}")
    xmin=$(awk -v tmin=99999999 ' { if ($1 < tmin) tmin=$1 } END { print tmin } ' "${FNAME}")
    ymax=$(awk -v tmax=0 ' { if ($2 > tmax) tmax=$2 } END { print tmax } ' "${FNAME}")
    ymin=$(awk -v tmin=99999999 ' { if ($2 < tmin) tmin=$2 } END { print tmin } ' "${FNAME}")

    xmax=$(echo "${xmax}"  | awk '{ print 10^(int(log($1)/2.302585093)+1) }' )
    xmin=$(echo "${xmin}"  | awk '{ print 10^(int(log($1)/2.302585093)) }' )
    ymax=$(echo "${ymax}"  | awk '{ print 10^(int(log($1)/2.302585093)) }' )
    ymin=$(echo "${ymin}"  | awk '{ print 10^(int(log($1)/2.302585093)-1) }' )

    xspan=$(echo "${xmax} ${xmin}"  | awk '{ print $1/$2 }' )
    xftic=$(echo "${xmax} ${xspan}" | awk '{ print $1/$2 }' )
    xatic=$(echo "${xftic}"         | awk '{ print $1/1  }' )

    yspan=$(echo "${ymax} ${ymin}"  | awk '{ print $1/$2 }' )
    yftic=$(echo "${ymax} ${yspan}" | awk '{ print $1/$2 }' )
    yatic=$(echo "${yftic}"         | awk '{ print $1/1  }' )
}

_get_med_text() # $FNAME $SIZE $POS $DEVS_LIS
{
    FNAME=${1}
    SIZE=${2}
    POS=${3}
    DEVS_LIS=${4}

    while read -r x
    do
        med=$(echo "${x}" | awk '{ print $3 }')
        Med=$(echo "${med}" | awk '{ printf "%7.2f", $1 }')
        xpos=$(echo "${xmin} ${xint}" | awk '{ print $1+2*$2 }')
        if [[ "x${POS}" = "xBOT" ]]; then
            ypos=$(echo "${ymin} ${ymax} ${yint}" | awk '{ print 1.*$1+0*$2+1*$3}')
        elif [[ "x${POS}" = "xTOP" ]]; then
            ypos=$(echo "${ymin} ${ymax} ${yint}" | awk '{ print 0.*$1+1*$2-1*$3}')
        fi
        if [[ -z ${cmd} ]]; then
            echo "${xpos} ${ypos} ${SIZE} 0 1 BL Median = ${Med} ns" >> "${DEVS_LIS}"
        else
            echo "${xpos} ${ypos} ${SIZE},1,blue BL Median = ${Med} ns" >> "${DEVS_LIS}"
        fi
    done <"${FNAME}"
}

_get_tdev_text() # $FNAME $SIZE $CODE $POS $DEVS_LIS
{
    FNAME=${1}
    SIZE=${2}
    CODE=${3}
    POS=${4}
    DEVS_LIS=${5}

    i=0
    while read -r tau tdev
    do
        Tau=$(echo "${tau}"   | awk '{ printf "%6.0f", $1 }')
        Tdev=$(echo "${tdev}" | awk '{ printf "%6.0f", $1*1.e12 }')
        if [[ "x$POS" = "xLEFT" ]]; then
            xpos=$(echo "${xmax} ${xint}" | awk '{ print $1+$2*5 }' )
        elif [[ "x$POS" = "xRIGHT" ]]; then
            xpos=$(echo "${xmax} ${xint}" | awk '{ print $1+$2*12 }' )
        fi
        ypos=$(echo "${ymax} ${yint} ${i}" | awk '{ printf "%9.2f", $1+$2*($3-4) }' )
        if [[ -z ${cmd} ]]; then
            echo "${xpos} ${ypos} ${SIZE} 0 1 LM ${Tau} s: ${CODE}= ${Tdev} ps" >> "${DEVS_LIS}"
        else
            echo "${xpos} ${ypos} ${SIZE},1,blue LM ${Tau} s: ${CODE}= ${Tdev} ps" >> "${DEVS_LIS}"
        fi
        i=$((i + 1))
    done <"${FNAME}"
}

_gmt4_clear_all()
{
    rm .gmt* 2> /dev/null
    rm gmt.conf 2> /dev/null
    rm gmt.history 2> /dev/null
}

_gmt4_set()
{
    gmtset N_HEADER_RECS 0
    gmtset BASEMAP_AXES WeSn
    gmtset LABEL_FONT_SIZE 16p
    gmtset HEADER_FONT_SIZE 20p
    gmtset PAPER_MEDIA a4
}

_gmt5_clear_all()
{
    ${cmd} clear history 2>/dev/null
    ${cmd} clear conf 2>/dev/null
    ${cmd} clear cache 2>/dev/null
}

_gmt5_set()
{
    ${cmd} gmtset IO_HEADER 0
    ${cmd} gmtset MAP_FRAME_AXES WeSn
    ${cmd} gmtset FONT_LABEL 16p
    ${cmd} gmtset FONT_TITLE 20p
    ${cmd} gmtset PS_MEDIA a4
}

_plot_code() # $COD1 $COD2
{
    COD1=${1}
    COD2=${2}

    ### diff plot
    if [[ -n "${COD2}" ]]; then bas=""; else bas="MJD"; fi

    # COD1 diff plot
    _get_lin_axes "${COD1}.tmp"
    ${cmd} psxy "${COD1}.tmp" -S-0.04 -Xa1.00i -Ya4.00i -R"${xmin}/${xmax}/${ymin}/${ymax}" -JX6.10i/2.50i -B"f${xftic}a${xatic}:${bas}:/f${yftic}a${yatic}:${COD1} /ns::.${TITLE}:" -W1,0/0/0 -K > "${OUTFILE}.ps"

    # COD2 diff plot
    if [[ -n "${COD2}" ]]; then
        _get_lin_axes "${COD2}.tmp"
        ${cmd} psxy -O "${COD2}.tmp" -S-0.04 -Xa1.00i -Ya1.00i -R"${xmin}/${xmax}/${ymin}/${ymax}" -JX6.10i/2.50i -B"f${xftic}a${xatic}:MJD:/f${yftic}a${yatic}:${COD2} /ns:" -W1,0/0/0 -K >> "${OUTFILE}.ps"
    fi


    ### Add text
    rm -f devs.lis
    # med
    _get_med_text "med${COD1,,}.tmp" 12 TOP devs.lis
    if [[ -n "${COD2}" ]]; then
        _get_med_text "med${COD2,,}.tmp" 12 BOT devs.lis
    fi
    # tdev
    _get_tdev_text "tdev${COD1,,}.tmp" 9 "${COD1}" LEFT  devs.lis
    if [[ -n "${COD2}" ]]; then
        _get_tdev_text "tdev${COD2,,}.tmp" 9 "${COD2}" RIGHT devs.lis
    fi
    #
    if [[ -n "${cmd}" ]]; then opt="-F+f+j"; else opt="-G0/0/255"; fi
    ${cmd} pstext -N -O -K -R"${xmin}/${xmax}/${ymin}/${ymax}" -Xa1.00i -Ya1.00i -JX6.20i/5.50i ${opt} devs.lis >> "${OUTFILE}.ps"
    rm -f devs.lis

    ### TDEV plot
    if [[ -n "${COD2}" ]]; then bas=""; kon="-K"; else bas="Tau /s"; kon=""; fi

    # COD1 TDEV plot
    _get_log_axes "tdev${COD1,,}.tmp" 2> "tdev${COD1,,}.err"
    ${cmd} psxy "tdev${COD1,,}.tmp" -O ${kon} -Xa8.50i -Ya3.70i -Sc0.05i -R"${xmin}/${xmax}/${ymin}/${ymax}" -JX2.50il/1.75il -B"f${xftic}a${xatic}:${bas}:/f1a1:Tdev${COD1} /s:" -W5,0/0/0 >> "${OUTFILE}.ps" 2>> "tdev${COD1,,}.err"

    if [[ -s "tdev${COD1,,}.err" ]]; then echo "warning tdev${COD1,,}.tmp"; rm -f "tdev${COD1,,}.err"; fi

    # COD2 TDEV plot
    if [[ -n "${COD2}" ]]; then
        _get_log_axes "tdev${COD2,,}.tmp" 2> "tdev${COD2,,}.err"
        ${cmd} psxy "tdev${COD2,,}.tmp" -O -Xa8.50i -Ya1.00i -Sc0.05i -R"${xmin}/${xmax}/${ymin}/${ymax}" -JX2.50il/1.75il -B"f${xftic}a${xatic}:Tau /s:/f1a1:Tdev${COD2} /s:" -W5,0/0/0 >> "${OUTFILE}.ps" 2>> "tdev${COD2,,}.err"

        if [[ -s "tdev${COD2,,}.err" ]]; then echo "warning tdev${COD2,,}.tmp"; rm -f "tdev${COD2,,}.err"; fi
    fi

    # transform to pdf
    ps2pdf "${OUTFILE}.ps"
    rm -f "${OUTFILE}.ps"
    mv "${OUTFILE}.pdf" "${OUTFILE}_${COD1}${COD2:1:1}.pdf"
}


### arguments
if [[ $# -lt 4 || $# -ge 6 ]]; then
    _usage
    exit 1
fi

### variables
doyd=$(echo "${3}" | cut -c1-5)
OUTFILE="${1:0:4}${2:0:4}${doyd}_${4}"
DATERUN=$(date '+%F')
TITLE="${DATERUN} ${OUTFILE}"

# check GMT software version 4/5
cmd=$(command -v gmt)

if [[ -z ${cmd} ]]; then
    _gmt4_clear_all
    _gmt4_set
else
    _gmt5_clear_all
    _gmt5_set
fi


### core

# Delete old *.tmp and *.err files
_tmp_cleanup


# Invoke dclrinex
dclrinex "${1}" "${2}" "${3}" "${4}" "${5}"


# GPS plot P1/P2
if [[ -s "P1.tmp" && -s "P2.tmp" ]]; then
    _plot_code P1 P2
elif [[ -s "P1.tmp" ]]; then
    _plot_code P1
elif [[ -s "P2.tmp" ]]; then
    _plot_code P2
fi

# GPS plot C1/C2
if [[ -s "C1.tmp" && -s "C2.tmp" ]]; then
    _plot_code C1 C2
elif [[ -s "C1.tmp" ]]; then
    _plot_code C1
elif [[ -s "C2.tmp" ]]; then
    _plot_code C2
fi

# GPS plot C5
if [[ -s "C5.tmp" ]]; then
    _plot_code C5
fi

# GLO plot D1/D2
if [[ -s "D1.tmp" && -s "D2.tmp" ]]; then
    _plot_code D1 D2
elif [[ -s "D1.tmp" ]]; then
    _plot_code D1
elif [[ -s "D2.tmp" ]]; then
    _plot_code D2
fi

# GLO plot R1/R2
if [[ -s "R1.tmp" && -s "R2.tmp" ]]; then
    _plot_code R1 R2
elif [[ -s "R1.tmp" ]]; then
    _plot_code R1
elif [[ -s "R2.tmp" ]]; then
    _plot_code R2
fi

# GAL plot E1/E5
if [[ -s "E1.tmp" && -s "E5.tmp" ]]; then
    _plot_code E1 E5
elif [[ -s "E1.tmp" ]]; then
    _plot_code E1
elif [[ -s "E5.tmp" ]]; then
    _plot_code E5
fi

# GAL plot E6/E7
if [[ -s "E6.tmp" && -s "E7.tmp" ]]; then
    _plot_code E6 E7
elif [[ -s "E6.tmp" ]]; then
    _plot_code E6
elif [[ -s "E7.tmp" ]]; then
    _plot_code E7
fi

# GAL plot E8
if [[ -s "E8.tmp" ]]; then
    _plot_code E8
fi

# BDS plot B1/B2
if [[ -s "B1.tmp" && -s "B2.tmp" ]]; then
    _plot_code B1 B2
elif [[ -s "B1.tmp" ]]; then
    _plot_code B1
elif [[ -s "B2.tmp" ]]; then
    _plot_code B2
fi

# BDS plot BC/B5
if [[ -s "BC.tmp" && -s "B5.tmp" ]]; then
    _plot_code BC B5
elif [[ -s "BC.tmp" ]]; then
    _plot_code BC
elif [[ -s "B5.tmp" ]]; then
    _plot_code B5
fi

# BDS plot B3/B7
if [[ -s "B3.tmp" && -s "B7.tmp" ]]; then
    _plot_code B3 B7
elif [[ -s "B3.tmp" ]]; then
    _plot_code B3
elif [[ -s "B7.tmp" ]]; then
    _plot_code B7
fi

# BDS plot B8
if [[ -s "B8.tmp" ]]; then
    _plot_code B8
fi



echo "Output files: ${OUTFILE}"

# clear GMT files
if [[ -z ${cmd} ]]; then
    _gmt4_clear_all
else
    _gmt5_clear_all
fi

### EOS
exit 0