module HistoRoutines
implicit none


private
public Histo
public FillHisto
public ConstructHisto
public DestructHisto
public PrintHisto
public DrawHisto

! definizione oggetto istogramma
type Histo
    integer                             :: nBin
    real                                :: xmin, xmax
    integer                             :: underflow, overflow
    integer, dimension(:), allocatable  :: counts
    real, dimension(:), allocatable     :: errors
    real, dimension(:), allocatable     :: BinCenters

    contains

    procedure           :: FillHisto
    procedure           :: ConstructHisto
    procedure           :: DestructHisto
    procedure           :: PrintHisto

end type

contains

! costruisce l'istogramma dando il numero di bin
! nBin e l'intervallo [xmin, xmax]
subroutine ConstructHisto(thisHisto,nBin,xmin,xmax)
implicit none
class(Histo), intent(inout) :: thisHisto
integer, intent(in)     :: nBin
real, intent(in) :: xmin, xmax
real :: dx, xx
integer :: iNumb

    thisHisto%nBin = nBin

    allocate(thisHisto%counts(nBin))
    allocate(thisHisto%errors(nBin))
    allocate(thisHisto%BinCenters(nBin))

    thisHisto%counts = 0
    thisHisto%xmin = xmin
    thisHisto%xmax = xmax

    dx = ( xmax - xmin ) / nBin
    do iNumb = 1, nBin
        xx = xmin + (iNumb - 0.5) * dx
        thisHisto%BinCenters(iNumb) = xx
    enddo

end subroutine ConstructHisto

! distruggere l'istogramma
subroutine DestructHisto(thisHisto)
class(Histo), intent(inout) :: thisHisto

    if( allocated(thisHisto%counts) ) then
        deallocate(thisHisto%counts)
    endif
    if( allocated(thisHisto%errors) ) then
        deallocate(thisHisto%errors)
    endif
    if( allocated(thisHisto%BinCenters) ) then
        deallocate(thisHisto%BinCenters)
    endif


end subroutine

! riempie l'istogramma, da chiamare in un ciclo do
subroutine fillHisto(thisHisto,xx)
implicit none
class(Histo), intent(inout) :: thisHisto
real, intent(in) 	    :: xx
integer                 :: iNumb, iBin
real                    :: dx

!integer, dimension(thisHisto%nBin) :: counts
real :: xmin, xmax
integer :: nBin
real :: underflow, overflow

xmin = thisHisto%xmin
xmax = thisHisto%xmax
nBin = thisHisto%nBin

    ! ampiezza del bin
	dx = ( xmax - xmin ) / nBin

    ! indice del bin dove cade xx
	iBin = int( ( xx - xmin ) / dx ) + 1

    ! riempimento dei conteggi
    if( iBin > 0 .and. iBin < Nbin + 1 ) then
        thisHisto%counts(iBin) = thisHisto%counts(iBin) + 1
        thisHisto%errors(iBin) = sqrt(1.0*thisHisto%counts(iBin))
    elseif( iBin < 1 ) then
        underflow = underflow + 1
    elseif( iBin > nBin ) then
        overflow = overflow + 1
    endif

end subroutine fillHisto


! printHisto deve essere chiamata alla fine del ciclo do
! che riempie l'istogramma.
! la routine chiede in input l'istogramma thisHisto
! e anche il nome del file 'fileName' dove scrivere
! i conteggi e i rispettivi errori
! nell'unità input/output specificato da 'fileUnit'.
! Nota: si assumono errori poissoniani.
subroutine printHisto(thisHisto,fileName,fileUnit)
implicit none
class(Histo), intent(inout) :: thisHisto
integer, intent(in)                     :: fileUnit
character(len = 128), intent(in)        :: fileName
integer                                 :: iNumb, ios
character(len = 128)                    :: fileNameTemp

fileNameTemp=trim(fileName)//".dat"

    open(unit=fileUnit,file=fileNameTemp,status='replace',iostat=ios)
    if (ios/=0) then
        print*,"Failed to open ",filename,"."
    end if

    
    do iNumb = 1, thisHisto%nBin
        write( unit = fileUnit , fmt = * )  thisHisto%BinCenters(iNumb),thisHisto%counts(iNumb),&
                                            thisHisto%errors(iNumb)
    enddo

    close(fileUnit)

end subroutine printHisto


