#ifndef _EVENTVARIABLE_H_
#define _EVENTVARIABLE_H_
///
///   @file	EventVar.h
///
///   @brief Header file for the EventVar Class, the main data container
///
///   @author John Burton <j.c.burton@gats-inc.com>
/// 
///   @date  Fri Mar 17 17:03:48 2006
///
///   @copyright
///   (&copy;) Copyright 2006 by GATS, Inc.,
///   11864 Canon Blvd, Suite 101, Newport News VA 23606
///
///   All Rights Reserved. No part of this software or publication may be
///   reproduced, stored in a retrieval system, or transmitted, in any form
///   or by any means, electronic, mechanical, photocopying, recording, or
///   otherwise without the prior written permission of GATS, Inc.
///
///   @class EventVar
///   @brief EventVar is the primary data container for the SOFIE level1 processing.
///   EventVar is derived from the std::valarray class and is used to contain
///   1-D data for the SOFIE level 1 processing. EventVar uses valarray's handling
///   of vector data (operators, methods & functions) and adds 6 properties
///   and several methods to help plot and store the vector data.
///

//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//
#include <math.h>
#include <iostream>
#include <fstream>
#include <valarray>
#include <numeric>
#include <algorithm>
#include <string>
#include <stdexcept>

#include "GracePlot.h"
//#include "GATS_Exception.h"
//#include "SOFIE_Exception.h"

//
//-----------------------------------------------------------------------
// Defines, Macros and Type Definitions:
//-----------------------------------------------------------------------
//

//using namespace std;


//
//-----------------------------------------------------------------------
// Class EventVar:
//-----------------------------------------------------------------------
//
template <class _Tp>
class EventVar: public std::valarray<_Tp> {
//
//-----------------------------------------------------------------------
// PublicProperties and Methods:
//-----------------------------------------------------------------------
//

public:
  /// default constructor - Constructs an empty EventVar object
  EventVar();

  /// Construct an EventVar object with \c n data elements
  /// @param n - number of data elements
  EventVar(int n);

  /// Construct a 2-D array EventVar object with \c rows * \c cols data elements
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(size_t rows, size_t cols);

  /// Copy constructor - creates a copy of \c d
  /// @param d - the EventVar object to be copied
  EventVar(const EventVar<_Tp>& d);

  /// Construct an EventVar object with \c n data elements initialized to 
  /// the value \c d
  /// @param d - initialize \c n data elements of EventVar object to this value
  /// @param n - number of data elements
  EventVar(const _Tp& d, int n);

