; docformat = 'rst'
;
; NAME:
;    genoidReadObs
; PURPOSE:
;    Plot the distribution of residuals in astrometry between an asteroid and its satellite
;
;+
; :Description:
;    Plot the distribution of residuals in astrometry between an asteroid and its satellite
;
; :Categories:
;    Disk I/O, Genoid
;    
; :Params:
;    file: in, required, type=string
;        The path to the VOTable to be read. Fields are::
;          .JD:  Date in Julian Days  
;          .ISO: Date in ISO format
;          .refName: Name of the 
;          .refSys: 
;          .X: .Obs: X Position - Observed
;              .Err: X Position - Uncertainty
;              .Cal: X Position - Computed by Genoide
;              .OMC: X Position - Observed minus Computed
;          .Y: similar to X structure
;          .timeScale: Time scale (TT, UTC...)
;          .centerFrame: Reference frame center
;          .typeFrame: Type of reference frame
;          .coordType: Type of coordinates
;          .refFrame: Reference plane (1:equator, 2:ecliptic)
;          .obsIAU: IAU code of the observatory
;          .method: Photometric method used to obtain x,y measurements
;          .Mag: Difference in magnitude between primary and satellite
;          .dMag: Standard deviation on mag
;          .QC: Quality note of the measure (A|A+|B|C|D)
;          .telescope: Telescope used to acquire the measure
;          .camera: Camera used to acquire the measure
;          .filter: Filter used to acquire the measure
;          .author: Name of the person who measured the image
;          .when: Epoch at which the image was measured the last time (ISO format)
;          .bib: Bibliographic reference
;
; :Keywords:
;    dump: in, optional, type=string, default='./genoid-obs-omc.eps'
;      The path to the figure file
;
; :Example:
;    IDL> file='genoide.obsdiff.41_Daphne-S2008-41-1.xml'
;    IDL> genoidPlotObs, file, dump='41Daphne_omc'
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;     Change History::
;        Written in April 2017, B. Carry (OCA)
;-
pro genoidPlotObs, file, dump=dump
  COMPILE_OPT hidden, idl2
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--I.1-- Syntax Validation  -------------------------------------------------------------------;
  if not keyword_set(file) then begin
    message, /ioError, 'Syntax: genoidPlotObs, file [,dump=, /omc ]'
    return
  endif


 ;--I.2-- File Accessibility  ------------------------------------------------------------------;
  if not file_test(file,/read) then begin
    message, /Info, 'File cannot be read: '+strTrim(file,2)
    return
  endif


 ;--I.3-- Default Output File  -----------------------------------------------------------------;
  if not keyword_set(dump) then dump='./genoid-obs-omc.eps'


 ;--I.4-- Read the Observations  ---------------------------------------------------------------;
  obs = genoidReadObs(file)
  nbObs = n_elements(obs)


 ;--I.5-- Usefull Constants, Arrays, ...  ------------------------------------------------------;
  theta= 2*!PI*findgen(361)/360.


;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Define Graphical Parameters                        -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--II.1-- Size of the Figure  -----------------------------------------------------------------;
  xSize=30.
  ySize=30.

 ;--II.2-- Graphics Corners in Normalized Coordinates  -----------------------------------------;
  xBot = 0.08
  yBot = xBot
  xTop = 0.75
  yTop = xTop
  xEnd = 0.98
  yEnd = xEnd
  xSep = 0.01
  ySep = xSep

 ;--II.3-- Colorbar Corners in Normalized Coordinates  -----------------------------------------;
  xBCol = 0.72
  xTCol = 0.74
  yBCol = 0.12
  yTCol = 0.70

 ;--II.4-- Axes Ranges and Ticks  --------------------------------------------------------------;
  xyRange   = [-1,1] * 3.5
  binSize   = 0.25
  xyTickInt = 1
  xyMinor   = 4
  void      = replicate(' ', 40)

  yearTick  = 2

 ;--II.5-- Symbol Size and Colors  -------------------------------------------------------------;
  ctIndex = 31
  symSize = 1.5
  barFill = 0.80



;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---  III -- Analyze Observations and Prepare Graphics          -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--III.1-- Normalize OMC  ---------------------------------------------------------------------;
  xOMC = obs.x.OMC / obs.x.Err
  yOMC = obs.y.OMC / obs.y.Err


 ;--III.2-- Determine Time Range  --------------------------------------------------------------;
 ;--III.2.1-- Min & Max year
  trash = min(obs.jd, pMin, subscript_max=pMax)
  minYear = strMid( obs[pMin].iso,0,4)
  maxYear = string( float(strMid( obs[pMax].iso,0,4))+1, format='(I4)')
  minJD = date_conv( minYear+'-01-01', 'JULIAN')
  maxJD = date_conv( maxYear+'-01-01', 'JULIAN')
  if maxJD-minJD le (3*365) then yearTick=1

 ;--III.2.2-- Color Gradient
  xyCol = bytscl( obs.jd, min=minJD, max=maxJD)

  
 ;--III.3-- Compute OMC Histograms  ------------------------------------------------------------;
 ;--III.3.1-- Fix Number of Bins
  nbBin = round( (xyRange[1]-xyRange[0])/binSize + 1 )
 
 ;--III.3.2-- Compute Histograms
  xHist = histogram( xOMC, loc=xArr, min=xyRange[0], max=xyRange[1], binSize=binSize )
  yHist = histogram( yOMC, loc=yArr, min=xyRange[0], max=xyRange[1], binSize=binSize )

 ;--III.3.3-- Adjust Normal Distribuions
  xFit = mpFitPeak( xArr, xHist, xParam, nTerms=3 )
  yFit = mpFitPeak( yArr, yHist, yParam, nTerms=3 )
  fine = xyRange[0] + (xyRange[1]-xyRange[0])*findgen(1000)/999.
  xGauss=evalGauss( fine, xParam)
  yGauss=evalGauss( fine, yParam)