! drawHisto disegna in formato istogramma i conteggi contenuti
! in fileHisto generando loadFile (un file .txt) nell'unità
! loadFileUnit. Il file loadFile contiene le istruzioni
! gnuplot per prepare il grafico.
! Le variabli xlabel e ylabel sono i titoli delle assi x e y
! mentre invece histoTitle è il titolo dell'istogramma.
! xmin e xmax specificano il range in x mentre ymin e ymax
! quello in y.
! questa routine chiama direttamente gnuplot e in output
! genera l'istogramma salvato nell'immagine loadFile.eps
! NOTA: la routine permette di disegnare un istogramma
!       alla volta.
subroutine drawHisto( fileHisto, loadFile, loadFileUnit,&
                      xlabel, xmin, xmax ,&
                      ylabel, ymin, ymax ,&
                      histoTitle )
implicit none
character(len = 128), intent(in)  :: fileHisto
character(len = 128), intent(out) :: loadFile
character(len = 128), intent(in)  :: xlabel
character(len = 128), intent(in)  :: ylabel
character(len = 128), intent(in)  :: histoTitle
real, intent(in)        :: xmin, xmax, ymin, ymax
integer, intent(in)     :: loadFileUnit
integer                 :: ios
character(len = 128)    :: fileHistoTemp, loadFileTemp, printFile
character(len = 128)    :: GnuPlotLoadFile

fileHistoTemp   = trim(fileHisto)//".dat"
loadFileTemp    = trim(loadFile)//".txt"
printFile       = trim(loadFile)//".eps"


open(unit=loadFileUnit,file=loadFileTemp,status='replace',iostat=ios)
    if (ios/=0) then
        print*,"Failed to open ",loadFile,"."
    end if
    !write(loadFileUnit,fmt=*) "set term wxt 3"

    ! titoli assi
    write(loadFileUnit,fmt=*) "set xlabel '",trim(xlabel),"'"
    write(loadFileUnit,fmt=*) "set ylabel '",trim(ylabel),"'"
    ! range in x
    write(loadFileUnit,fmt=*) "set xrange [",xmin,":",xmax,"]"
    ! range in y
    ! ymax è lasciato come libero, se lo si vuole utilizzare
    ! commentare la riga seguente e togliere il commento
    ! a quella successiva
    write(loadFileUnit,fmt=*) "set yrange [",ymin,":]" !",ymax,"]"
    !write(loadFileUnit,fmt=*) "set yrange [",ymin,":",ymax,"]"
    !personalizzazione istogramma
    write(loadFileUnit,fmt=*) "set style fill solid 0.5 "
    write(loadFileUnit,fmt=*) "set palette rgb 7,5,15"
    write(loadFileUnit,fmt=*) "set linetype 1 lc rgb 'black'"
    write(loadFileUnit,fmt=*) "set xtics out nomirror"
    write(loadFileUnit,fmt=*) "set ytics out nomirror"
    write(loadFileUnit,fmt=*) "set xtics font ',20'"
    write(loadFileUnit,fmt=*) "set ytics font ',20'"
    write(loadFileUnit,fmt=*) "set xlabel offset 18,-1"
    write(loadFileUnit,fmt=*) "set ylabel offset 0,5.8"
    write(loadFileUnit,fmt=*) "set xlabel font ',25'"
    write(loadFileUnit,fmt=*) "set ylabel font ',25'"
    write(loadFileUnit,fmt=*) "set term postscript color blacktext 'Helvetica' 24"
    write(loadFileUnit,fmt=*) "set key"
    ! togliere il commento alla riga successiva per avere una scala
    ! logaritmica in y
    !write(loadFileUnit,fmt=*) "set logscale y"
    ! per personalizzazione scala logaritmica in y togliere il commento alla
    ! riga successiva
    !write(loadFileUnit,fmt=*) "set format y '%2.0t{/Symbol \327}10^{%L}'"

    ! scrittura del file interpretabile da gnuplot
    write(loadFileUnit,fmt=*) "set output '",trim(printFile)
    write(loadFileUnit,fmt=*) "plot '",trim(fileHistoTemp),"' u 1:2 w boxes fs solid 0.5&
        fillcolor 'orange' title '",trim(histoTitle),"',&
        '' u 1:2:3 w error pt 7 ps 0.3 lt 1 notitle"
    write(loadFileUnit,fmt=*) "unset term"
    ! una volta eseguite le istruzioni, si esce da gnuplot.
    ! il file lo trovate già nella cartella dove si trova
    ! il programma.
    write(loadFileUnit,fmt=*) "q"

close(loadFileUnit)

! chiamata a gnuplot
GnuPlotLoadFile = "gnuplot "//trim(loadFileTemp)
call system(trim(GnuPlotLoadFile))

end subroutine drawHisto



end module HistoRoutines