  /// Construct a 2-D array EventVar object with \c rows * \c cols data elements 
  /// initialized to the value \c d
  /// @param d - initialize \c n data elements of EventVar object to this value
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const _Tp& d, size_t rows, size_t cols);

  /// Construct an eventVar object with \c n data elements initialized to
  /// the first \c n values pointed to by \c d
  /// @param d - pointer to the first of \c n data elements to be copied  
  /// @param n - number of data elements  
  EventVar(const _Tp* d, int n);

  /// Construct an eventVar object with data elements initialized to
  /// the values from \c t1 to \c t2 with an incerment of \c t3
  /// @param t1 - starting value 
  /// @param t2 - ending value 
  /// @param t3 - increment value 
  EventVar(const _Tp t1, const _Tp t2, const _Tp t3);

  /// Construct an eventVar object with \c rows*cols data elements 
  /// initialized to the first \c row*cols values pointed to by \c d
  /// @param d - pointer to the first of \c n data elements to be copied  
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const _Tp* d, size_t rows, size_t cols);

  /// Construct an empty EventVar object and initializes the property Name
  /// @param name - a string containing the name of the EventVar object
  EventVar(const std::string name);

  /// Construct an EventVar object with \c n elements and set the Name property to
  /// \c name
  /// @param name - a string containing the name of the EventVar object
  /// @param n - number of data elements
  EventVar(const std::string name, int n);

  /// Construct an EventVar object with \c rows * \c cols elements and 
  /// set the Name property to \c name
  /// @param name - a string containing the name of the EventVar object
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const std::string name, size_t rows, size_t cols);

  /// Copy constructor - creates a copy of \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param d  - The EventVar object to be copied
  EventVar(const std::string name, const EventVar<_Tp>& d);

  /// Copy constructor - creates a copy of \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param d  - The valarray object to be copied
  EventVar(const std::string name, const std::valarray<_Tp>& d);

  /// Copy constructor - creates a copy of \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param d  - The valarray object to be copied
  EventVar(const std::string name, const std::mask_array<_Tp>& d);

  /// Construct an EventVar object with \c n data elements initialized to 
  /// the value \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param d - initialize \c n data elements of EventVar object to this value
  /// @param n - number of data elements
  EventVar(const std::string name, const _Tp& d, int n);

  /// Construct an eventVar object with data elements initialized to
  /// the values from \c t1 to \c t2 with an incerment of \c t3
  /// @param name - a string containing the name of the EventVar object
  /// @param t1 - starting value 
  /// @param t2 - ending value 
  /// @param t3 - increment value 
  EventVar(const std::string name, const _Tp t1, const _Tp t2, const _Tp t3);

  /// Construct an EventVar object with \c rows * \c cols data elements 
  /// initialized to the value \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param d - initialize \c n data elements of EventVar object to this value
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const std::string name, const _Tp& d, size_t rows, size_t cols);

  /// Construct an eventVar object with \c n data elements 
  /// initialized to the first \c n values pointed to by \c d
  /// @param name - a string containing the name of the data in EventVar
  /// @param d - pointer to the first of \c n data elements to be copied  
  /// @param n - number of data elements
  EventVar(const std::string name, const _Tp* d, int n);

  /// Construct an eventVar object with \c rows * \c cols data elements 
  /// initialized to the first \c n values pointed to by \c d
  /// @param name - a string containing the name of the data in EventVar
  /// @param d - pointer to the first of \c n data elements to be copied  
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const std::string name, const _Tp* d, size_t rows, size_t cols);

  /// Construct an empty EventVar object and initializes the property Name
  /// @param name - a string containing the name of the data in EventVar
  /// @param desc - long (multi-word) description of the data in EventVar
  EventVar(const std::string name, const std::string desc);

  /// Construct an EventVar object with \c n elements and set the Name property to
  /// \c name
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param n - number of data elements
  EventVar(const std::string name, const std::string desc, int n);

  /// Construct an EventVar object with \c rows * \c cols elements and 
  /// set the Name property to \c name
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const std::string name, const std::string desc, size_t rows, size_t cols);

  /// Copy constructor - creates a copy of \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param d  - The EventVar object to be copied
  EventVar(const std::string name, const std::string desc, const EventVar<_Tp>& d);

  /// Construct an EventVar object with \c n data elements initialized to 
  /// the value \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param d - initialize \c n data elements of EventVar object to this value
  /// @param n - number of data elements
  EventVar(const std::string name, const std::string desc, const _Tp& d, int n);

  /// Construct an EventVar object with \c rows*cols data elements initialized to 
  /// the value \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param d - initialize \c n data elements of EventVar object to this value
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const std::string name, const std::string desc, const _Tp& d, size_t rows, size_t cols);

  /// Construct an EventVar object with \c n data elements initialized to 
  /// the value \c d
  /// @param name - a string containing the name of the data in EventVar
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param d - pointer to the first of \c n data elements to be copied  
  /// @param n - number of data elements
  EventVar(const std::string name, const std::string desc, const _Tp* d, int n);

  /// Construct an EventVar object with \c rows*cols data elements initialized to 
  /// the value \c d
  /// @param name - a string containing the name of the data in EventVar
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param d - pointer to the first of \c n data elements to be copied  
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const std::string name, const std::string desc, const _Tp* d, size_t rows, size_t cols);

  /// Construct an EventVar object with \c n elements and set the Name property to
  /// \c name
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param units - engineering units of the data
  EventVar(const std::string name, const std::string desc, 
	   const std::string units);

  /// Construct an EventVar object with \c n elements and set the Name property to
  /// \c name
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param units - engineering units of the data
  /// @param n - number of data elements
  EventVar(const std::string name, const std::string desc, 
	   const std::string units, int n);

  /// Construct an EventVar object with \c rows*cols elements and set the Name property to
  /// \c name
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param units - engineering units of the data
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const std::string name, const std::string desc, 
	   const std::string units, size_t rows, size_t cols);

  /// Copy constructor - creates a copy of \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param units - engineering units of the data
  /// @param d  - The EventVar object to be copied
  EventVar(const std::string name, const std::string desc, 
	   const std::string units, const EventVar<_Tp>& d);

  /// Construct an EventVar object with \c n data elements initialized to 
  /// the value \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param units - engineering units of the data
  /// @param d - initialize \c n data elements of EventVar object to this value
  /// @param n - number of data elements
  EventVar(const std::string name, const std::string desc, 
	   const std::string units, const _Tp& d, int n);

  /// Construct an EventVar object with \c rows*cols data elements initialized to 
  /// the value \c d
  /// @param name - a string containing the name of the EventVar object
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param units - engineering units of the data
  /// @param d - initialize \c rows*cols data elements of EventVar object to this value
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const std::string name, const std::string desc, 
	   const std::string units, const _Tp& d, size_t rows, size_t cols);

  /// Construct an EventVar object with \c n data elements initialized to 
  /// the value \c d
  /// @param name - a string containing the name of the data in EventVar
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param units - engineering units of the data
  /// @param d - pointer to the first of \c n data elements to be copied  
  /// @param n - number of data elements
  EventVar(const std::string name, const std::string desc, 
	   const std::string units, const _Tp* d, int n);

  /// Construct an EventVar object with \c rows*cols data elements initialized to 
  /// the value \c d
  /// @param name - a string containing the name of the data in EventVar
  /// @param desc - long (multi-word) description of the data in EventVar
  /// @param units - engineering units of the data
  /// @param d - pointer to the first of \c rows*cols data elements to be copied  
  /// @param rows - number of rows in the 2-D array
  /// @param cols - number of columns in the 2-D array
  EventVar(const std::string name, const std::string desc, 
	   const std::string units, const _Tp* d, size_t rows, size_t cols);

  /// EventVar destructor
  virtual ~EventVar();

  /// Set the value of the \c Name property
  /// @param name - a string containing the name of the data in EventVar
  void setName(std::string name);

  /// Set the value of the \c Description property
  /// @param desc - long (multi-word) description of the data in EventVar
  void setDescription(std::string desc);

  /// Set the value of the \c Units property
  /// @param units - engineering units of the data
  void setUnits(std::string units);

  /// Set the value of the \c AxisRange property based on the max and min
  /// data values
  void setAxisRange();

  /// Set the value of the \c AxisRange property
  /// @param lo - lower axis limit for plotting (AxisRange[0])
  /// @param hi - upper axis limit for plotting (AxisRange[1])
  void setAxisRange(_Tp lo, _Tp hi);

  /// Set the value of the \c FillValue property
  /// @param val - special value to indicate no data for this data element
  void setFillValue(_Tp val);

  /// Set the value of the \c MissingValue property
  /// @param val - special value to indicate no data for this data element
  void setMissingValue(_Tp val);

  /// Sets the number of rows in a 2-D array EventVar. Returns
  /// 0 (zero) if the EventVar is 1-D
  void setRows(size_t rows);

  /// Sets the number of columns in a 2-D array EventVar. Returns
  /// 0 (zero) if the EventVar is 1-D
  void setCols(size_t cols);

  /// Sets all the data elements to zero
  void setValueZero();

  /// Sets all the data elements to one
  void setValueOne();

  /// Sets all the data elements to one
  void setValueToMissing();

  /// Sets all the data elements to zero
  void setValueConst(_Tp d);


  /// Returns the value of the \c Name property
  std::string getName();

  /// Returns the value of the \c Description property
  std::string getDescription();

  /// Returns the value of the \c Units property
  std::string getUnits();

  /// Returns the value of the \c AxisRange property
  _Tp    *getAxisRange();

  /// Returns the value of the \c FillValue property
  _Tp     getFillValue();

  /// Returns the value of the \c MissingValue property
  _Tp     getMissingValue();

  /// Returns the number of rows in a 2-D array EventVar. Returns
  /// 0 (zero) if the EventVar is 1-D
  int   getRows();

  /// Returns the number of columns in a 2-D array EventVar. Returns
  /// 0 (zero) if the EventVar is 1-D
  int   getCols();

  /// Copies the data and the properties of one EventVar object to another
  /// @param ev - source EventVar object.
  void copy(EventVar<_Tp>& ev);

  /// Prints the properties and data elements fo EventVar to std::cerr
  void PrintVar();

  /// Prints the properties for EventVar to std::cerr
  void PrintVarProperties();

  /// Writes the properties and data elements for EventVar to file
  /// /c filename
  void writeVar(std::string filename);

  /// Reads the properties and data elements for EventVar from file
  /// /c filename
  void readVar(std::string filename);

  /// Assignment operator that assigns the values in a valarray to 
  /// the data elements of EventVar
  EventVar<_Tp>& operator = (const std::valarray<_Tp>& a);

  /// Assignment operator that assigns the values in a vector to 
  /// the data elements of EventVar
  EventVar<_Tp>& operator = (const std::vector<_Tp>& a);

  /// Assignment operator that assigns the values in a mask_array to 
  /// the data elements of EventVar
  EventVar<_Tp>& operator = (const std::mask_array<_Tp>& a);

  /// Assignment operator that assigns the values in a valarray to 
  /// the data elements of EventVar
  EventVar<_Tp>& operator = (const EventVar<_Tp>& a);

  /// Assignment operator that assigns the values in a valarray to 
  /// the data elements of EventVar
  template<class _Dom>
  EventVar<_Tp>& operator = (const std::_Expr<_Dom,_Tp> &a);

  /// 2-D Array index operator - only valid for 2-D array EventVars
  _Tp&  operator() (int r, int c) {return (*this)[r*getCols()+c];}

  /// 2-D Array index operator - only valid for 2-D array EventVars
  _Tp   operator() (int r, int c) const {return (*this)(r)[c];}

  /// Plots the data elements for EventVar using xmgrace. The EventVar 
  /// values are plotted as the Y (dependent) variable vs the index of 
  /// the data element as the X (independent) variable)
  /// @param Title - the title to be placed at the top of the chart
  /// @param SubTitle - the subtitle to be place at the top of the chart 
  /// below the title
  void Plot(std::string Title, std::string SubTitle);

  /// Plots the data elements for EventVar using xmgrace. The EventVar 
  /// values are plotted as the Y (dependent) variable vs the index of 
  /// the data element as the X (independent) variable)
  /// @param Title - the title to be placed at the top of the chart
  /// @param SubTitle - the subtitle to be place at the top of the chart 
  /// below the title
  /// @param graph - the graph object container for multiple plots.
  void Plot(std::string Title, std::string SubTitle, gracesc::XmgraceGraph &graph );

  /// Plots the data elements for EventVar using xmgrace. The EventVar 
  /// values are plotted as the Y (dependent) variable vs the values of 
  //!  the xvar EventVar as the X (independent) variable)
  /// @param Title - the title to be placed at the top of the chart
  /// @param SubTitle - the subtitle to be place at the top of the chart 
  /// below the title
  /// @param xvar - the X (independent) variable values to be used
  void Plot(std::string Title, std::string SubTitle, const EventVar<_Tp>& xvar);
  
  /// Interpolates within the \c yin EventVar to determine the output value 
  /// that corresponds to \c xout within the EventVar object. The eventVar object
  /// and EventVar \c yin must be the same length.The EventVar object (\c this) 
  /// represents a series of \c X values. The EventVar object \c yin represents
  /// the corresponding series of \c Y values. Interpol determines the Y value 
  /// corresponding to \c xout.
  /// @param yin - the EventVar object containing a series of \c Y values
  /// @param xout - determine the Y value that corresponds to this.
  _Tp    Interpol(EventVar<_Tp>& yin, _Tp xout);

  /// Interpolates within the index range of the EventVar to determine the output value 
  /// that corresponds to \c xout within the EventVar object. .The EventVar object (\c this) 
  /// represents a series of \c Y values, with the indices of the EventVar object
  /// representing the corresponding \c X values. Interpol determines the Y value 
  /// corresponding to \c xout.
  /// @param xout - determine the Y value that corresponds to this.
  _Tp    Interpol_ndx(const _Tp xout);

  EventVar<_Tp> VInterpol(EventVar<_Tp>& yin, EventVar<_Tp>& xout);
  EventVar<_Tp> VInterpol(EventVar<_Tp>& yin, EventVar<_Tp>& xout, int ext);

  EventVar<_Tp> VInterpol_ndx(EventVar<_Tp>& xout);

  int   Vlocate(_Tp x);
  int   Nearest(_Tp x);

  EventVar<_Tp> Smooth(int M);
  double Mean(void);
  double Stddev(void);
  EventVar<_Tp> Stddev(int M);

