#---------------------------------------------------------------------------#
#           Group on Data Assimilation Development -GDAD/CPTEC/INPE         #
#---------------------------------------------------------------------------#
#BOP
#
# !ISCRIPT: total_count.gnu
#
# !DESCRIPTON: Este script é utilizado para gerar o scatterplot entre
#              observacao e simulações
#
#              Para a produção destes graficos são utilizados os
#              seguintes arquivos [ no formato ascii ]:
#
#              Cada arquivo possui a informação correspondente as
#              variaveis convencionais processadas pelo GSI, sendo a
#              primeira coluna referente as observações, a segunda
#              referente ao modelo(GSI). Os tipos de observação
#              processados são os seguintes:
#
#                 * SST : Temperatura da Superfície do Mar
#                 * PS  : Pressão à Superfície
#                 * PW  : Água Precipitável
#                 * UV  : Vento
#                 * T   : Temperatura do Ar
#                 * Q   : Umidade Específica
#
#
#\\
#\\

#
# !INTERFACE:
#

VarPlot=${1}
TimeFile=${2}
DataType=${3} # tipo de informacao: 1 - First Guess, 2 Analise
DataLev=${4}
PathInp=${5}
PathOut=${6}


#
#
# !REVISION HISTORY:
#  17 May 2012 - J. G. de Mattos - Initial Version
#
#
#
# !SEE ALSO:
#   
#   diag_gsi_conf.f90 : programa fortran que pos-processa as
#                       informações de diagnostico do GSI.
#
#EOP
#---------------------------------------------------------------------------//
#BOC

#
# BEGIN: Configuração
#


#
#Nomes Para titulo e Arquivo de Saida
#

NTitle[1]="Temp. da Sup. do Mar"
NTitle[2]="Pressão em Superfície"
NTitle[3]="Água Precipitável"
NTitle[4]="Vento"
NTitle[5]="Temperatura do ar"
NTitle[6]="Umidade Específica"
NTitle[7]="Radiocultação GPS"

varname[1]="sst"
varname[2]="ps"
varname[3]="pw"
varname[4]="uv"
varname[5]="t"
varname[6]="q"
varname[7]="gps"

#
# tipo de informacao
#

type[1]="fgs"
type[2]="anl"

#
# Pegando Ano, Mes, dia e Hora
#

ano=${TimeFile:0:4}
mes=${TimeFile:4:2}
dia=${TimeFile:6:2}
hor=${TimeFile:8:2}

#
# Titulo dos eixos
#

TitleXAxis="Observação"
TitleYAxis="Modelo"

#
# Intervalos pre-definidos
#

minval[1]=200; maxval[1]=330   # sst [K]
minval[2]=500; maxval[2]=1100  # ps  [hPa]
minval[3]=200; maxval[3]=330   # pw  []
minval[4]=-20; maxval[4]=100   # uv  []
minval[5]=150; maxval[5]=330   # t   []
minval[6]=000; maxval[6]=020   # q   []
minval[7]=000; maxval[7]=350   # gps []

#
# Arquivos de entrada/saida
#

FNameInp="${PathInp}/sso_${varname[$VarPlot]}.${TimeFile}.${type[DataType]}.dat"
FNameOut="${PathOut}/scatter_${varname[$VarPlot]}_${DataLev}.${TimeFile}.${type[DataType]}"

#
# Iniciando gnuplot
#


gnuplot << EOF
#cat << EOF >> teste.teste

#
# Definindo titulo
#

set title "${NTitle[$VarPlot]} : ${dia}/${mes}/${ano} ${hor}UTC em ${DataLev} hPa"

#
# Identificando os Eixos
#

set xlabel "${TitleXAxis}"
set ylabel "${TitleYAxis}"


#
# Gerando Estatistica para imprimir no grafico
#

stats "${FNameInp}" using 1:(\$3==${DataLev}?\$2:1/0) prefix "A"
stats "${FNameInp}" using (\$3==${DataLev}?\$2-\$1:1/0) prefix "B"

