MODULE xg
  USE constants
  USE basic, ONLY : nr, nz, nrun, zgrid, rgrid, nplasma, itgraph, resfile
  USE beam, ONLY: partslist
  USE fields, ONLY: nrank, vec1, vec2
  IMPLICIT NONE
!
  REAL(kind=db) :: tstep=0.0_db, dnorm=1.0_db
  REAL(kind=db) :: zb1, zb2, rb1, rb2
  REAL(kind=db), ALLOCATABLE :: pzxg(:), prxg(:), pthetxg(:)
  INTEGER :: nhist, maxhist
  REAL(kind=db), ALLOCATABLE :: thist(:), tepot(:), tekin(:), tetot(:), terror(:)
  INTEGER :: x0=250, xsize = 410, y0=0, ysize=335
!
CONTAINS
  SUBROUTINE initw
!
!   Initialize Xgrafix
!
!  Particle phase space
!
    ALLOCATE (pzxg(partslist(1)%Nploc))
    ALLOCATE (prxg(partslist(1)%Nploc))
    ALLOCATE (pthetxg(partslist(1)%Nploc))

!
! Bounds of space domain
!
    zb1 = zgrid(0)
    zb2 = zgrid(nz)
    rb1 = rgrid(0)
    rb2 = rgrid(nr)
!
!  Time history
!
    nhist = 0
    maxhist=nrun/itgraph+1
    ALLOCATE(thist(0:maxhist))
    ALLOCATE(tepot(0:maxhist))
    ALLOCATE(tekin(0:maxhist))
    ALLOCATE(tetot(0:maxhist))
    ALLOCATE(terror(0:maxhist))
!
    CALL updt_xg_var
!
    CALL xginit(3,'fennecs:'//TRIM(resfile),'fennecs'//TRIM(resfile)//'.o','fennecs'//TRIM(resfile)//'.dump','N4','N5','N6',tstep)
!
!   Particle phase space plots
!
    pzxg(:) = partslist(1)%u(3,1:partslist(1)%Nploc)
    prxg(:) = partslist(1)%u(1,1:partslist(1)%Nploc)
    pthetxg(:) = partslist(1)%u(2,1:partslist(1)%Nploc)
!   Es potential
!
    CALL xgset2d('linlin','Z','Z-R phase space','open',x0, y0, dnorm, dnorm, &
         &        .FALSE., .FALSE., ZB1, ZB2, RB1, RB2)
    CALL xgscat2d(partslist(1)%pos(3,:),partslist(1)%pos(1,:),partslist(1)%Nploc,1)
!
    CALL xgset2d('linlin','PZ','PZ-PR phase space','open',x0+2*xsize, y0, 1._db, 1._db, &
         &        .TRUE., .TRUE., 0._db, 1._db, 0._db, 1._db)
    CALL xgscat2d(pzxg, prxg, partslist(1)%Nploc,1)
!
    CALL xgset2d('linlin','PZ','PZ-PTHET phase space','open',x0+2*xsize, y0+ysize, 1._db, 1._db, &
         &        .TRUE., .TRUE., 0._db, 1._db, 0._db, 1._db)
    CALL xgscat2d(pzxg, pthetxg, partslist(1)%Nploc,1)
!
    CALL xgset2d('linlin','Z','Z-PZ phase space','open',x0+xsize, y0, dnorm, 1.0_db, &
         &        .FALSE., .TRUE., ZB1, ZB2, 0._db, 1._db)
    CALL xgscat2d(partslist(1)%pos(3,:), pzxg, partslist(1)%Nploc,1)
!
!   Time history
!
    CALL xgset2d('linlin','T','Energies', 'open', x0, y0+ysize, 1._db, 1._db, &
         &       .TRUE., .TRUE., 0._db, 1._db, 0._db, 1._db)
    CALL xgcurve(thist, tepot, nhist, 1)
    CALL xgcurve(thist, tekin, nhist, 2)
    CALL xgcurve(thist, tetot, nhist, 3)
!
    CALL xgset2d('linlog','T','Relative error in Etot', 'open', x0+xsize, y0+ysize, 1._db, 1._db, &
         &       .TRUE., .TRUE., 0._db, 1._db, 0._db, 1._db)
    CALL xgcurve(thist, terror, nhist, 1)

!
!   Xgrafix control
!
    CALL xgevent
    CALL xgupdate
  END SUBROUTINE initw
!--------------------------------------------------------------------------------
  SUBROUTINE updt_xg_var
!
!   Update XG variables
!
    USE basic, ONLY : cstep, dt, tnorm
    USE beam, ONLY : ekin, epot, etot, etot0
!
!  Phase space
    tstep = REAL(cstep,db)
    pzxg(:) = partslist(1)%u(3,1:partslist(1)%Nploc)
    prxg(:) = partslist(1)%u(1,1:partslist(1)%Nploc)
    pthetxg(:) = partslist(1)%u(2,1:partslist(1)%Nploc)
!
!  Time history
!
    thist(nhist) = nhist*itgraph*dt*tnorm
    tekin(nhist) = ekin
    tepot(nhist) = epot
    tetot(nhist) = etot
    terror(nhist) = ABS((etot-etot0)/etot) ! To be dispalyed on log scale!
    nhist = nhist+1
!   

  END SUBROUTINE updt_xg_var
END MODULE xg