//  void resize(int n) {this->resize(n);}

//
//-----------------------------------------------------------------------
// Private Properties and Methods:
//-----------------------------------------------------------------------
//

private:
  std::string Name;        //!< Name of the Variable - will be used by Event
                           //!<  object
  std::string Description; //!< Long description of Variable - will be used 
                           //!<  for plotting
  std::string Units;       //!< Units for the variable
  _Tp      AxisRange[2];      //!< Range of values to plot
  _Tp      FillValue;         //!< Fill the data area with this value
  _Tp      MissingValue;      //!< This value to be used for missing data elements
  size_t Rows;             //!< Number of Rows in a 2-D EventVar
  size_t Cols ;             //!< Number of Columns in a 2-D EventVar
//  gracesc::Xmgrace* graph;
  void  init();
  _Tp     linterp(_Tp x, _Tp xa, _Tp xb, _Tp ya, _Tp yb);
//  int   Vlocate(_Tp x);

}; // End of EventVar Class definition


//
//-----------------------------------------------------------------------
// Inline Methods:
//-----------------------------------------------------------------------
//

template <class _Tp>
inline void EventVar<_Tp>::init()
{
  Rows         = 0;
  Cols         = 0;
  Name         = "";
  Description  = "";
  Units        = "";
  AxisRange[0] = 0;
  AxisRange[1] = 0;
  FillValue    = (_Tp)-1e24;
  MissingValue = (_Tp)-1e24;
//  FillValue    = (_Tp)-9999;
//  MissingValue = (_Tp)-9999;
}