#
# definindo ranges e formato do grafico
#

set xrange[${minval[$VarPlot]}:${maxval[$VarPlot]}]
set yrange[${minval[$VarPlot]}:${maxval[$VarPlot]}]
set size ratio -1



#
# Plotando informacoes no grafico
#

# labels
set label 1 at graph 0.05 , graph 0.95 gprintf("a = %.3f",A_slope) left front
set label 2 at graph 0.05 , graph 0.90 gprintf("b = %.3f",A_intercept) left front
set label 3 at graph 0.05 , graph 0.85 gprintf("count = %g",B_records) left front
set label 4 at graph 0.05 , graph 0.80 gprintf("bias = %.3f",B_mean) left front
set label 5 at graph 0.05 , graph 0.75 gprintf("stddev = %.3f",B_stddev) left front
set label 6 at graph 0.05 , graph 0.70 gprintf("xycorr = %.3f",A_correlation) left front
# retangulo

set object 1 rect from graph 0.03, graph 0.67 to graph 0.32, graph 0.98 front
set object 1 rect fc rgb "white" fillstyle transparent solid 0.5 front

#
# desenhando linhas dos desvios padroes
#

#set arrow 1 from A_mean_x-A_stddev_x,graph 0 to A_mean_x-A_stddev_x,graph 1 nohead lt 2
#set arrow 2 from A_mean_x+A_stddev_x,graph 0 to A_mean_x+A_stddev_x,graph 1 nohead lt 2

#set arrow 3 from graph 0, second A_mean_x-A_stddev_x  to graph 1,second A_mean_x-A_stddev_x nohead lt 2
#set arrow 4 from graph 0, second A_mean_x+A_stddev_x  to graph 1,second A_mean_x+A_stddev_x nohead lt 2

#
# linha diagonal
#

set style line 1 lt 3 lc rgb "black" lw 1
set style arrow 1 nohead ls 1
set arrow from graph 0, graph 0 to graph 1, graph 1 as 1 front

#
# adicionando background
#

set object 2 rect from graph 0, graph 0 to graph 1, graph 1 behind
set object 2 rect fc rgb "#FFF9D9" fillstyle solid 1.0
#set object 2 rect fc rgb "gray" fillstyle solid 1.0

set style fill pattern 2 bo 1


#
# adicionando grade
#
#EPV
#set grid lc rgb "black" front

#
# Arquivo de saida
#

#set termoption dash

set encoding utf8
set term postscript enhanced color
set output "${FNameOut}.eps"

#
# posicionando a legenda
#

set key right bottom

#
# funcao para plotar regiao de desvio padrao
#

f(x) = A_slope*x + A_intercept
#up(x) = f(x) + A_stddev_x
#do(x) = f(x) - A_stddev_x

#
# estilo da regiao demarcada com o desvio padrao
#

set style line 2 lc rgb "red"
set style fill solid

#
# estilo da linha de ajuste
# 

set style line 1 lt 1 lc rgb "blue" lw 3

#
# plotando 
#

#plot '+' u 1:(up(\$1)):(do(\$1)) w filledcurves lt 2 lc rgb "#EDECEB" title 'Stddev Obs',\
#"${FNameInp}" u 1:(\$3==${DataLev}?\$2:1/0) notitle  w p pt 7 lt 1 ps 1, \
#(f(x)) ls 1 title 'Ajuste Linear',\
#(up(x)) ls 2 notitle ,\
#(do(x)) ls 2 notitle

plot "${FNameInp}" u 1:(\$3==${DataLev}?\$2:1/0) notitle  w p pt 7 lt 1 ps 1, \
(f(x)) ls 1 title 'Ajuste Linear'

#
# saindo
#

exit

EOF

convert -rotate 90 -density 100 ${FNameOut}.eps -flatten ${FNameOut}.png
#convert -rotate 90 -density 100 ${FNameOut}.eps -flatten tmp.png
#/bin/bash ../bin/autotrim tmp.png ${FNameOut}.png
rm -fr ${FNameOut}.eps