;  print, xParam
;  print, yParam




;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   IV -- Create Graphics                                    -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
 ;--II.1-- Open EPS  ---------------------------------------------------------------------------;
  cgPS_open, Filename=dump, /Metric, /Decomposed, /Portrait, $
             xSize=xSize, ysize=ySize, language_level=2, /Quiet
  
  
 ;--II.2-- Plot 2-D Distribution of OMC  -------------------------------------------------------;
 ;--II.2.1-- Plot Area
  cgPlot, xOMC, yOMC, /NoData, /Isotrop, /Normal, $
          xRange=xyRange, xTickInt=xyTickInt, xMinor=xyMinor, xStyle=1, $
          yRange=xyRange, yTickInt=xyTickInt, yMinor=xyMinor, yStyle=1, $
          position=[xBot, yBot, xTop, yTop], $
          xTitle = '( X!Do!N - X!Dc!N ) / '+cgGreek('sigma')+'!DX!N', $
          yTitle = '( Y!Do!N - Y!Dc!N ) / '+cgGreek('sigma')+'!DY!N'

 ;--II.2.2-- One-Two-Three Sigma Contours
  cgPlot, /OverPlot, 1*cos(theta), 1*sin(theta), color='Dark Gray'
  cgPlot, /OverPlot, 2*cos(theta), 2*sin(theta), color='Dark Gray'
  cgPlot, /OverPlot, 3*cos(theta), 3*sin(theta), color='Dark Gray'

 ;--II.2.3-- OMC Data Points
  cgLoadCt, ctIndex, /Brewer
  for kO=0, nbObs-1 do begin
    cgPlot, /OverPlot, xOMC[kO], yOMC[kO], psym='Filled circle', color=xyCol[kO], symSize=symSize
    cgPlot, /OverPlot, xOMC[kO], yOMC[kO], psym='Open circle', symSize=symSize
;    cgText, xOMC[kO], yOMC[kO], strMid(obs[kO].iso,0,10), charSize=1
  endfor


 ;--II.3-- Plot Color Bar  ---------------------------------------------------------------------;
  cgColorBar, position=[xBCol, yBCol, xTCol, yTCol], ctIndex=ctIndex, /Vertical, /Brewer, $
              yRange=float([minYear,maxYear]), yTickInt=yearTick
  


 ;--II.4-- Plot Vertical Histogram  ------------------------------------------------------------;
 ;--II.4.1-- Plot Area
  cgPlot, yHist, yArr, /NoData, /Normal, /NoErase, $
          yRange=xyRange, xTickInt=5, xMinor=5, $
          yTickName=void, yTickInt=xyTickInt, yMinor=xyMinor, $
          position=[xTop+xSep, yBot, xEnd, yTop]

 ;--II.4.2-- Reference sigma=1 Gaussian
  gParam=yParam
  gParam[0] = max(yHist)*0.8
;  gParam[0] /= 2.
  gParam[1:2] = [0, 1.]
  gGauss=evalGauss( fine, gParam)
  cgPlot, /OverPlot, gGauss, fine, lineStyle=2, color='Gray'
  cgColorFill, [gGauss,gGauss[0]], [fine,fine[0]], color='Light Gray'




 ;--II.4.3-- Draw Bars
  for kB=0, nbBin-1 do begin
    xx = yHist[kB] * [0,0,1,1,0]
    yy = yArr[kB] + binSize/2. + [-1,1,1,-1,-1]*BarFill*binSize/2.
    cgColorFill, xx, yy, color='Cornflower Blue'
    cgPlot, /OverPlot, xx, yy
  endfor

 ;--II.4.4-- Overplot Gaussian Approxmation
; cgPlot, /OverPlot, yGauss, fine, lineStyle=0, color='Dark Gray'



  
 ;--II.5-- Plot Horizontal Histogram  ----------------------------------------------------------;
 ;--II.5.1-- Plot Area
  cgPlot, xArr, xHist, /NoData, /Normal, /NoErase, $
          xStyle=1, xRange=xyRange, yTickInt=5, yMinor=5, $
          xTickInt=xyTickInt, xMinor=xyMinor, xTickName=void, $
          position=[xBot, yTop+xSep, xTop, yEnd]
    
 ;--II.5.2-- Reference sigma=1 Gaussian
  gParam = xParam
  gParam[0] = max(xHist)*0.8
  gParam[1:2] = [0, 1.]
  gGauss=evalGauss( fine, gParam)
  cgPlot, /OverPlot, fine, gGauss, lineStyle=2, color='Gray'
  cgColorFill, [fine,fine[0]], [gGauss,gGauss[0]], color='Light Gray'

 ;--II.5.3-- Draw Bars
  for kB=0, nbBin-1 do begin
    xx = xArr[kB] + binSize/2. + [-1,1,1,-1,-1]*BarFill*binSize/2.
    yy = xHist[kB] * [0,0,1,1,0]
    cgColorFill, xx, yy, color='Cornflower Blue'
    cgPlot, /OverPlot, xx, yy
  endfor

 ;--II.5.4-- Overplot Gaussian Approxmation
; cgPlot, /OverPlot, fine, xGauss, lineStyle=0, color='Dark Gray'
  
  

;---test gaussian sigma=1



 ;--II.6-- Close Graphic File  -----------------------------------------------------------------;
  cgPS_close, /png, Delete_PS=0
  

end