//
// Constructors 
//

template <class _Tp>
inline EventVar<_Tp>::EventVar() : std::valarray<_Tp>()
{ 
  init(); 
}


template <class _Tp>
inline EventVar<_Tp>::EventVar(int n) : std::valarray<_Tp>(n)
{ 
  init(); 
}


template <class _Tp>
inline EventVar<_Tp>::EventVar(size_t r, size_t c) : std::valarray<_Tp>(r*c)
{ 
  init(); 
  Rows  = r;
  Cols  = c;
}


template <class _Tp>
inline EventVar<_Tp>::EventVar(const EventVar<_Tp>& d) : std::valarray<_Tp>(d)
{ 
  Rows         = d.Rows;
  Cols         = d.Cols;
  Name         = d.Name;
  Description  = d.Description;
  Units        = d.Units;
  AxisRange[0] = d.AxisRange[0];
  AxisRange[1] = d.AxisRange[1];
  FillValue    = d.FillValue;
  MissingValue = d.MissingValue;
}


//
// Constructors for  Name 
//

template <class _Tp>
inline EventVar<_Tp>::EventVar(const _Tp& d, int n) : std::valarray<_Tp>(d,n)
{ 
  init(); 
}


template <class _Tp>
inline EventVar<_Tp>::EventVar(const _Tp* d, int n) : std::valarray<_Tp>(d,n)
{ 
  init(); 
}


template <class _Tp>
inline EventVar<_Tp>::EventVar(const _Tp t1, const _Tp t2, const _Tp t3) : std::valarray<_Tp>((int)((t2-t1)/t3)+1)
{
  int i;
  _Tp t;
  int n = (int)((t2-t1)/t3)+1;

  init();
  t = t1;
  for(i=0;i<n;i++)
  {
    (*this)[i] = t;
    t += t3;
  }
}



template <class _Tp>
inline EventVar<_Tp>::EventVar(const _Tp* d, size_t rows, size_t cols) : std::valarray<_Tp>(d,rows*cols)
{ 
  init(); 
  Rows  = rows;
  Cols  = cols;
}



template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name) : std::valarray<_Tp>()
{
  init();
  Name = name;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, int n) : std::valarray<_Tp>(n)
{
  init();
  Name = name;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			       size_t rows, size_t cols) : std::valarray<_Tp>(rows*cols)
{
  init();
  Rows  = rows;
  Cols  = cols;
  Name  = name;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name,
			       const _Tp t1, const _Tp t2, const _Tp t3) : std::valarray<_Tp>(((int)((t2-t1)/t3)+1))
{
  int i;
  _Tp t;
  int n = (int)((t2-t1)/t3)+1;

  init();
  t = t1;
  for(i=0;i<n;i++)
  {
    (*this)[i] = t;
    t += t3;
  }
  Name = name;
}


template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			       const _Tp& d, int n) : std::valarray<_Tp>(d,n)
{
  init();
  Name = name;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			       const _Tp& d, size_t rows, size_t cols) : std::valarray<_Tp>(d,rows*cols)
{
  init();
  Rows  = rows;
  Cols  = cols;
  Name  = name;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const _Tp* d, int n) : std::valarray<_Tp>(d,n)
{
  init();
  Name = name;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const _Tp* d, size_t rows, size_t cols) : std::valarray<_Tp>(d,rows*cols)
{
  init();
  Rows  = rows;
  Cols  = cols;
  Name  = name;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const EventVar<_Tp>& d) : std::valarray<_Tp>(d)
{
  Name         = name;
  Rows         = d.Rows;
  Cols         = d.Cols;
  Description  = d.Description;
  Units        = d.Units;
  AxisRange[0] = d.AxisRange[0];
  AxisRange[1] = d.AxisRange[1];
  FillValue    = d.FillValue;
  MissingValue = d.MissingValue;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::valarray<_Tp>& d) : std::valarray<_Tp>(d)
{
  init();
  Name         = name;
}


template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::mask_array<_Tp>& d) : std::valarray<_Tp>(d)
{
  init();
  Name         = name;
}


//
// Constructors for  Name and Description
//

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc) : std::valarray<_Tp>()
{
  init();
  Name        = name;
  Description = desc;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     int n) : std::valarray<_Tp>(n)
{
  init();
  Name        = name;
  Description = desc;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     size_t rows, size_t cols) : std::valarray<_Tp>(rows*cols)
{
  init();
  Rows        = rows;
  Cols        = cols;
  Name        = name;
  Description = desc;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const _Tp& d, int n) : std::valarray<_Tp>(d,n)
{
  init();
  Name        = name;
  Description = desc;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const _Tp& d, size_t r, size_t c) : std::valarray<_Tp>(d,r*c)
{
  init();
  Rows        = r;
  Cols        = c;
  Name        = name;
  Description = desc;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const _Tp* d, int n) : std::valarray<_Tp>(d,n)
{
  init();
  Name        = name;
  Description = desc;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const _Tp* d, size_t rows, size_t cols) : std::valarray<_Tp>(d,rows*cols)
{
  init();
  Rows        = rows;
  Cols        = cols;
  Name        = name;
  Description = desc;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const EventVar<_Tp>& d) : std::valarray<_Tp>(d)
{
  Name         = name;
  Description  = desc;
  Rows         = d.Rows;
  Cols         = d.Cols;
  Units        = d.Units;
  AxisRange[0] = d.AxisRange[0];
  AxisRange[1] = d.AxisRange[1];
  FillValue    = d.FillValue;
  MissingValue = d.MissingValue;
}


//
// Constructors for Name, Description, and Units
//

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const std::string units) : std::valarray<_Tp>()
{
  init();
  Name        = name;
  Description = desc;
  Units       = units;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const std::string units, 
			     int n) : std::valarray<_Tp>(n)
{
  init();
  Name        = name;
  Description = desc;
  Units       = units;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const std::string units, 
			     size_t rows, size_t cols) : std::valarray<_Tp>(rows*cols)
{
  init();
  Rows        = rows;
  Cols        = cols;
  Name        = name;
  Description = desc;
  Units       = units;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const std::string units, 
			     const _Tp& d, int n) : std::valarray<_Tp>(d,n)
{
  init();
  Name        = name;
  Description = desc;
  Units       = units;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const std::string units, 
			     const _Tp& d, size_t rows, size_t cols) : std::valarray<_Tp>(d,rows*cols)
{
  init();
  Rows        = rows;
  Cols        = cols;
  Name        = name;
  Description = desc;
  Units       = units;
}

template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const std::string units, 
			     const _Tp* d, int n) : std::valarray<_Tp>(d,n)
{
  init();
  Name        = name;
  Description = desc;
  Units       = units;
}



template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const std::string units, 
			     const _Tp* d, size_t rows, size_t cols) : std::valarray<_Tp>(d,rows*cols)
{
  init();
  Rows        = rows;
  Cols        = cols;
  Name        = name;
  Description = desc;
  Units       = units;
}



template <class _Tp>
inline EventVar<_Tp>::EventVar(const std::string name, 
			     const std::string desc, 
			     const std::string units, 
			     const EventVar<_Tp>& d) : std::valarray<_Tp>(d)
{
  Name         = name;
  Description  = desc;
  Units        = units;
  AxisRange[0] = d.AxisRange[0];
  AxisRange[1] = d.AxisRange[1];
  FillValue    = d.FillValue;
  MissingValue = d.MissingValue;
}


//
// Destructor
//

template <class _Tp>
inline EventVar<_Tp>::~EventVar()
{
//  Name.clear();
//  Description.clear();
//  Units.clear();
//  Name        = "";
//  Description = "";
//  Units       = "";
//  this->std::valarray<_Tp>::~valarray<_Tp>();
}


//
// Set Property Methods
//

template <class _Tp>
inline void EventVar<_Tp>::setName(std::string name)
{
  Name = name;
}


template <class _Tp>
inline void EventVar<_Tp>::setDescription(std::string desc)
{
  Description = desc;;
}


template <class _Tp>
inline void EventVar<_Tp>::setUnits(std::string units)
{
  Units = units;
}


template <class _Tp>
inline void EventVar<_Tp>::setAxisRange()
{
  AxisRange[0] = this->min();
  AxisRange[1] = this->max();
}


template <class _Tp>
inline void EventVar<_Tp>::setAxisRange(_Tp lo, _Tp hi)
{
  AxisRange[0] = lo;
  AxisRange[1] = hi;
}


template <class _Tp>
inline void EventVar<_Tp>::setFillValue(_Tp val)
{
  FillValue = val;
}


template <class _Tp>
inline void EventVar<_Tp>::setMissingValue(_Tp val)
{
  MissingValue = FillValue = val;
}

template <class _Tp>
inline void EventVar<_Tp>::setRows(size_t rows)
{
  Rows = rows;
}

template <class _Tp>
inline void EventVar<_Tp>::setCols(size_t cols)
{
  Cols = cols;
}

template <class _Tp>
inline void EventVar<_Tp>::setValueZero()
{
  int i;
  for(i=0;i<this->size();i++)
    (*this)[i] = (_Tp)0;
}

template <class _Tp>
inline void EventVar<_Tp>::setValueOne()
{
  int i;
  for(i=0;i<this->size();i++)
    (*this)[i] = (_Tp)1;
}

template <class _Tp>
inline void EventVar<_Tp>::setValueToMissing()
{
  int i;
  for(i=0;i<this->size();i++)
    (*this)[i] = MissingValue;
}

template <class _Tp>
inline void EventVar<_Tp>::setValueConst(_Tp d)
{
  int i;
  for(i=0;i<this->size();i++)
    (*this)[i] = (_Tp)d;
}



//
// Get Property Methods
//


template <class _Tp>
inline std::string EventVar<_Tp>::getName()
{
  return Name;
}


template <class _Tp>
inline std::string EventVar<_Tp>::getDescription()
{
  return Description;
}


template <class _Tp>
inline std::string EventVar<_Tp>::getUnits()
{
  return Units;
}


template <class _Tp>
inline _Tp *EventVar<_Tp>::getAxisRange()
{
  return AxisRange;
}


template <class _Tp>
inline _Tp EventVar<_Tp>::getFillValue()
{
  return FillValue;
}


template <class _Tp>
inline _Tp EventVar<_Tp>::getMissingValue()
{
  return MissingValue;
}

template <class _Tp>
inline int EventVar<_Tp>::getRows()
{
  return Rows;
}

template <class _Tp>
inline int EventVar<_Tp>::getCols()
{
  return Cols;
}


//
// Utility Methods
//

template <class _Tp>
void EventVar<_Tp>::copy(EventVar<_Tp> &ev)
{
  *this        = ev;
  Name         = ev.Name;
  Description  = ev.Description;
  Units        = ev.Units;
  AxisRange[0] = ev.AxisRange[0];
  AxisRange[1] = ev.AxisRange[1];
  FillValue    = ev.FillValue;
  MissingValue = ev.MissingValue;
}


template <class _Tp>
inline void EventVar<_Tp>::PrintVarProperties()
{
  int i;
  int last = this->size() - 1;

  std::cerr << "EventVar::PrintVarProperties\n";
  std::cerr << "Name         = " << this->getName() << "\n";
  std::cerr << "Description  = " << this->getDescription() << "\n";
  std::cerr << "Units        = " << this->getUnits() << "\n";
  std::cerr << "AxisRange[0] = " << this->getAxisRange()[0] << "\n";
  std::cerr << "AxisRange[1] = " << this->getAxisRange()[1] << "\n";
  std::cerr << "FillValue    = " << this->getFillValue() << "\n";
  std::cerr << "MissingValue = " << this->getMissingValue() << "\n";

}


template <class _Tp>
inline void EventVar<_Tp>::PrintVar()
{
  int i;
  int last = this->size() - 1;

  std::cerr << "EventVar::PrintVar\n";
  std::cerr << "Name         = " << this->getName() << "\n";
  std::cerr << "Description  = " << this->getDescription() << "\n";
  std::cerr << "Units        = " << this->getUnits() << "\n";
  std::cerr << "AxisRange[0] = " << this->getAxisRange()[0] << "\n";
  std::cerr << "AxisRange[1] = " << this->getAxisRange()[1] << "\n";
  std::cerr << "FillValue    = " << this->getFillValue() << "\n";
  std::cerr << "MissingValue = " << this->getMissingValue() << "\n";

  std::cerr << " " << this->size() << " ";
  std::cerr << " [";
  for(i=0;i<last;i++)
    std::cerr << (*this)[i] << ",";
  std::cerr << (*this)[last] << "]\n";
}



template <class _Tp>
inline void EventVar<_Tp>::writeVar(std::string filename)
{
  std::ofstream outfile(filename.c_str());

  if(outfile.bad()){
    std::string errstr = "Cannot Open Output FIle " + std::string(filename);
    throw std::runtime_error(errstr);
  }

  int i;

  outfile << this->getName() << "\n";
  outfile << this->getDescription() << "\n";
  outfile << this->getUnits() << "\n";
  outfile << this->getAxisRange()[0] << "\n";
  outfile << this->getAxisRange()[1] << "\n";
  outfile << this->getFillValue() << "\n";
  outfile << this->getMissingValue() << "\n";

  outfile << this->size() << "\n";
  outfile << this->getRows() << "\n";
  outfile << this->getCols() << "\n";

  outfile.setf(std::ios::scientific, std::ios::floatfield);
  outfile.precision(22);

  for(i=0;i<this->size();i++)
    outfile << (*this)[i] << "\n";
  outfile.close();
}


template <class _Tp>
inline void EventVar<_Tp>::readVar(std::string filename)
{
  std::ifstream infile(filename.c_str());

  if(infile.bad()){
    std::string errstr = "Cannot open input file: " + std::string(filename);
    throw std::runtime_error(errstr);
  }

  std::string buff;
  _Tp temp;

  int i;
  int size,rows,cols;

  getline(infile,buff);
  this->setName(buff);
  getline(infile,buff);
  this->setDescription(buff);
  getline(infile,buff);
  this->setUnits(buff);
  infile >> AxisRange[0];
  infile >> AxisRange[1];
  infile >> FillValue;
  infile >> MissingValue;
  infile >> size;
  infile >> Rows;
  infile >> Cols;
  this->resize(size);
  for(i=0;i<size;i++){
    infile >> (*this)[i];
  }
  infile.close();

}


template <class _Tp>
EventVar<_Tp>  Range(_Tp t1, _Tp t2, _Tp t3)
{
  int i;
  _Tp t;
  int n = (int)((t2-t1)/t3)+1;

  _Tp ret[n];
//  std::valarray<_Tp> ret(n);
  t = t1;
  for(i=0;i<n;i++)
  {
    ret[i] = t;
    t += t3;
  }
  return EventVar<_Tp>(ret,n);
}


template <class _Tp>
inline void EventVar<_Tp>::Plot(std::string Title, std::string SubTitle)
{
  int i;
  int nl = this->size();
  std::string xlabel = "Index";
  std::string ylabel = Description + "(" + Units + ")";
  _Tp delta;

  gracesc::Xmgrace graph(1);
  for(i=0;i<nl;i++)
  {
    delta = fabs((*this)[i] - MissingValue);
    if(delta > 0.001)
      graph[0] << i << (*this)[i] << gracesc::endl;
  }
  graph[0].draw();
//  graph[0].Scale_xAxis(0,0,nl);
//  graph[0].Scale_yAxis(0,AxisRange[0],AxisRange[1]);
  graph[0].Display_titles(Title,SubTitle,xlabel,ylabel);
  graph[0].Set_legend_location("1.0","0.8");
  graph[0].SetLegend(0,Name);
  graph[0].autoscale();
  graph[0].timestamp();
  graph.redraw();
}


template <class _Tp>
inline void EventVar<_Tp>::Plot(std::string Title, std::string SubTitle,gracesc::XmgraceGraph &graph )
{
  int i;
  int nl = this->size();
  std::string xlabel = "Index";
  std::string ylabel = Description + "(" + Units + ")";
  _Tp delta;

//  gracesc::Xmgrace graph(1);
  for(i=0;i<nl;i++)
  {
    delta = fabs((*this)[i] - MissingValue);
    if(delta > 0.001)
      graph << i << (*this)[i] << gracesc::endl;
  }
  graph.draw();
//  graph.Scale_xAxis(0,0,nl);
//  graph.Scale_yAxis(0,AxisRange[0],AxisRange[1]);
  graph.Display_titles(Title,SubTitle,xlabel,ylabel);
  graph.Set_legend_location("1.0","0.8");
  graph.SetLegend(0,Name);
  graph.autoscale();
  graph.timestamp();
//  graph.redraw();
}

  
template <class _Tp>
inline void EventVar<_Tp>::Plot(std::string Title, std::string SubTitle, 
				 const EventVar<_Tp>& xvar)
{
  int i;
  int nl = this->size(); 
  std::string xlabel = xvar.Description + " (" + xvar.Units + ")";
  std::string ylabel = Description + " (" + Units + ")";
  _Tp delta;
//  if(nl != xvar.size())
//    error_return();

  gracesc::Xmgrace graph(1);
  for(i=0;i<nl;i++)
  {
    delta = fabs((*this)[i] - MissingValue);
    if(delta > 0.001)
      graph[0] << xvar[i] << (*this)[i] << gracesc::endl;
  }
  graph[0].draw();
//  graph[0].Scale_xAxis(0,xvar.AxisRange[0],xvar.AxisRange[1]);
//  graph[0].Scale_yAxis(0,AxisRange[0],AxisRange[1]);
  graph[0].Display_titles(Title,SubTitle,xlabel,ylabel);
  graph[0].Set_legend_location("1.0","0.8");
  graph[0].SetLegend(0,Name);
  graph[0].autoscale();
  graph[0].timestamp();
  graph.redraw();
}


    
template <class _Tp>
inline EventVar<_Tp>& EventVar<_Tp>::operator = (const std::valarray<_Tp>& a) 
{

  if(this->size() != a.size())
    this->resize(a.size());
  std::valarray<_Tp>::operator=(a);
  return *this;
}
    
template <class _Tp>
inline EventVar<_Tp>& EventVar<_Tp>::operator = (const std::vector<_Tp>& a) 
{
  if(this->size() != a.size())
    this->resize(a.size());
  for(int i = 0; i < a.size(); i++)
    (*this)[i] = a[i];
  return *this;
}
    

template <class _Tp>
inline EventVar<_Tp>& EventVar<_Tp>::operator = (const std::mask_array<_Tp>& a) 
{
  std::valarray<_Tp> tmp = a;

  if(this->size() != tmp.size())
    this->resize(tmp.size());
  std::valarray<_Tp>::operator=(tmp);
  return *this;
}
    

template <class _Tp>
inline EventVar<_Tp>& EventVar<_Tp>::operator = (const EventVar<_Tp>& a) 
{
  Name         = a.Name;
  Description  = a.Description;
  Units        = a.Units;
  AxisRange[0] = a.AxisRange[0];
  AxisRange[1] = a.AxisRange[1];
  FillValue    = a.FillValue;
  MissingValue = a.MissingValue;

  if(this->size() != a.size())
    this->resize(a.size());
  std::valarray<_Tp>::operator=(a);
  return *this;
}
    

template <class _Tp> template<class _Dom>
inline EventVar<_Tp>& EventVar<_Tp>::operator = (const std::_Expr<_Dom, _Tp>& e) 
{
  if(this->size() != e.size())
    this->resize(e.size());
  std::valarray<_Tp>::operator=(e);
  return *this;
}
    


template <class _Tp>
inline int EventVar<_Tp>::Vlocate(_Tp x)
{
  int asc;
  int hi, mid, lo;

  lo  = -1;
  hi  = this->size();
  asc = ( (*this)[hi-1] > (*this)[0]);

  while((hi - lo) > 1)
  {
    mid = (hi + lo) / 2;
    if ((x > (*this)[mid]) == asc)
      lo = mid;
    else
      hi = mid;
  }
  return lo;
}


template <class _Tp>
inline int EventVar<_Tp>::Nearest(_Tp x)
{
  int ndx,len;
  _Tp d1,d2;
  int hi, lo;

  lo  = this->Vlocate(x);
  len = this->size();
  if(lo <= 0)
    lo = 0;
  if(lo >= len-1)
    lo = len-1;
//    lo = len-2;
//  hi = lo + 1;
//  d1 = abs((*this)[lo] - x);
//  d2 = abs((*this)[hi] - x);
//  if (d1 < d2)
//    ndx = lo;
//  else
//    ndx = hi;
  
//  return ndx;
  return lo;
}


template <class _Tp>
inline _Tp EventVar<_Tp>::linterp(_Tp x, _Tp xa, _Tp xb, _Tp ya, _Tp yb)
{
  _Tp dx, y;
  _Tp a;

  if((xa > MissingValue) && (xb > MissingValue) &&
     (ya > MissingValue) && (yb > MissingValue))
  {
    dx = xb - xa;
    if(dx != 0)
    {
      a = (x - xa) / dx;
      y = ya + a * (yb - ya);
    }else
      y = (ya * yb) / 2.0;
  }
  else
    y = MissingValue;
  return y;
}

template <class _Tp>
inline _Tp EventVar<_Tp>::Interpol(EventVar<_Tp>& xin, _Tp xout)
{
  int ndx1, ndx2, len;
  _Tp   yout;

  len  = xin.size();
  ndx1 = xin.Vlocate(xout);

  if(ndx1 <= 0)
    ndx1 = 0;
  if(ndx1 >= len-1)
    ndx1 = len-2;

  ndx2 = ndx1 + 1;

  yout = linterp(xout,xin[ndx1],xin[ndx2],(*this)[ndx1],(*this)[ndx2]);

  return yout;
}


template <class _Tp>
inline EventVar<_Tp> EventVar<_Tp>::VInterpol(EventVar<_Tp>& xin, EventVar<_Tp>& xout)
{
  int i;
  int cnt = xout.size();
  EventVar<_Tp> yout(cnt);
  _Tp xmin,xmax;
//  std::cerr << "cnt = " << cnt << "\n";

  yout.Name         = xout.Name;
  yout.Description  = xout.Description;
  yout.Units        = xout.Units;
  yout.AxisRange[0] = xout.AxisRange[0];
  yout.AxisRange[1] = xout.AxisRange[1];
  yout.FillValue    = xout.FillValue;
  yout.MissingValue = xout.MissingValue;
  xmin = xin.min();
  xmax = xin.max();

  for(i=0;i<cnt;i++)
  {
    if((xout[i] < xmin) || (xout[i] > xmax))
      yout[i] = MissingValue;
    else
      yout[i] = this->Interpol(xin,xout[i]);
  }
  return yout;
}


template <class _Tp>
inline EventVar<_Tp> EventVar<_Tp>::VInterpol(EventVar<_Tp>& xin, EventVar<_Tp>& xout, int ext)
{
  int i;
  int cnt = xout.size();
  EventVar<_Tp> yout(cnt);
  EventVar<_Tp> xin_tmp,this_tmp;
  _Tp xmin,xmax;
  _Tp mv = xin.getMissingValue();

  std::valarray<bool> valid = (((*this) > MissingValue) && (xin > mv));

  xin_tmp  = xin[valid];
  this_tmp = (*this)[valid];

  yout.Name         = xout.Name;
  yout.Description  = xout.Description;
  yout.Units        = xout.Units;
  yout.AxisRange[0] = xout.AxisRange[0];
  yout.AxisRange[1] = xout.AxisRange[1];
  yout.FillValue    = xout.FillValue;
  yout.MissingValue = xout.MissingValue;
  xmin = xin_tmp.min();
  xmax = xin_tmp.max();

  if(ext > 0)
  {
    for(i=0;i<cnt;i++)
    {
//      if((xout[i] < xmin) || (xout[i] > xmax))
//	yout[i] = MissingValue;
//      else
	yout[i] = this_tmp.Interpol(xin_tmp,xout[i]);
    }
  }
  else
  {
    for(i=0;i<cnt;i++)
    {
      if((xout[i] < xmin) || (xout[i] > xmax))
	yout[i] = MissingValue;
      else
	yout[i] = this_tmp.Interpol(xin_tmp,xout[i]);
    }
  }

  return yout;
}

template <class _Tp>
inline _Tp EventVar<_Tp>::Interpol_ndx(const _Tp xout)
{
  int ndx1, ndx2, len;
  _Tp   yout;

  len  = this->size();
  ndx1 = (int)xout;

  if(ndx1 <= 0)
    ndx1 = 0;
  if(ndx1 >= len-1)
    ndx1 = len-2;

  ndx2 = ndx1 + 1;

  yout = linterp(xout,(_Tp)ndx1,(_Tp)ndx2,(*this)[ndx1],(*this)[ndx2]);

  return yout;
}



template <class _Tp>
inline EventVar<_Tp> EventVar<_Tp>::VInterpol_ndx(EventVar<_Tp>& xout)
{
  int i;
  int cnt = xout.size();
  EventVar<_Tp> yout(xout);
  _Tp xmin,xmax;

//  yout.Name         = xout.Name;
//  yout.Description  = xout.Description;
//  yout.Units        = xout.Units;
//  yout.AxisRange[0] = xout.AxisRange[0];
//  yout.AxisRange[1] = xout.AxisRange[1];
//  yout.FillValue    = xout.FillValue;
//  yout.MissingValue = xout.MissingValue;
//  xmin = this->min();
//  xmax = this->max();
  xmin = (_Tp)0;
  xmax = (_Tp)this->size();

//  std::cerr << "Name = " << Name << "\n";
//  std::cerr << "cnt  = " << cnt  << "\n";
//  std::cerr << "xmin = " << xmin << "\n";
//  std::cerr << "xmax = " << xmax << "\n";
  for(i=0;i<cnt;i++)
  {
    if((xout[i] < xmin) || (xout[i] > xmax))
      yout[i] = MissingValue;
    else
      yout[i] = this->Interpol_ndx(xout[i]);
  }
  return yout;
}



template <class _Tp>
inline EventVar<_Tp> EventVar<_Tp>::Smooth(int M)
{
  _Tp  acc;
  int i,j,k;
  int M2;
  int N  = this->size();
  EventVar<_Tp> yout(N);

  M  = ((M%2) == 0) ? (M-1) : M;
  M2 = M / 2;
  for(i = 0; i < N; i++)
  {
    acc = (_Tp)0;
    for(j = -M2; j<= M2; j++)
    {
      k = (j + i);
      k = (k < 0) ? 0 : ((k > (N-1)) ? (N-1) : k);
      acc += (*this)[k];
    }
    yout[i] = acc / (_Tp)M;
  }
  return yout;
}


template <class _Tp>
inline double EventVar<_Tp>::Mean(void)
{
  double sum = 0.0;
  double mean;
  int len = this->size();
  int i;

  for(i=0;i<len; i++)
    sum += (double)(*this)[i];

  mean = sum / (double)len;
  return mean;
}

template <class _Tp>
inline EventVar<_Tp> EventVar<_Tp>::Stddev(int M)
{
  _Tp  acc;
  double sum,mean,tmp,stddev;
  int i,j,k;
  int M2,cnt;
  int N  = this->size();
  EventVar<_Tp> yout(N);
  std::valarray<bool> valid = ((*this) > MissingValue);

  M  = ((M%2) == 0) ? (M-1) : M;
  M2 = M / 2;
  for(i = 0; i < N; i++)
  {
    cnt = 0;
    sum = 0.0;
    for(j = -M2; j<= M2; j++)
    {
      k = (j + i);
      k = (k < 0) ? 0 : ((k > (N-1)) ? (N-1) : k);
      if(valid[k])
      {
	sum += (double)(*this)[k];
	cnt++;
      }
    }
    mean = sum / (double) cnt;
    sum = 0.0;
    cnt = 0;
    for(j = -M2; j<= M2; j++)
    {
      k = (j + i);
      k = (k < 0) ? 0 : ((k > (N-1)) ? (N-1) : k);
      if(valid[k])
      {
	tmp = (double)(*this)[k] - mean;
	sum = sum + (tmp * tmp);
	cnt++;
      }
    }
    tmp = sum / (double) cnt;
    yout[i] = sqrt(tmp);
  }
  return yout;
}


template <class _Tp>
inline double EventVar<_Tp>::Stddev(void)
{
  double sum = 0.0;
  double mean, tmp, stddev;
  int len = this->size();
  int i;

  for(i=0;i<len; i++)
    sum += (double)(*this)[i];

  mean = sum / (double)len;
  sum  = 0.0;
  for(i=0;i<len; i++)
  {
    tmp = (double)(*this)[i] - mean;
    sum = sum + (tmp * tmp);
  }
  tmp = sum / (double)len;
  stddev = sqrt(tmp);
  return stddev;
}


#endif