 module HandleTimeSeries         
   

  use selectrealkind , only : kpr
  use Incapsulated , only : GetUnit,GetRows    
  use interpolation
  use timeconversions 
  use HandleTimeStrings
  use CuteStrings
  use statistics
  use Numbers
  use Search 
  use DigitalFilters
  use linearalgebra 
  use linearregression
  use GeodeticConversions
  use PolynomialRoots
  use SingularSpectrumAnalysis 
  use iso_fortran_env
  use ifport

  implicit none 

  type series
   real(kind=kpr),allocatable,dimension(:) :: x ! 
   real(kind=kpr),allocatable,dimension(:) :: y ! f(x) 
   real(kind=kpr),allocatable,dimension(:) :: model ! .. used in LS fit problems 
   real(kind=kpr),allocatable,dimension(:) :: sy ! formal errors related to y 
   contains 
     procedure,nopass :: Allocate
     procedure,pass :: ConvertTimeTag
     procedure,pass :: SingularSpectrumReconstruction
     procedure,pass :: FillTheGapsWithZeros 
     procedure,pass :: FillTheGapsWithTheReconstructedValues 
     procedure,pass :: GapDetector
     procedure,pass :: Sort 
     procedure,pass :: GetTheWRMS 
     procedure,pass :: GetTheLSQLinearTrend 
     procedure,pass :: NormalizeTheSeries 
     procedure,pass :: UnNormalizeTheSeries 
     procedure,pass :: FindInBetween
     procedure,pass :: PremoliRecursiveFiltering
     procedure,pass :: Remove
     procedure,pass :: write => writets
     procedure,pass :: read => readts
     procedure,pass :: fit => fitTimeSeries
     procedure,pass :: bin => binTimeSeries 
     procedure,pass :: butterworth => butterworthTimeSeries  
     procedure,pass :: diff => diffTimeSeries
     procedure,pass :: trim => trimTimeSeries
  end type series

  type, extends(series) :: PositionTimeSeries
    character(len=1) :: technique 
    character(len=:),allocatable :: id
    real(kind=kpr) :: lat,lon
    real(kind=kpr),dimension(:),allocatable :: time
    integer :: nobservations
    type(series) :: e,n,u
    type(series) :: cen,ceu,cnu
    real(kind=kpr),dimension(:,:),allocatable :: covariance
    contains 
      !procedure,nopass :: allocate => allocatePositionTimeSeries
      procedure,pass :: writePositionTimeSeries
      procedure,pass :: writePositionTimeSeriesHeader
      procedure,pass :: writeSREFPositionTimeSeriesHeader
      procedure,pass :: scale => scalePositionTimeSeries
      procedure,pass :: enu2xyz 
      procedure,pass :: xyz2enu
      procedure,pass :: binPositionTimeSeries 
      procedure,pass :: butterworthPositionTimeSeries
      procedure,pass :: fitPositionTimeSeries
      procedure,pass :: diffPositionTimeSeries
      procedure,pass :: trimPositionTimeSeries
      procedure,pass :: getStdev 
      procedure,pass :: stdev => stdevPositionTimeSeries
      procedure      :: add => addPositionTimeSeries
      generic        :: operator(+) => add
  end type PositionTimeSeries 

  private 

  character(len=80),dimension(3),parameter :: & 
    formatArray = ['(i1,1x,a1,15x)',& 
                   '(f24.12,6x,2(f24.12,2x))', &
                   '(i5,1x,a10,1x)']

  public :: & 
             PositionTimeSeries,& 
             allocatePositionTimeSeries,&
             findObservationsAtTime,&
             series, & 
             AllocateTimeSeries, &
             DeallocateTimeSeries, & 
             ReadTimeSeries, &
             SelectTimeSeries, &
             !WriteTimeSeries,&
             WriteResidualTimeSeries, &
             WriteModelandResidualsTimeSeries, &
             SelectCommonObs,&
             tch, & 
             GenTCH, & 
             DetectGaps, &
             Difference, & 
             FilterTime, &
             DeMean,&
             MakeTemporalMask, & 
             ApplyTemporalMask, & 
             PadTimeSeriesWithZeros, &
             WhereToInterpolate, &
             InterpolateTimeSeries, &
             PadWithZerosAndInterpolateTimeSeries,&
             readPositionTimeSeries,& 
             takeTheDifferences

  contains 

  function stdevPositionTimeSeries(self) result(stdeviations)

    implicit none 
    class(PositionTimeSeries),intent(in) :: self 
    real(kind=kpr),dimension(:),allocatable :: stdeviations 

    stdeviations = [ stdev(self.e.y),stdev(self.n.y),stdev(self.u.y) ] 

  end function stdevPositionTimeSeries

  subroutine getStdev(self,code,domes,filename) 

    implicit none 
    class(PositionTimeSeries),intent(in) :: self 
    character(len=*),intent(in) :: code,domes 
    character(len=*),intent(in) :: filename 

    !Local Variables 
    integer :: iError,iUnit 
    real(kind=kpr),dimension(:),allocatable :: stdeviations 

    stdeviations = [ stdev(self.e.y),stdev(self.n.y),stdev(self.u.y) ] 

    open(newunit=iUnit,&
         file=filename,&
         action='write',&
         iostat=iError) 
    if(iError/=0)stop' [getStdev] file not opened '

    write(*,*) ' ... Writing ',trim(filename)

    write(unit=iUnit,fmt=1) code,domes,self.lat,self.lon,&
                            stdeviations

    1 format(a4,2x,a9,2(2x,f12.5),3(2x,f12.3))

  end subroutine getStdev 

  function addPositionTimeSeries(a,b) result(output) 

    implicit none 
    class(PositionTimeSeries),intent(in) :: a,b
    class(PositionTimeSeries),allocatable :: output 

    class(series),allocatable :: sum 

    if(.not.allocated(a.time)) stop ' [addPositionTimeSeries] a.time not allocated ' 
    if(.not.allocated(b.time)) stop ' [addPositionTimeSeries] b.time not allocated ' 
    if((size(a.time)-size(b.time))/=0) stop ' [addPositionTimeSeries] a anb b have different sizes'

    allocate(output)

    output.e = addTimeSeries( a.e,b.e )
    output.n = addTimeSeries( a.n,b.n )
    output.u = addTimeSeries( a.u,b.u )

    output.time = a.time
    output.e.x  = a.e.x
    output.n.x  = a.n.x
    output.u.x  = a.u.x
    output.lat  = a.lat
    output.lon  = a.lon


  end function addPositionTimeSeries 

  function addTimeSeries(a,b) result(output) 

    implicit none 
    class(series),intent(in) :: a,b
    class(series),allocatable :: output 

    real(kind=kpr),dimension(:),allocatable :: vector 

    if(.not.allocated(a.x)) stop ' [addTimeSeries] a.x not allocated ' 
    if(.not.allocated(b.x)) stop ' [addTimeSeries] b.x not allocated ' 
    if((size(a.x)-size(b.x))/=0) stop ' [addTimeSeries] a anb b have different sizes'

    allocate(output) 
    output.x = a.x
    output.y = (a.y+b.y) 

  end function addTimeSeries 

  function trimPositionTimeSeries(self,before,after) result(output) 

    implicit none 

    class(PositionTimeSeries),intent(in) :: self 
    real(kind=kpr),intent(in),optional :: before,after 
    class(PositionTimeSeries),allocatable :: output 

    !Local Variables 
    real(kind=kpr) :: rAfter,rBefore 
    logical :: lDebug 

    lDebug=.false.

    if(.not.allocated(self.time))stop' [trimTimeSeries] self.time not allocated '

    select case (present(before))
     case (.true.)
        rBefore = before 
     case (.false.)
        rBefore = minval(self.time)
    end select 

    select case (present(after))
     case (.true.)
        rafter = after 
     case (.false.)
        rafter = minval(self.time)
    end select 

    allocate(output)
    output.e    = self.e.trim(before=rBefore,after=rAfter) 
    output.n    = self.n.trim(before=rBefore,after=rAfter)
    output.u    = self.u.trim(before=rBefore,after=rAfter)
    output.lat  = self.lat 
    output.lon  = self.lon 
    output.time = output.e.x 

    if(.not.allocated(output.e.x)) output.e.x = self.e.x 
    if(.not.allocated(output.n.x)) output.n.x = self.n.x 
    if(.not.allocated(output.u.x)) output.u.x = self.u.x 

  end function trimPositionTimeSeries

  function diffPositionTimeSeries(self,order) result(output)

    implicit none 

    class(PositionTimeSeries),intent(in) :: self 
    integer,intent(in) :: order 
    class(PositionTimeSeries),allocatable :: output 

    allocate(output)

    output.e = self.e.diff(order=order) 
    output.n = self.n.diff(order=order)
    output.u = self.u.diff(order=order)

    output.lat  = self.lat 
    output.lon  = self.lon 
    output.time = output.e.x 

    if(.not.allocated(output.e.x)) output.e.x = self.e.x 
    if(.not.allocated(output.n.x)) output.n.x = self.n.x 
    if(.not.allocated(output.u.x)) output.u.x = self.u.x 

  end function diffPositionTimeSeries 

  function trimTimeSeries(self,before,after) result(output) 

    implicit none 

    class(series),intent(in) :: self 
    real(kind=kpr),intent(in),optional :: before,after 
    class(series),allocatable :: output 

    !Local Variables 
    integer :: i,nElements 
    real(kind=kpr) :: rBefore,rAfter

    if(.not.allocated(self.x))stop' [trimTimeSeries] self.x not allocated '

    if(.not.(present(before)))then
        rBefore = minval(self.x)
    else 
        rBefore = before 
    endif 

    if(.not.(present(after)) )then 
        rAfter = maxval(self.x)
    else 
        rAfter = after 
    endif 

    allocate(output)
    nElements=0
    do i=1,size(self.x) 
       if( (self.x(i)>=rBefore).and.(self.x(i)<=rAfter) ) nElements=nElements+1
    enddo 
    allocate(output.x(nElements))
    allocate(output.y(nElements))
    if(allocated(self.sy)) allocate(output.sy(nElements))

    nElements=0
    do i=1,size(self.x) 
       if( (self.x(i)>=rBefore).and.(self.x(i)<=rAfter) ) then 
            nElements=nElements+1
            output.x(nElements) = self.x(i)
            output.y(nElements) = self.y(i) 
            if(allocated(self.sy)) output.sy(nElements) = self.sy(i) 
       endif 
    enddo 

  end function trimTimeSeries 

  function diffTimeSeries(self,order) result(output) 

    implicit none 

    class(series),intent(in) :: self 
    integer,intent(in) :: order 
    class(series),allocatable :: output 

    allocate(output)
    allocate(output.x(size(self.x)-int(order)))
    allocate(output.y(size(self.x)-int(order)))

    output.x=self.x(order+1:) 
    output.y=diff(a=self.y,n=order) 

  end function diffTimeSeries 

  function butterworthPositionTimeSeries(self,cutoffFrequency,samplingFrequency,order) result(output)

    implicit none 

    class(PositionTimeSeries),intent(in) :: self 
    real(kind=kpr),intent(in) :: cutoffFrequency,samplingFrequency
    integer,intent(in) :: order 
    class(PositionTimeSeries),allocatable :: output 

    allocate(output)

    output.e = self.e.butterworth(cutoffFrequency=cutoffFrequency,samplingFrequency=samplingFrequency,order=order)
    output.n = self.n.butterworth(cutoffFrequency=cutoffFrequency,samplingFrequency=samplingFrequency,order=order)
    output.u = self.u.butterworth(cutoffFrequency=cutoffFrequency,samplingFrequency=samplingFrequency,order=order)

    output.lat = self.lat 
    output.lon = self.lon 
    output.time = self.e.x 
    if(.not.allocated(output.e.x)) output.e.x = self.e.x 
    if(.not.allocated(output.n.x)) output.n.x = self.n.x 
    if(.not.allocated(output.u.x)) output.u.x = self.u.x 

  end function butterworthPositionTimeSeries 

  function butterworthTimeSeries(self,cutoffFrequency,samplingFrequency,order) result(output) 

    implicit none 

    class(series),intent(in) :: self 
    real(kind=kpr),intent(in) :: cutoffFrequency,samplingFrequency
    integer,intent(in) :: order 
    class(series),allocatable :: output 
    
    integer :: iAllocate
    real(kind=kpr),dimension(:),allocatable :: filtered

    if(.not.allocated(self.x))stop' [butterworthTimeSeries] self.x not allocated '
    if(.not.allocated(self.y))stop' [butterworthTimeSeries] self.y not allocated '

    !Output Allocation
    allocate(output,stat=iAllocate)
    if(iAllocate/=0)stop' [butterworthTimeSeries] output not allocated ' 
    !Let's copy self.x that contains t-time
    allocate(output.x,source=self.x)
    !Let's copy self.sy
    if(allocated(self.sy)) allocate(output.y,source=self.sy)

    !Filter self.y
    filtered=ButterLowPassFilter(&
                 x=self.y,&
                 cf=cutoffFrequency,&
                 fs=SamplingFrequency,&
                 order=order) 
    allocate(output.y,source=filtered)

  end function butterworthTimeSeries

  function binPositionTimeSeries(self,sampling,timeunit) result(output) 

    implicit none 

    class(PositionTimeSeries),intent(in) :: self 
    character(len=*),intent(in) :: sampling 
    character(len=*),intent(in) :: timeunit
    class(PositionTimeSeries),allocatable :: output 

    allocate(output)

    output.e = self.e.bin(sampling=sampling,timeunit=timeunit)
    output.n = self.n.bin(sampling=sampling,timeunit=timeunit)
    output.u = self.u.bin(sampling=sampling,timeunit=timeunit)

    output.lat  = self.lat 
    output.lon  = self.lon 
    output.time = output.e.x 

  end function binPositionTimeSeries 

  function binTimeSeries(self,sampling,timeunit) result(output) 

    implicit none 

    class(series),intent(in) :: self 
    character(len=*),intent(in) :: sampling 
    character(len=*),intent(in) :: timeunit
    class(series),allocatable :: output 

    !Local Variables 
    integer :: i,iCounter
    integer,dimension(:),allocatable :: inx 
    real(kind=kpr) :: meanValue
    real(kind=kpr),dimension(:),allocatable :: timetag,uniqueTimeTag
    logical :: lError 

    if (.not.allocated(self.x)) stop ' [binTimeSeries] self.x not allocated '
    allocate(output) 
    allocate(output.x(0))
    allocate(output.y(0))

    select case ( lower(sampling) ) 
     case ('weekly')

       select case ( lower(timeunit) ) 
        case ('mjd')
          !convert to gpsw
          allocate(timeTag,mold=self.x)
          timeTag=zero 
          do i=1,size(timeTag)
             timeTag(i)=real(fromMJDToGPSw(mjd=self.x(i)),kind=kpr) 
          enddo
          uniqueTimeTag=unique_real(timeTag)
          do i=1,size(uniqueTimeTag) 
             inx=indexMatch_real(source=timeTag,target=uniqueTimeTag(i))
             !inx=findIndex(source=timeTag,target=uniqueTimeTag(i))
             if(size(inx)>0)then
                meanValue=mean(self.x(inx)) 
                output.x=[output.x,meanValue]
                meanValue=mean(self.y(inx))
                output.y=[output.y,meanValue]  
             elseif(size(inx)==0)then
               write(*,*) ' ### Time Tag: ', uniqueTimeTag(i)
               stop' [binTimeSeries] indexMatch_real found no match '
             endif
          enddo
        case ('decy','decyr','yr')
          !convert to gpsw
          stop' [binTimeSeries] case decy - tbc '
        case default 
          stop' [binTimeSeries] select case timeunit - Option not Implemented '
       end select

     case ('daily')

       select case ( lower(timeunit) ) 
        case ('mjd') 
          !do nothing other than removing decimals
          timeTag=real(int(self.x),kind=kpr) 
          uniqueTimeTag=unique_real(timetag)
          inx=[0]
          do i=1,size(uniqueTimeTag)-1 
             !inx=indexMatch_real(source=timeTag,target=uniqueTimeTag(i))
             inx=binarySearch(source=timeTag,target=uniqueTimeTag(i),startIndex=inx(size(inx))+1)
             if(size(inx)>0)then
                meanValue=mean(self.x(inx)) 
                output.x=[output.x,meanValue]
                meanValue=mean(self.y(inx))
                output.y=[output.y,meanValue]  
             elseif(size(inx)==0)then
               write(*,*) ' ### Time Tag: ', uniqueTimeTag(i)
               stop' [binTimeSeries] indexMatch_real found no match '
             endif
          enddo
        case ('decy','decyr','yr')
          !convert to mjd
          allocate(timeTag,mold=self.x)
          timeTag=zero 
          do i=1,size(timeTag)
             call yr2mjd(yr=self.x(i),mjd=timeTag(i),err=lError) 
          enddo 
          stop' [binTimeSeries] case decy - tbc '
        case default 
          stop'[binTimeSeries] select case timeunit - Option not Implemented '
       end select 

     case default
         stop'[binTimeSeries] select case sampling - Option not Implemented '

    end select 

    if(size(output.x)==0)stop' [binTimeSeries] output.x is empty ' 
    if(size(output.y)==0)stop' [binTimeSeries] output.y is empty ' 
    if(size(output.y)/=size(output.x))stop' [binTimeSeries] ODD: output.x and output.y have different dimensions ' 

  end function binTimeSeries

  function xyz2enu(self) result(output) 

    implicit none 

    class(positionTimeSeries),intent(in) :: self 
    class(positionTimeSeries),allocatable :: output 

    integer :: i,nElements 
    real(kind=kpr),dimension(:),allocatable :: vector

    if( (self.lat==zero).and.(self.lon==zero) )stop' [xyz2enu] self.lat and self.lon are undefined ' 
    nElements = size(self.time)

    allocate(output) 
    output.lat=self.lat 
    output.lon=self.lon
    allocate(output.time,source=self.time)
    allocate(output.e.x,source=self.time)
    allocate(output.n.x,source=self.time)
    allocate(output.u.x,source=self.time)
    allocate(output.e.y,source=self.e.y)
    allocate(output.n.y,source=self.e.y)
    allocate(output.u.y,source=self.e.y)
    output.e.y=zero 
    output.n.y=zero
    output.u.y=zero

    do i=1,nElements 
       vector = dx2dh( lat  = self.lat,&
                       lon  = self.lon,&
                       dx   = self.e.y(i),&
                       dy   = self.n.y(i),& 
                       dz   = self.u.y(i),&
                       unit = 'deg')
       output.e.y(i) = vector(1) ! -e-
       output.n.y(i) = vector(2) ! -n-
       output.u.y(i) = vector(3) ! -u-
    enddo

  end function xyz2enu

  function enu2xyz(self) result(output) 

    class(positionTimeSeries),intent(in) :: self 
    class(positionTimeSeries),allocatable :: output 

    integer :: i,nElements,iAllocate 
    real(kind=kpr),dimension(:),allocatable :: vector

    if( (self.lat==zero).and.(self.lon==zero) )stop' [enu2xyz] self.lat and self.lon are undefined ' 
    nElements = size(self.time)

    allocate(output,stat=iAllocate)
    if(iAllocate/=0) stop '[enu2xyz] output Allocation Error '

    output.lat=self.lat 
    output.lon=self.lon
    allocate(output.time,source=self.time)
    allocate(output.e.x,source=self.time)
    allocate(output.n.x,source=self.time)
    allocate(output.u.x,source=self.time)
    allocate(output.e.y,source=self.e.y)
    allocate(output.n.y,source=self.e.y)
    allocate(output.u.y,source=self.e.y)
    output.e.y=zero 
    output.n.y=zero
    output.u.y=zero

    do i=1,nElements 
       vector = dh2dx( lat  = self.lat,&
                       lon  = self.lon,&
                       de   = self.e.y(i),&
                       dn   = self.n.y(i),& 
                       du   = self.u.y(i),&
                       unit = 'deg')
       output.e.y(i) = vector(1) ! -x-
       output.n.y(i) = vector(2) ! -y-
       output.u.y(i) = vector(3) ! -z-
    enddo

  end function enu2xyz

  subroutine allocatePositionTimeSeries(self,nel)

    implicit none 
    class(positionTimeSeries),allocatable,intent(inout) :: self
    integer,intent(in) :: nel

    integer :: iAllocate
    
    if(allocated(self)) deallocate(self)
    allocate(self,stat=iAllocate)
    if(iAllocate/=0)stop' [allocatePositionTimeSeries] self allocation error '

    allocate(self.n.x(nel),stat=iAllocate)
    if(iAllocate/=0)stop' [allocatePositionTimeSeries] self.n.x allocation error '
    allocate(self.n.y(nel),stat=iAllocate)
    if(iAllocate/=0)stop' [allocatePositionTimeSeries] self.n.y allocation error '
    allocate(self.n.sy(nel),stat=iAllocate)
    if(iAllocate/=0)stop' [allocatePositionTimeSeries] self.n.sy allocation error '
    allocate(self.e.x(nel),stat=iAllocate)
    if(iAllocate/=0)stop' [allocatePositionTimeSeries] self.e.x allocation error '
    allocate(self.e.y(nel),stat=iAllocate)
    if(iAllocate/=0)stop' [allocatePositionTimeSeries] self.e.y allocation error '
    allocate(self.e.sy(nel),stat=iAllocate)
    if(iAllocate/=0)stop' [allocatePositionTimeSeries] self.e.sy allocation error '
    allocate(self.u.x(nel),stat=iAllocate)
    if(iAllocate/=0)stop' [allocatePositionTimeSeries] self.u.x allocation error '
    allocate(self.u.y(nel),stat=iAllocate)
    if(iAllocate/=0)stop' [allocatePositionTimeSeries] self.u.y allocation error '
    allocate(self.u.sy(nel),stat=iAllocate)
    if(iAllocate/=0)stop' [allocatePositionTimeSeries] self.u.sy allocation error '

  end subroutine allocatePositionTimeSeries 

  function takeTheDifferences(a,b,timeUnit) result(output) 

    implicit none 
    class(positionTimeSeries),intent(inout) :: a,b 
    character(len=*),intent(in) :: timeUnit
    class(positionTimeSeries),allocatable :: output 

    !Local Variables
    integer :: i,j,k
    integer :: iAllocate,iCounter
    integer,dimension(:),allocatable :: iVectorA,iVectorB
    integer :: nElementsInA,nElementsInB
    real(kind=kpr) :: timeMJD
    real(kind=kpr) :: timeTarget
    real(kind=kpr),dimension(:),allocatable :: formalErrors
    real(kind=kpr),dimension(:),allocatable :: timeCopyA,timeCopyB
    character(len=:),allocatable :: fromTo
    logical :: lDebug,lError,lCovariance

    lCovariance=.true.
    lDebug=.false. 
    if(lDebug) write(*,*) ' [takeTheDifferences] debug mode activated '

    ! output allocation 
    allocate(output,stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] output not allocated '

    ! Save original time tags 
    allocate(timeCopyA,source=a.e.x,stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] timeCopyA - Allocation Error ' 
    allocate(timeCopyB,source=b.e.x,stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] timeCopyB - Allocation Error ' 

    !Time Conversions
    fromTo=timeUnit//'mjd'
    nElementsInA=size(a.e.x)
    nElementsInB=size(b.e.x)
    if(lDebug)write(*,*) ' [takeTheDifferences] nElementsInA = ',nElementsInA
    if(lDebug)write(*,*) ' [takeTheDifferences] nElementsInB = ',nElementsInB
    select case (lower(fromTo))
     case ('decymjd','yrmjd')
       do i=1,nElementsinA
          call yr2mjd(a.e.x(i),& 
                      timeMJD,& 
                      lError)
          timeMJD=real(int(timeMJD),kind=kpr)  
          a.e.x(i)=timeMJD
          if(lDebug)write(*,*)a.e.x(i) 
          a.e.x(i)=timeMJD
          a.n.x(i)=timeMJD
          a.u.x(i)=timeMJD
          a.cen.x(i)=timeMJD
          a.ceu.x(i)=timeMJD
          a.cnu.x(i)=timeMJD
       enddo 
       do i=1,nElementsinB
          call yr2mjd(b.e.x(i),& 
                      timeMJD,& 
                      lError)
          timeMJD=real(int(timeMJD),kind=kpr)  
          b.e.x(i)=timeMJD
          if(lDebug)write(*,*)b.e.x(i) 
          b.n.x(i)=timeMJD
          b.u.x(i)=timeMJD
          b.cen.x(i)=timeMJD
          b.ceu.x(i)=timeMJD
          b.cnu.x(i)=timeMJD
       enddo 
     case ('mjdmjd')
       !do nothing
     case default 
       stop' [TakeDifference] select case option undefined '
    end select
    !Find Common Observations
    iCounter=0
    loopOverA: do i=1,nElementsInA
        timeTarget=a.e.x(i) 
        loopOverB: do j=1,nElementsInB
          if(timeTarget.eq.b.e.x(j)) iCounter=iCounter+1
        enddo loopOverB  
    enddo loopOverA
    if(lDebug) write(*,*)' iCounter = ',iCounter  
    if(iCounter.eq.0) stop ' [takeTheDifferences] a and b have no common observations ' 

    allocate(iVectorA(iCounter),stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] iVectorA not Allocated' 
    iVectorA=0
    allocate(iVectorB(iCounter),stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] iVectorB not Allocated' 
    iVectorB=0
    !Fill in the iVectors      
    iCounter=0
    loopOverA_2: do i=1,nElementsInA
        timeTarget=a.e.x(i) 
        loopOverB_2: do j=1,nElementsInB
          if(timeTarget.eq.b.e.x(j))then
            iCounter=iCounter+1
            iVectorA(iCounter)=i 
            iVectorB(iCounter)=j 
          endif 
        enddo loopOverB_2  
    enddo loopOverA_2
    if(lDebug)then 
      do i=1,size(iVectorA) 
         write(*,*), ' i = ',i,' iVectorA ', iVectorA(i),' iVectorB = ',iVectorB(i) 
      enddo
    endif
    !Allocate Time Series
    allocate(output.e.x,source=a.e.x(iVectorA),stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] output.e.x allocation error '
    allocate(output.e.y,source=a.e.y(iVectorA)-b.e.y(iVectorB),stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] output.e.y allocation error '
    allocate(output.n.x,source=a.n.x(iVectorA),stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] output.n.x allocation error '
    allocate(output.n.y,source=a.n.y(iVectorA)-b.n.y(iVectorB),stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] output.n.y allocation error '
    allocate(output.u.x,source=a.u.x(iVectorA),stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] output.u.x allocation error '
    allocate(output.u.y,source=a.u.y(iVectorA)-b.u.y(iVectorB),stat=iAllocate)
    if(iAllocate/=0)stop' [takeTheDifferences] output.n.y allocation error '
    if(lCovariance)then  
      allocate(formalErrors(size(output.e.x)),stat=iAllocate)
      if(iAllocate/=0)stop' [takeTheDifferences] formalErrors - allocation error '
      formalErrors=one 
      allocate(output.e.sy,source=formalErrors,stat=iAllocate)
      if(iAllocate/=0)stop' [takeTheDifferences] output.e.sy allocation error '
      allocate(output.n.sy,source=formalErrors,stat=iAllocate)
      if(iAllocate/=0)stop' [takeTheDifferences] output.e.ny allocation error '
      allocate(output.u.sy,source=formalErrors,stat=iAllocate)
      if(iAllocate/=0)stop' [takeTheDifferences] output.e.ny allocation error '
      allocate(output.cen.y,mold=output.e.x,stat=iAllocate)
      if(iAllocate/=0)stop' [takeTheDifferences] output.cen.y allocation error '
      output.cen.y=zero
      allocate(output.ceu.x,source=a.e.x(iVectorA),stat=iAllocate)
      if(iAllocate/=0)stop' [takeTheDifferences] output.ceu.x allocation error '
      allocate(output.ceu.y,mold=output.e.x(iVectorA),stat=iAllocate)
      if(iAllocate/=0)stop' [takeTheDifferences] output.ceu.y allocation error '
      output.ceu.y=zero
      allocate(output.cnu.x,source=a.e.x(iVectorA),stat=iAllocate)
      if(iAllocate/=0)stop' [takeTheDifferences] output.cnu.x allocation error '
      allocate(output.cnu.y,mold=output.e.x(iVectorA),stat=iAllocate)
      if(iAllocate/=0)stop' [takeTheDifferences] output.cnu.y allocation error '
      output.cnu.y=zero
    endif

    !Transform back to decy, if need be
    select case (lower(fromTo))
     case ('decymjd','yrmjd')
          a.e.x=timeCopyA
          a.n.x=timeCopyA
          a.u.x=timeCopyA
          a.cen.x=timeCopyA
          a.ceu.x=timeCopyA
          a.cnu.x=timeCopyA
          b.e.x=timeCopyB
          b.n.x=timeCopyB
          b.u.x=timeCopyB
          b.cen.x=timeCopyB
          b.ceu.x=timeCopyB
          b.cnu.x=timeCopyB
          output.e.x=timeCopyA(iVectorA) 
          output.n.x=timeCopyA(iVectorA) 
          output.u.x=timeCopyA(iVectorA) 
          output.cen.x=timeCopyA(iVectorA) 
          output.ceu.x=timeCopyA(iVectorA) 
          output.cnu.x=timeCopyA(iVectorA) 
    end select

    !Construct the metadata in the class output 
    output.lon=a.lon 
    output.lat=a.lat  
    output.nObservations=size(iVectorA)

  end function takeTheDifferences

  function findObservationsAtTime(input,time,unit) result(output) 

    implicit none 
    class(positionTimeSeries),dimension(:),intent(inout) :: input
    real(kind=kpr),intent(in) :: time
    real(kind=kpr),dimension(:,:),allocatable :: output
    character(len=*),intent(in) :: unit

    !Local Variables 
    integer :: i,j,k
    integer :: nSites 
    integer :: iAllocate
    integer :: iCounter
    real(kind=kpr) :: targetTime
    real(kind=kpr) :: tempMJD,tempDecy
    logical :: lDebug,lError
    logical :: lConversion 

    lDebug=.false.
    lConversion=.true.

    select case (lower(unit)) 
     case ('mjd')
       targetTime=time
     case ('decy')
       call yr2mjd(time,targetTime,lError)
     case default 
       stop' [findObservationsAtTime] select case unit - option undefined '
    end select 

    if(lDebug) write(*,*)' [findObservationsAtTime] lDebug Activated '

    if(lDebug)then 
      write(*,*)' target Time input  = ',time 
      write(*,*)' target Time (mjd)  = ',targetTime 
    endif

    nSites=size(input)
    if(lDebug) write(*,*)' nSites = ', nSites
    if(lConversion)then 
      do i=1,nSites
       do j=1,size(input(i).e.x)
          call yr2mjd(input(i).e.x(j),tempMJD,lError)
          input(i).e.x(j)=tempMJD
       enddo
      enddo  
    endif 
    !Search for targetTime 
    iCounter=0
    do i=1,nSites
     do j=1,size(input(i).e.x)
        if(input(i).e.x(j).eq.targetTime)then
          iCounter=iCounter+1
        endif 
     enddo
    enddo  
    if(lDebug)write(*,*)' iCounter = ',iCounter
    
    allocate(output(iCounter,11),stat=iAllocate)
    if(iAllocate/=0)stop' [findObservationsAtTime] output - allocation error '

    !Fill in output 
    iCounter=0
    do i=1,nSites
     do j=1,size(input(i).e.x)
        if(input(i).e.x(j).eq.targetTime)then
          iCounter=iCounter+1
          output(iCounter,1)=input(i).lat
          output(iCounter,2)=input(i).lon
          output(iCounter,3)=input(i).e.y(j)
          output(iCounter,4)=input(i).n.y(j)
          output(iCounter,5)=input(i).u.y(j)
          output(iCounter,6)=input(i).e.sy(j)
          output(iCounter,7)=input(i).n.sy(j)
          output(iCounter,8)=input(i).u.sy(j)
          output(iCounter,9)=input(i).cen.y(j)
          output(iCounter,10)=input(i).ceu.y(j)
          output(iCounter,11)=input(i).cnu.y(j)
        endif 
     enddo
    enddo  
    if(lDebug)then 
      write(*,*) ' Print Output '
      do i=1,size(output,dim=1)
         write(*,*) (output(i,j),j=1,5)
      enddo 
    endif
    if(lDebug)then 
      write(*,*) ' Print Output -- Variances '
      do i=1,size(output,dim=1)
         write(*,*) (output(i,j),j=6,8)
      enddo 
    endif
    if(lDebug)then 
      write(*,*) ' Print Output -- Covariances '
      do i=1,size(output,dim=1)
         write(*,*) (output(i,j),j=9,11)
      enddo 
    endif

    ! Very Last Step - we Reconvert time tag  
    if(lConversion)then 
      do i=1,nSites
       do j=1,size(input(i).e.x)
          call mjd2yr(input(i).e.x(j),tempDecy,lError)
          input(i).e.x(j)=tempDecy
       enddo
      enddo  
    endif 

  end function findObservationsAtTime

  subroutine readPositionTimeSeries(output,fname,timeunit) 

    implicit none 
    class(positionTimeSeries),intent(out),allocatable :: output
    character(len=*),intent(in) :: fname 
    character(len=*),intent(in),optional :: timeunit 

    !Local Variables
    integer :: i,ifile
    integer :: nElements,yr
    real(kind=kpr) :: timeMJD,doy
    logical :: lError
    logical :: lDebug 
    class(positionTimeSeries),allocatable :: elasticDisplacements 
    character(len=:),allocatable :: units,tunit

    lDebug=.false.
    if(lDebug) write(*,*)' [readPositionTimeSeries] debugging mode activated '  

    !Handle Optional
    select case (present(timeunit))
     case (.true.)
      tunit=timeunit
     case (.false.) 
      tunit='decy'
    end select 
    if(lDebug) write(*,*) ' tUnit = ' , tunit 

    call readENU(fname=fname,output=elasticDisplacements,covariance=.true.)
    call move_alloc(from=elasticDisplacements,to=output)
    nElements=size(output.e.x)

    select case (lower(tUnit)) 
     case ('mjd')
       do i=1,nElements 
        call yr2mjd(output.e.x(i),timeMJD,lError)
        output.e.x(i)=timeMJD
        output.n.x(i)=timeMJD
        output.u.x(i)=timeMJD
        output.cen.x(i)=timeMJD
        output.ceu.x(i)=timeMJD
        output.cnu.x(i)=timeMJD
       enddo 
     case ('doy')
       do i=1,nElements 
        call yr2jd(output.e.x(i),timeMJD,lError)
        call jd2doy(timeMJD,yr,doy,lError)
        output.e.x(i)=doy
        output.n.x(i)=doy
        output.u.x(i)=doy
        output.cen.x(i)=doy
        output.ceu.x(i)=doy
        output.cnu.x(i)=doy
       enddo 
     case ('decy')
       !Does Nothing
     case default 
       stop' [ReadPositionTimeSeries] select case option undefined ' 
    end select 

    if(lDebug) write(*,*) ' Number of Elements in positionTimeSeries = ' , nElements
    if(lDebug) write(*,*) ' Time Tags = ' , output.e.x 
    if(lDebug) write(*,*) ' e.y  = ' , output.e.y
    if(lDebug) write(*,*) ' e.sy = ' , output.e.sy
    if(lDebug) write(*,*) ' n.y  = ' , output.n.y
    if(lDebug) write(*,*) ' n.sy = ' , output.n.sy
    if(lDebug) write(*,*) ' u.y  = ' , output.u.y
    if(lDebug) write(*,*) ' u.sy = ' , output.u.sy

  end subroutine readPositionTimeSeries 

  subroutine readENU(fname,output,covariance) 

    implicit none 
    character(len=*),intent(in) :: fname
    class(PositionTimeSeries),allocatable,intent(out) :: output
    logical,intent(in),optional :: covariance

    !Local Variables 
    integer :: i 
    integer :: iUnit,iError,iAllocate
    integer :: nRows,nElements
    real(kind=kpr),dimension(:),allocatable :: formalErrors
    real(kind=kpr),dimension(:),allocatable :: time,e,n,u
    logical :: lDebug,lCovariance

    lDebug = .false. 
    
    !Handle Optional Argument
    select case (present(covariance))
     case (.true.)
       lCovariance=covariance
     case (.false.)
       lCovariance=.false.
    end select 

    !Allocate Output 
    allocate(output,stat=iAllocate)
    if(iAllocate/=0)stop' [readENU] Allocation Error ' 

    open(newunit=iUnit,&
         file=fname,&
         action='read',&
         status='unknown',& 
         iostat=iError,&
         access='sequential'&
         ) 
    if(iError/=0)stop' [readENU] File Not Found '

    if(lDebug) write(*,*) ' [readENU] fname = ' , trim(fname) 
    call GetRows(fname=fname,nRow=nRows) 
    if(lDebug) write(*,*) ' [readENU] nRows = ' , nRows 
    nElements=nRows-1
    allocate(time(nElements),stat=iAllocate) 
    if(iAllocate/=0)stop' [readENU] Time Allocation Error ' 
    allocate(e,n,u,source=time,stat=iAllocate)
    if(iAllocate/=0)stop' [readENU] e,n,u Allocation Error ' 
    !Read Lat and Lon from First Line 
    read(unit=iUnit,fmt='(6x,f12.4,10x,f12.5)') output.lon,output.lat 
    if(lDebug)then 
      write(*,*)' output.lon =  ' , output.lon 
      write(*,*)' output.lat =  ' , output.lat 
    endif 
    output.nobservations=nElements
    do i=1,nElements 
       read(unit=iUnit,fmt='(f12.6,3(1x,f30.25),f18.4)')time(i),e(i),n(i),u(i)
       if(time(i).eq.zero)then 
         print*,' found zero value in time(i) i = ',i,trim(fname)
         stop
       endif
    enddo
    !Allocate Time Series
    allocate(output.e.x,source=time,stat=iAllocate)
    if(iAllocate/=0)stop' [readENU] output.e.x allocation error '
    allocate(output.e.y,source=e,stat=iAllocate)
    if(iAllocate/=0)stop' [readENU] output.e.y allocation error '
    allocate(output.n.x,source=time,stat=iAllocate)
    if(iAllocate/=0)stop' [readENU] output.n.x allocation error '
    allocate(output.n.y,source=n,stat=iAllocate)
    if(iAllocate/=0)stop' [readENU] output.n.y allocation error '
    allocate(output.u.x,source=time,stat=iAllocate)
    if(iAllocate/=0)stop' [readENU] output.u.x allocation error '
    allocate(output.u.y,source=u,stat=iAllocate)
    if(iAllocate/=0)stop' [readENU] output.u.y allocation error '
    if(lCovariance)then  
      allocate(formalErrors(size(output.e.x)),stat=iAllocate)
      if(iAllocate/=0)stop' [readENU] formalErrors - allocation error '
      formalErrors=one 
      allocate(output.e.sy,source=formalErrors,stat=iAllocate)
      if(iAllocate/=0)stop' [readENU] output.e.sy allocation error '
      allocate(output.n.sy,source=formalErrors,stat=iAllocate)
      if(iAllocate/=0)stop' [readENU] output.e.ny allocation error '
      allocate(output.u.sy,source=formalErrors,stat=iAllocate)
      if(iAllocate/=0)stop' [readENU] output.e.ny allocation error '
      allocate(output.cen.x,source=time,stat=iAllocate)
      if(iAllocate/=0)stop' [readENU] output.cen.x allocation error '
      allocate(output.cen.y,mold=output.e.x,stat=iAllocate)
      if(iAllocate/=0)stop' [readENU] output.cen.y allocation error '
      output.cen.y=zero
      allocate(output.ceu.x,source=time,stat=iAllocate)
      if(iAllocate/=0)stop' [readENU] output.ceu.x allocation error '
      allocate(output.ceu.y,mold=output.e.x,stat=iAllocate)
      if(iAllocate/=0)stop' [readENU] output.ceu.y allocation error '
      output.ceu.y=zero
      allocate(output.cnu.x,source=time,stat=iAllocate)
      if(iAllocate/=0)stop' [readENU] output.cnu.x allocation error '
      allocate(output.cnu.y,mold=output.e.x,stat=iAllocate)
      if(iAllocate/=0)stop' [readENU] output.cnu.y allocation error '
      output.cnu.y=zero
    endif 
    close(unit=iUnit)

  end subroutine readENU

  subroutine scalePositionTimeSeries(self,scalingFactor) 

    implicit none 
    class(positionTimeSeries),intent(inout) :: self
    real(kind=kpr),intent(in) :: scalingFactor

    !Local Variables 
    self.e.y=self.e.y*scalingFactor
    self.e.sy=self.e.sy*scalingFactor
    self.n.y=self.n.y*scalingFactor
    self.n.sy=self.n.sy*scalingFactor
    self.u.y=self.u.y*scalingFactor
    self.u.sy=self.u.sy*scalingFactor

  end subroutine scalePositionTimeSeries

  subroutine writePositionTimeSeries(self,fname,timeunit,headertype,covariance,scalingfactor,fmt)

    implicit none 
    class(positionTimeSeries),intent(in) :: self
    character(len=*),intent(in) :: fname 
    character(len=*),intent(in),optional :: timeunit,headertype 
    logical,intent(in),optional :: covariance
    real(kind=kpr),intent(in),optional :: scalingFactor
    character(len=*),intent(in),optional :: fmt  

    !Local Variables
    integer :: i,ifile
    real(kind=kpr) :: temp,ConversionFactor,scaleFactor 
    logical :: lError,lCovariance
    character(len=:),allocatable :: units,tunit,tHeader,format

    !Handle Optional
    select case (present(timeunit))
     case (.true.)
      tunit=timeunit
     case (.false.) 
      tunit='decy'
    end select 
    select case (present(headertype))
     case (.true.)
      tHeader=headertype
     case (.false.) 
      tHeader='verbose'
    end select 
    select case (present(covariance))
     case (.true.)
       lCovariance=covariance
     case (.false.) 
       lCovariance=.true.
    end select
    select case (present(scalingfactor))
     case (.true.)
       scaleFactor=scalingFactor
     case (.false.)
       scaleFactor=one
    end select   
    select case (present(fmt))
     case (.false.)
       format='(f12.6,3(1x,f30.24),3(3x,e8.3),3(3x,e8.3))'
     case (.true.)
       format=fmt
    end select

    open(newunit=ifile,file=fname,action='write',&
         access='sequential',status='unknown')

    !Write Header 
    select case (lower(tHeader))
      case ('verbose') 
        call self.writePositionTimeSeriesHeader(iUnit=ifile)
      case ('sref')
        call self.writeSREFPositionTimeSeriesHeader(iUnit=ifile)
      case ('short') 
        call writeShortHeader(iUnit=ifile,& 
                                  lat=self.lat,& 
                                  lon=self.lon)
      case default
        stop' [writePositionTimeSeriesHeader] select case tHeader - option undefined '
    end select 

    do i=1,size(self.e.x)  
       select case (lower(tunit))
        case ('mjd')
          call mjd2yr(self.e.x(i),temp,lError) 
        case ('decy')
          temp=self.e.x(i) 
       end select 
       select case (lCovariance) 
        case (.false.) 
              write(unit=ifile,fmt=format) temp,self.e.y(i)*scaleFactor,&
                                           self.n.y(i)*scaleFactor,&
                                           self.u.y(i)*scaleFactor
        case (.true.) 
              write(unit=ifile,fmt=2) temp,self.e.y(i)*scaleFactor,&
                                           self.n.y(i)*scaleFactor,&
                                           self.u.y(i)*scaleFactor,&
                                           self.e.sy(i)*scaleFactor,& 
                                           self.n.sy(i)*scaleFactor,&
                                           self.u.sy(i)*scaleFactor,&
                                           self.cen.y(i),&
                                           self.ceu.y(i),&
                                           self.cnu.y(i)
       end select
    end do 

    close(unit=ifile) 

  1 format(f12.4,2x,e23.15)  
  2 format(f12.6,3(1x,f30.24),3(3x,e8.3),3(3x,e8.3))  

  end subroutine writePositionTimeSeries

  subroutine WriteShortHeader(iUnit,lat,lon)
 
    implicit none 
    integer,intent(in) :: iUnit
    real(kind=kpr),intent(in) :: lat,lon
 
    write(unit=iUnit,fmt=1) lon,lat  
    1 format (' rlon=',f12.4,5x,'rlat=',f12.5 )
 
  end subroutine WriteShortHeader  

  function fitPositionTimeSeries(self,model,kind,xshift,timeunit,logfile) result(output) 

   implicit none 
   class(PositionTimeSeries),intent(in) :: self
   character(len=*),intent(in) :: model 
   character(len=*),intent(in) :: kind 
   real(kind=kpr),intent(in),optional :: xshift 
   character(len=*),intent(in),optional :: timeunit 
   character(len=*),intent(in),optional :: logfile
   class(PositionTimeSeries),allocatable :: output 

    allocate(output) 
    allocate(output.time,source=self.e.x)
    output.lat = self.lat
    output.lon = self.lon

    output.e = self.e.fit(model=model,kind=kind,xshift=xshift,timeunit=timeunit,logfile=logfile//'.e')
    output.n = self.n.fit(model=model,kind=kind,xshift=xshift,timeunit=timeunit,logfile=logfile//'.n')
    output.u = self.u.fit(model=model,kind=kind,xshift=xshift,timeunit=timeunit,logfile=logfile//'.u')

  end function fitPositionTimeSeries 

  function fitTimeSeries(self,model,kind,xshift,timeunit,logfile) result(output)

   implicit none 
   class(series),intent(in) :: self
   character(len=*),intent(in) :: model 
   character(len=*),intent(in) :: kind 
   real(kind=kpr),intent(in),optional :: xshift 
   character(len=*),intent(in),optional :: timeunit 
   character(len=*),intent(in),optional :: logfile
   class(series),allocatable :: output 

   !Local Variables 
   integer :: i,j,iAllocate
   real(kind=kpr) :: shift
   real(kind=kpr),dimension(:),allocatable :: x,y,sy
   character(len=:),allocatable :: time,log
   logical :: lDebug
   type(NormalSystem) :: N

   lDebug=.false. 
   if(lDebug) write(*,*)' [fitTimeSeries] debugging mode activated '

   if(allocated(self.x).eqv..false.)stop' [fitTimeSeries] self.x Not Allocated '

   !Handle Optional 
   select case (present(xshift))
    case (.true.)
      shift=xshift
    case (.false.)
      shift=mean(self.x)
   end select     
   select case (present(timeunit))
    case (.true.)
      time=timeunit
    case (.false.)
      time='mjd'
   end select     
   select case (present(logfile))
    case (.true.)
      log=logfile
    case (.false.)
      log='./log.lsq'
   end select     

   !Copy and Allocate per Component 
   allocate(x,source=self.x,stat=iAllocate)
   if(iAllocate/=0)stop' [fitTimeSeries] x allocation error'
   allocate(y,source=self.y,stat=iAllocate)
   if(iAllocate/=0)stop' [fitTimeSeries] y allocation error'
   if(allocated(self.sy)) allocate(sy,source=self.sy,stat=iAllocate)
   if(iAllocate/=0)stop' [fitTimeSeries] sy allocation error'

   select case (lower(model)) 
    case ('mean')
      allocate(N.obs,source=y)
      allocate(N.res,source=(y-mean(y)))
      allocate(N.sol(1),stat=iAllocate)
      if(iAllocate/=0)stop' [fitTimeSeries] N.sol(1) not allocated '
      N.sol=mean(y)
      call N.WriteSolution(filename=log,model='mean')
    case default  
      N=BuildAndSolveNormalEquations(x=x,y=y,xshift=shift,model=model,timeunit=time,logfile=log)
   end select
   !Allocate output
   allocate(output)
   allocate(output.x,source=self.x)
   allocate(output.y,source=N.res)
   if(lDebug)then 
    do i=1,size(x)
       write(*,*) x(i),y(i),sy(i),output.y(i) 
    enddo
   endif
   if(allocated(x)) deallocate(x)
   if(allocated(y)) deallocate(y)
   if(allocated(sy)) deallocate(sy)

  end function fitTimeSeries

  subroutine writets(self,fname,kind,timeunit,fmt)

    implicit none 
    class(series),intent(in) :: self
    character(len=*),intent(in) :: fname 
    character(len=*),intent(in) :: kind 
    character(len=*),intent(in),optional :: timeunit 
    character(len=*),intent(in),optional :: fmt  

    !Local Variables
    integer :: i,ifile
    real(kind=kpr) :: temp,ConversionFactor 
    logical :: lError
    character(len=:),allocatable :: units,tunit

    !Handle Optional
    select case (present(timeunit))
     case (.true.)
      tunit=timeunit
     case (.false.) 
      tunit='decy'
    end select 

    if(allocated(self.x).eqv..false.)then 
      stop ' [writets] self.x not allocated '
    endif 
    if(allocated(self.y).eqv..false.)then 
      stop' [writets] self.y not Allocated '
    endif 

    open(newunit=ifile,file=fname,action='write',&
         access='sequential',status='unknown')

    !Write Header 
    call WriteTimeSeriesHeader(iUnit=ifile,unit=' ',nElements=size(self.x))
    do i=1,size(self.x)  
       select case (lower(tunit))
        case ('mjd')
          call mjd2yr(self.x(i),temp,lError) 
        case ('decy')
          temp=self.x(i) 
       end select 
       select case (kind) 
        case ('xy')
          write(unit=ifile,fmt=1) temp,self.y(i)
          !write(unit=ifile,fmt=1) temp,self.x(i),self.y(i)
        case ('xye')
          select case (present(fmt))
           case (.false.)
            write(unit=ifile,fmt=2) self.x(i),self.y(i),self.sy(i)
           case (.true.)
            write(unit=ifile,fmt=fmt) self.x(i),self.y(i),self.sy(i)
          end select
        case default 
          stop ' [writets] kind undefined '
       end select
    end do 
    close(unit=ifile)

  1 format(f12.4,2x,e23.15)  
  2 format(f12.4,2(3x,e23.15))  
  !1 format(f12.4,2(3x,e23.15),f18.6)  
  !2 format(f12.4,3(3x,e23.15),f18.6)  

  end subroutine writets  

  function Remove(self,model,xshift,timeunit,weighting,logfile,texfile) result(output) 

    implicit none 
    class(series),intent(inout) :: self
    character(len=*),intent(in) :: model 
    real(kind=kpr),intent(in),optional :: xshift 
    character(len=*),intent(in),optional :: timeunit,logfile,texfile 
    logical,optional :: weighting

    class(series), allocatable :: output 

   !Local Variables 
   type(NormalSystem) :: N
   real(kind=kpr) :: shift
   real(kind=kpr),dimension(:),allocatable :: x,y,sy
   character(len=:),allocatable :: time,lfile,tfile
   logical :: lWeighting

   !Sanity Check 
   if(allocated(self.x).eqv..false.)stop' [Remove] self.x Not Allocated '
   if(allocated(self.y).eqv..false.)stop' [Remove] self.y Not Allocated '
   if(.not. allocated(self.sy))then 
     allocate(self.sy,mold=self.x)
     self.sy=one 
   endif  

   !Handle Optional 
   select case (present(xshift))
    case (.true.)
      shift=xshift
    case (.false.)
      shift=mean(self.x)
   end select     
   select case (present(timeunit))
    case (.true.)
      time=timeunit
    case (.false.)
      time='mjd'
   end select     

   !Default is unweighted fit
   select case (present(weighting))
    case (.true.)
      lWeighting=weighting
    case (.false.)
      lWeighting=.false.
   end select     
   select case (present(logfile))
    case (.true.)
      lfile=logfile
    case (.false.)
      lfile='timeseries.lsq'
   end select   
   select case (present(texfile))
    case (.true.)
      tfile=texfile
    case (.false.)
      tfile='timeseries.tex'
   end select   

   ! Temporary Variables allocations
   allocate(x,source=self.x)
   allocate(y,source=self.y)

   ! Weighting 
   select case (lWeighting)
    case (.true.)
      allocate(sy,source=self.sy)
    case (.false.)
      allocate(sy,mold=self.sy)
      sy=one
   end select

   select case (model) 
    case ('Mean','MEAN','mean')
      allocate(N.res,source=(y-mean(y)))
    case default  
      N=BuildAndSolveNormalEquations(x=x,y=y,sy=sy,xshift=shift,model=model,timeunit=time,logfile=lfile,texfile=tfile)
   end select
   !Output 
   allocate(output)
   allocate(output.x,source=self.x)
   allocate(output.y,source=N.res)
   allocate(output.sy,source=self.sy)
   allocate(output.model,source=N.model)
   deallocate(x,y,sy)

  end function Remove  

  function UnNormalizeTheSeries(self,mean,stdev) result (output)

    implicit none 
    class(series),intent(in) :: self 
    real(kind=kpr),intent(in) :: mean
    real(kind=kpr),intent(in) :: stdev 
    class(series),allocatable :: output

    allocate(output,source=self)
    output%y=(output%y*stdev)+mean

  end function UnNormalizeTheSeries

  function NormalizeTheSeries(self) result(output)  

    implicit none 
    class(series),intent(in) :: self 
    class(series),allocatable :: output

    !Local Variables
    class(series),allocatable :: temp 
    real(kind=kpr) :: mn,st

    allocate(temp,source=self)
    mn=mean(self.y)
    st=stdev(self.y)
    temp.y=(temp.y-mn)/st 
    allocate(output,source=temp)
    deallocate(temp)

  end function NormalizeTheSeries

  real(kind=kpr) function GetTheWRMS(self) result(WRMSts) 

    implicit none 
    class(series),intent(in) :: self 

    !Local variables
    WRMSts=wrms(self%y,self%sy)

  end function GetTheWRMS

  function GetTheLSQLinearTrend(self,ReferenceEpoch) result(Solution)

    implicit none 
    class(series),intent(in) :: self 
    real(kind=kpr),intent(in),optional :: ReferenceEpoch
    real(kind=kpr),dimension(2) :: Solution

    ! Local Variables 
    real(kind=kpr) :: Offset,Trend 
     
    call GetLinearTrend(self%x,self%y,self%sy,Offset,Trend,ReferenceEpoch)
    Solution(1)=Trend
    Solution(2)=Offset 

  end function GetTheLSQLinearTrend 

  function PremoliRecursiveFiltering(self,kind) result(filt)

    implicit none 

    class(series),intent(in) :: self 
    type(series) :: filt

    integer :: i 
    character(len=*),intent(in),optional :: kind
    character(len=:),allocatable :: typeOfTimeSeries
    real(kind=kpr),allocatable :: temp(:) 

  ! Handle Optional
    if(present(kind).eqv..false.)then 
     allocate(character(len=2)::typeOfTimeSeries)
     typeOfTimeSeries='xy' 
    elseif(present(kind).eqv..true.)then
     allocate(character(len=3)::typeOfTimeSeries)
     typeOfTimeSeries='xye' 
    endif 
    call AllocateTimeSeries(ts=filt,nel=size(self%x),kind=typeOfTimeSeries)
    
    select case (typeOfTimeSeries)
     case ('xy')
      filt%x=self%x
      filt%y=self%y
     case ('xye')
      filt%x=self%x
      filt%y=self%y
      filt%y=self%sy
     case ('default')
       stop' [PremoliRecursiveFiltering] select case Error'
    end select   

    ! Allocate temp 
    allocate(temp,source=filt%y)
    do i=3,size(self%x)
       filt%y(i)=temp(i)-temp(i-1)&
                -(temp(i-1)-temp(i-2))
    enddo
    ! Deallocate tempo
    deallocate(temp)

  end function PremoliRecursiveFiltering

  function FindInBetween(self,beg,end) result(chunk)

    implicit none 

    class(series),intent(in) :: self
    real(kind=kpr),intent(in) :: beg,end
    real(kind=kpr),dimension(:),allocatable :: chunk 

    !Local Variables
    integer :: i,j,icounter

    icounter=0
    do i=1,size(self%x)
     if (self%x(i).ge.beg.and.self%x(i).lt.end) then
        icounter=icounter+1
     endif
    enddo
    !Allocate
    allocate(chunk(icounter))
    chunk=zero 
    !Fill-in
    icounter=0
    do i=1,size(self%x)
     if (self%x(i).ge.beg.and.self%x(i).lt.end) then
         icounter=icounter+1
         chunk(icounter)=self%y(i) 
     endif
     if(icounter.eq.size(chunk)) exit 
    enddo

  end function FindInBetween   

  function FilterTime(tsIn,timethreshold) result(tsOut) 

   implicit none 
   type(series),intent(in) :: tsIn
   type(series) :: tsOut
   real(kind=kpr),intent(in) :: timethreshold 

   !Local Variables 
   integer :: ife,i,j,k,n,ierr,icounter 

   !Handle Formal Errors
   ife=0
   if(allocated(tsIn%sy)) ife=1

   icounter=0
   do i=1,size(tsIn%x)
    if(tsIn%x(i).ge.timethreshold)then
     icounter=icounter+1
    endif
   enddo 
   select case (ife)
    case (0)
      call AllocateTimeSeries(ts=tsOut,nel=icounter,kind='xy')
    case (1)
      call AllocateTimeSeries(ts=tsOut,nel=icounter,kind='xye')
   end select   
   icounter=0
   do i=1,size(tsIn%x)
    if(tsIn%x(i).ge.timethreshold)then
     icounter=icounter+1
     tsOut%x(icounter)=tsIn%x(i)
     tsOut%y(icounter)=tsIn%y(i)
     if(ife.eq.1) tsOut%sy(icounter)=tsIn%sy(i)
    endif
   enddo 

  end function FilterTime

  subroutine PadWithZerosAndInterpolateTimeSeries(BiTime,tsIn,tsOut)

   implicit none 
   real(kind=kpr),intent(in) :: BiTime(:,:)
   type(series),intent(in) :: tsIn
   type(series),intent(out) :: tsOut

   !Local Variables 
   real(kind=kpr),dimension(:),allocatable :: xp,yp

   call AllocateTimeSeries(tsOut,size(BiTime,dim=1))

   call ZeroPadAndInterpolation(tsIn%x,tsIn%y,xp,yp,BiTime)

   tsOut%x=xp
   tsOut%y=yp

   deallocate(xp,yp)

  end subroutine PadWithZerosAndInterpolateTimeSeries

  subroutine PadTimeSeriesWithZeros(BiTime,tsIn,tsOut)

   implicit none 
   real(kind=kpr),intent(in) :: BiTime(:,:)
   type(series),intent(in) :: tsIn
   type(series),intent(out) :: tsOut

   !Local Variables 
   real(kind=kpr),dimension(:),allocatable :: xp,yp

   call AllocateTimeSeries(tsOut,size(BiTime,dim=1))

   call ZeroPadding(tsIn%x,tsIn%y,xp,yp,BiTime)

   tsOut%x=xp
   tsOut%y=yp

   deallocate(xp,yp)

  end subroutine PadTimeSeriesWithZeros 

  subroutine ZeroPadAndInterpolation(x,y,xp,yp,BiTime)

   implicit none 
   real(kind=kpr),intent(in) :: x(:),y(:),BiTime(:,:)
   real(kind=kpr),intent(inout),allocatable :: xp(:),yp(:)

   ! Local Variables
   integer :: i,j,k,l,m,icounter 
   integer :: ierr 
   real(kind=kpr),dimension(:),allocatable :: timeinterp
   real(kind=kpr) :: r
   type(series) :: tsIn,tsInterpolated
   
   allocate(xp(size(BiTime,dim=1)), &
            yp(size(BiTime,dim=1)),stat=ierr)
   if(ierr/=0)stop' [ZeroPadAndInterpolation] xp,yp Allocation Error '         
   xp=zero
   yp=zero

   ! Want to Interpolate?
   allocate(timeinterp(size(WhereToInterpolate(BiTime))))
   timeinterp=WhereToInterpolate(BiTime)
   call AllocateTimeSeries(tsIn,size(x))
   tsIn%x=x
   tsIn%y=y
   call InterpolateTimeSeries(tsIn,tsInterpolated,timeinterp)
   call DeallocateTimeSeries(tsIn)

   do i=1,size(BiTime,dim=1)
      if(BiTime(i,2).eq.two)then 
         xp(i)=BiTime(i,1) 
         yp(i)=zero
      endif    
      if(BiTime(i,2).eq.one.or.& 
         BiTime(i,2).eq.zero)then
         xp(i)=BiTime(i,1) 
         timeobs: do j=1,size(tsInterpolated%x)
          if(BiTime(i,1).eq.tsInterpolated%x(j))then 
             yp(i)=tsInterpolated%y(j)
             exit 
          endif    
         enddo timeobs
      endif    
   enddo
   call DeallocateTimeSeries(tsInterpolated)

  end subroutine ZeroPadAndInterpolation

  subroutine ZeroPadding(x,y,xp,yp,BiTime)

   implicit none 
   real(kind=kpr),intent(in) :: x(:),y(:),BiTime(:,:)
   real(kind=kpr),intent(inout),allocatable :: xp(:),yp(:)

   ! Local Variables
   integer :: i,j,k,l,m,icounter 
   integer :: ierr 
   real(kind=kpr) :: r
   
   allocate(xp(size(BiTime,dim=1)), &
            yp(size(BiTime,dim=1)),stat=ierr)
   if(ierr/=0)stop' [ZeroPadding] xp,yp Allocation Error '         
   xp=zero
   yp=zero

   do i=1,size(BiTime,dim=1)
      if(BiTime(i,2).eq.two)then 
         xp(i)=BiTime(i,1) 
         yp(i)=zero
      endif    
      if(BiTime(i,2).eq.one.or.& 
         BiTime(i,2).eq.zero)then
         xp(i)=BiTime(i,3) 
         timeobs: do j=1,size(x)
          if(BiTime(i,3).eq.x(j))then 
             yp(i)=y(j)
             exit 
          endif    
         enddo timeobs
      endif    
   enddo

  end subroutine ZeroPadding
  
  subroutine InterpolateTimeSeries(tsIn,tsOut,time)

   implicit none 

   real(kind=kpr),intent(in) :: time(:)
   type(series),intent(in) :: tsIn
   type(series),intent(out) :: tsOut

   !Local Variables 
   real(kind=kpr),dimension(:),allocatable :: splinecoeff, &
                                               xvar, & 
                                               polynomials, &
                                               interpvalue 
   integer :: i,j,h,n,ierr

   call AllocateTimeSeries(tsOut,size(time))
   tsOut%x=time

   ! Get the Cubic Splines coefficients
   n=size(tsIn%x)-1
   allocate(splinecoeff(0:(n-1)), &
            xvar(0:(n-1)), & 
            polynomials(0:(n-1)),stat=ierr)
   if(ierr/=0)stop' [InterpolateTimeSeries] Allocation Error '         
   xvar=tsIn%x
   splinecoeff=zero
   polynomials=tsIn%y
   call spline3_coef((n-1),xvar,polynomials,splinecoeff)

   ! Get the Interpolations
   h=0
   do while (h.lt.size(time))
      h=h+1
      tsOut%y(h)=spline3_eval((n-1),xvar,polynomials,splinecoeff,tsOut%x(h)) 
   enddo

  end subroutine InterpolateTimeSeries 

  function WhereToInterpolate(BiTime) result(time)

   implicit none 

   real(kind=kpr),dimension(:),allocatable :: time 
   real(kind=kpr),dimension(:,:),intent(in) :: BiTime

   !Local Variables 
   integer :: i,j,k,ierr,icounter
   
   icounter=0
   do i=1,size(BiTime,dim=1)
      if(BiTime(i,2).eq.one) icounter=icounter+1
   enddo
   allocate(time(icounter),stat=ierr)
   if(ierr/=0)stop' [WhereToInterpolate] time Allocation Error '
   icounter=0
   do i=1,size(BiTime,dim=1)
      if(BiTime(i,2).eq.one) then
         icounter=icounter+1
         time(icounter)=BiTime(i,1)
      endif   
   enddo

  end function WhereToInterpolate

  subroutine DetectGaps(x,timeunit,SamplingPeriod,tol,BiTime,x0)

   ! x, SamplingPeriod, tol are expressed in decy
   implicit none 
   real(kind=kpr),intent(in) :: x(:),SamplingPeriod,tol
   character(len=*),intent(in) :: timeunit 
   real(kind=kpr),dimension(:,:),intent(inout),allocatable :: BiTime 
   real(kind=kpr),intent(in),optional :: x0

   ! Local Variables 
   integer :: ifound,icounter,iselector,i,j,n,ierr 
   integer :: igap
   real(kind=kpr) :: k,threshold,lag,jd,et,mjd
   real(kind=kpr),dimension(:),allocatable :: diff,timeconv 
   logical :: lerr 

   select case (timeunit)
    case ('secj2000','Secj200','SECJ2000')
        iselector=1
        lag=SamplingPeriod*(real(60*60*24*365.25,kind=kpr))
        threshold=tol*(real(60*60*24*365.25,kind=kpr))
    case ('mjd','MJD')
        iselector=2
        lag=SamplingPeriod*365.25_kpr 
        threshold=tol*365.25_kpr
   end select     

   ! Time Tag Conversion 
   n=size(x)
   allocate(timeconv(n),stat=ierr)
   timeconv=zero 
   if(ierr/=0)stop' [DetectGaps] timeconv Allocation Error '
   do i=1,n
    select case (iselector) 
     case (1) ! to secJ2000
        call yr2jd(x(i),jd,lerr) 
        call jd_to_j2000(jd,timeconv(i))
     case (2) ! to MJD 
        call yr2mjd(x(i),timeconv(i),lerr) 
    end select
   enddo 

   ! Evenly Sampled Array
   icounter=0
   if(present(x0))then
    k=x0   
   else 
    k=x(1)
   endif
   do while (k.ge.x(1).and.k.le.x(size(x)))
    icounter=icounter+1
    k=k+SamplingPeriod
   end do  
   if(allocated(BiTime)) deallocate(BiTime)
   allocate(BiTime(icounter,3),stat=ierr)
   if(ierr/=0)stop' [DetectGaps] 2DTime Allocation Error'
   BiTime=zero
   BiTime(:,2)=two
   k=x(1)-SamplingPeriod
   do i=1,size(BiTime,dim=1)
    k=k+SamplingPeriod
    BiTime(i,1)=k
   enddo

   ! Compare BiTime and timeconv
   uneven: do i=1,size(timeconv)
     ifound=0
     even: do j=1,size(BiTime,dim=1)
      if(abs(BiTime(j,1)-x(i)).le.tol)then   
        if(abs(BiTime(j,1)-x(i)).eq.zero)then
           BiTime(j,3)=x(i) 
           BiTime(j,2)=zero ! No need to Interp
           ifound=1
        else
           BiTime(j,3)=x(i) 
           BiTime(j,2)=one ! Gentle Interpolation 
           ifound=1
        endif
      endif  
      if(ifound.eq.1) cycle uneven 
     enddo even 
   enddo uneven
   do i=1,size(BiTime,dim=1)
    write(unit=*,fmt='(3(2x,f12.4))')BiTime(i,1),BiTime(i,2),BiTime(i,3)
   enddo 

   ! Form Differences
   allocate(diff(n-1),stat=ierr)
   if(ierr/=0)stop' [DetectGaps] diff Allocation Error '
   diff=zero
   do i=1,size(diff) 
    diff(i)=timeconv(i+1)-timeconv(i)
   end do

   ! Detect Gaps 
   igap=0
   do i=1,size(diff)
    if(diff(i)-lag.gt.threshold) igap=igap+1
   enddo  

  end subroutine DetectGaps  

  subroutine MakeTemporalMask(mask,iBeg,iStop,ts)

   implicit none 
   integer,intent(in) :: iBeg,iStop
   type(series),intent(in) :: ts
   logical,dimension(:),allocatable,intent(inout) :: mask

   ! Local Variables 
   integer :: i,j,ierr
   integer,dimension(:),allocatable :: BegVector,EndVector

   ! Allocate mask
   if(allocated(mask)) deallocate(mask)
   allocate(mask(size(ts%x)),stat=ierr)
   if(ierr/=0)stop" mask Alloc Error "
   ! Allocate Time Tag arrays
   allocate(BegVector(size(ts%x)),EndVector(size(ts%x)),stat=ierr)
   if(ierr/=0)stop" BegVector Alloc Error "
   ! Initialize Arrays
   mask=.false.
   BegVector=iBeg
   EndVector=iStop
   ! Define Mask 
   mask=ts%x.ge.iBeg.and.ts%x.le.iStop
   ! Count 
   write(*,'(a,1x,i10,/)')' Number of Elements Selected = ',count(mask)
   ! Destroy Arrays
   deallocate(BegVector,EndVector,stat=ierr)
   if(ierr/=0)stop" BegVector Dealloc Error "

  end subroutine MakeTemporalMask 

  subroutine ApplyTemporalMask(mask,ts)

   type(series),intent(inout) :: ts
   logical,dimension(:),intent(in) :: mask

   ! Local Variables
   type(series) :: tempts
   integer :: i,j,ierr

   call AllocateTimeSeries(tempts,count(mask))
   do i=1,count(mask)
    if(mask(i).eqv..true.)then
       tempts%x(i)=ts%x(i) 
       tempts%y(i)=ts%y(i)
    endif 
   enddo
   call DeAllocateTimeSeries(ts)
   call AllocateTimeSeries(ts,count(mask))
   ts%x=tempts%x
   ts%y=tempts%y
   call DeAllocateTimeSeries(tempts)

  endsubroutine ApplyTemporalMask

  function Allocate(nel,kind) result(output)

  implicit none 
  integer,intent(in) :: nel
  character(len=*),intent(in),optional :: kind
  type(series) :: output 

  integer :: ierr 
  character(len=:),allocatable :: typeOfTimeSeries

  !Handle Optional
   if(present(kind).eqv..false.)then 
     allocate(character(len=2)::typeOfTimeSeries)
     typeOfTimeSeries='xy' 
   elseif(present(kind).eqv..true.)then
     allocate(character(len=3)::typeOfTimeSeries)
     typeOfTimeSeries='xye' 
   endif 

  call  AllocateTimeSeries(ts=output,nel=nel,kind=typeOfTimeSeries)

  end function Allocate

  subroutine AllocateTimeSeries(ts,nel,kind)

  implicit none 
  integer,intent(in) :: nel
  type(series),intent(out) :: ts
  character(len=*),intent(in),optional :: kind

  !Local Variables
  integer :: ierr 
  character(len=:),allocatable :: typeOfTimeSeries

  !Handle Optional Values
  if(present(kind).eqv..false.)then 
    allocate(character(len=2)::typeOfTimeSeries)
    typeOfTimeSeries='xy' 
  elseif(present(kind).eqv..true.)then
    allocate(character(len=3)::typeOfTimeSeries)
    typeOfTimeSeries='xye' 
  endif 

  if(allocated(ts%x))then 
     call DeAllocateTimeSeries(ts)
  endif   

  select case (typeOfTimeSeries)
   case ('xy')
    allocate(ts%x(nel),ts%y(nel),stat=ierr)
    if(ierr/=0)stop' [AllocateTimeSeries] Allocation Error '
   case ('xye')  
    allocate(ts%x(nel),ts%y(nel),ts%sy(nel),stat=ierr)
    if(ierr/=0)stop' [AllocateTimeSeries] Allocation Error '
    ts%sy=zero
  end select

  !Initialize
  ts%x=zero
  ts%y=zero 

  end subroutine AllocateTimeSeries 

  subroutine DeallocateTimeSeries(t)

  implicit none 
  type(series),intent(inout) :: t

  integer :: ierr 

  deallocate(t%x,t%y,t%sy,stat=ierr)
  if(ierr/=0)stop' [Allocatet] Allocation Error '

  end subroutine DeallocateTimeSeries 

  subroutine readts(self,fname,kind,changeSign)
   
   implicit none 

   class(series),intent(inout) :: self 
   character(len=*),intent(in) :: fname 
   character(len=*),intent(in),optional :: kind,changeSign

   integer :: i,ifn,ierr 
   integer :: nel 

  character(len=:),allocatable :: sChangeSign,typeOfTimeSeries

  logical :: lDebug

  lDebug=.false.

  !Handle Optional Values
  select case (present(changeSign))
   case (.true.)
     sChangeSign=changeSign
   case (.false.)
     sChangeSign='000'
  end select     

   if(present(kind).eqv..false.)then 
     allocate(character(len=2)::typeOfTimeSeries)
     typeOfTimeSeries='xy' 
   elseif(present(kind).eqv..true.)then
     allocate(character(len=3)::typeOfTimeSeries)
     typeOfTimeSeries='xye' 
   endif 

   call GetUnit(ifn)
   open(unit=ifn,file=fname,status='unknown',action='read',iostat=ierr)
   if(ierr/=0)stop' [readts] I/O Error '
   
   ! Allocate time series
   call GetRows(fname,nel)
   call AllocateTimeSeries(self,nel,typeOfTimeSeries)

   ! Read time series
   select case (typeOfTimeSeries)
   case ('xy') 
     do i=1,nel
      read(unit=ifn,fmt=*)self%x(i),self%y(i)
     enddo 
   case ('xye') 
     do i=1,nel
      read(unit=ifn,fmt=*)self%x(i),self%y(i),self%sy(i)
     enddo 
   end select    
   
   ! Flip Sign 
   select case (sChangeSign)
    case ('000')  
    case ('010')
      if(lDebug) write(*,*) '[readts] ... sign flipped on self.y ... '
      self%y=-self%y
    case default 
      stop'[readts select case sChangeSign option undefined]'
   end select

  end subroutine readts

  subroutine ReadTimeSeries(fname,ts,kind)
   
   implicit none 

   type(series),intent(inout) :: ts 
   character(len=*),intent(in) :: fname 
   character(len=*),intent(in),optional :: kind

   integer :: i,ifn,ierr 
   integer :: nel 

  character(len=:),allocatable :: typeOfTimeSeries

  !Handle Optional Values
   if(present(kind).eqv..false.)then 
     allocate(character(len=2)::typeOfTimeSeries)
     typeOfTimeSeries='xy' 
   elseif(present(kind).eqv..true.)then
     allocate(character(len=3)::typeOfTimeSeries)
     typeOfTimeSeries='xye' 
   endif 

   call GetUnit(ifn)
   open(unit=ifn,file=fname,status='unknown',action='read',iostat=ierr)
   if(ierr/=0)stop' [ReadTimeSeries] I/O Error '
   
   ! Allocate time series
   call GetRows(fname,nel)
   call AllocateTimeSeries(ts,nel,typeOfTimeSeries)

   ! Read time series
   select case (typeOfTimeSeries)
   case ('xy') 
     do i=1,nel
      read(unit=ifn,fmt=*)ts%x(i),ts%y(i)
     enddo 
   case ('xye') 
     do i=1,nel
      read(unit=ifn,fmt=*)ts%x(i),ts%y(i),ts%sy(i)
     enddo 
   end select    

  end subroutine ReadTimeSeries

  subroutine WriteTimeSeries(self,fname,timearray,from,to)

   implicit none 

   class(series),intent(in) :: self
   character(len=*),intent(in) :: fname 
   real(kind=kpr),intent(in),optional :: timearray(:,:)
   real(kind=kpr),intent(in),optional :: from,to

   integer :: i,ifn,ierr 
   real(kind=kpr) :: fromtime,totime

   !Handle Optional Argumenself 
   select case (present(from))
    case (.true.)
      fromtime=from
    case (.false.)
      fromtime=self%x(1)
   end select    
   select case (present(to))
    case (.true.)
      totime=to
    case (.false.)
      totime=self%x(size(self%x))
   end select    

   call GetUnit(ifn)
   open(unit=ifn,file=fname,status='unknown',action='write',iostat=ierr)
   call WriteTimeSeriesHeader(iUnit=ifn,unit='mm',nElements=size(self%x))
   do i=1,size(self%x)
    if (self%x(i).ge.fromtime.and.&
        self%x(i).le.totime) then  
         if(allocated(self%sy))then 
            select case (present(timearray))
             case (.false.)
               write(unit=ifn,fmt=10)self%x(i),self%y(i),self%sy(i)
             case (.true.)
               if (timearray(i,2).eq.two) then
                  write(unit=ifn,fmt=12)self%x(i),self%y(i),self%sy(i),'R'
               else   
                  write(unit=ifn,fmt=12)self%x(i),self%y(i),self%sy(i),'O'
               endif
            end select
         else
            write(unit=ifn,fmt=11)self%x(i),self%y(i)
         endif   
    endif 
   enddo
   close(unit=ifn)

   10 format(f20.12,2(2x,f28.16))
   11 format(f20.12,1(2x,f28.16))
   12 format(f20.12,2(2x,f28.16),8x,a1)

  end subroutine WriteTimeSeries 

  subroutine WriteModelandResidualsTimeSeries(ts,model,residuals,fname)

   implicit none 
   real(kind=kpr),dimension(:),intent(in) :: model,residuals
   type(series),intent(in) :: ts
   character(len=*),intent(in) :: fname 

   integer :: i,ifn,ierr 

   call GetUnit(ifn)
   open(unit=ifn,file=fname,status='unknown',action='write',iostat=ierr)
   do i=1,size(ts%x)
    if(allocated(ts%sy))then 
       write(unit=ifn,fmt=201)ts%x(i),ts%y(i),residuals(i),model(i),ts%sy(i)
    else
       write(unit=ifn,fmt=202)ts%x(i),ts%y(i),residuals(i),model(i)
    endif   
   enddo
   close(unit=ifn)

   201 format(f26.19,2x,e26.19,2x,e26.19,2x,e26.19,2x,e26.19)
   202 format(f26.19,2x,e26.19,2x,e26.19,2x,e26.19)

  end subroutine WriteModelandResidualsTimeSeries

  subroutine WriteResidualTimeSeries(ts,residuals,fname)

   implicit none 
   real(kind=kpr),dimension(:),intent(in) :: residuals
   type(series),intent(in) :: ts
   character(len=*),intent(in) :: fname 

   integer :: i,ifn,ierr 

   call GetUnit(ifn)
   open(unit=ifn,file=fname,status='unknown',action='write',iostat=ierr)
   do i=1,size(ts%x)
    if(allocated(ts%sy))then 
       write(unit=ifn,fmt=201)ts%x(i),ts%y(i),residuals(i),ts%sy(i)
    else
       write(unit=ifn,fmt=202)ts%x(i),ts%y(i),residuals(i)
    endif   
   enddo
   close(unit=ifn)

   201 format(f26.19,2x,e26.19,2x,e26.19,2x,e26.19)
   202 format(f26.19,2x,e26.19,2x,e26.19)

  end subroutine WriteResidualTimeSeries

  subroutine SelectTimeSeries(ts,beg,stp)

   implicit none
   real(kind=kpr),intent(in) :: beg,stp 
   type(series),intent(inout) :: ts

   !Local Variables
   type(series) :: temp
   integer :: i,counter

   counter=0
   do i=1,size(ts%x) 
     if(ts%x(i).ge.beg.and.ts%x(i).le.stp)then 
       counter=counter+1
     endif    
   enddo

   call AllocateTimeSeries(temp,counter)

   ! Clip the time series 
   counter=0
   do i=1,size(ts%x) 
     if(ts%x(i).ge.beg.and.ts%x(i).le.stp)then 
       counter=counter+1
       temp%x(counter)=ts%x(i)
       temp%y(counter)=ts%y(i)
       temp%sy(counter)=ts%sy(i)
     endif    
   enddo

   call DeAllocateTimeSeries(ts)
   call AllocateTimeSeries(ts,counter)
   ts%x=temp%x
   ts%y=temp%y
   ts%sy=temp%sy 
   call DeAllocateTimeSeries(temp)

  end subroutine SelectTimeSeries

  function DeMean(tsIn) result(tsOut)
   
   implicit none

   type(series),intent(in) :: tsIn 
   type(series) :: tsOut 

   !Local Variables 
   real(kind=kpr) :: m 
   
   call AllocateTimeSeries(tsOut,size(tsIn%x)) 
   m=mean(tsIn%y)
   tsOut%x=tsIn%x
   tsOut%y=tsIn%y-m

  end function DeMean

  function GenTCH(ts1,ts2,ts3) result(tch) 

   implicit none 

   type(series),intent(in) :: ts1,ts2,ts3
   real(kind=kpr),dimension(6) :: tch

   !Local Variables
   type(series) :: d12,d13,d23
   integer :: i
   real(kind=kpr) :: fmin,a,b,s11,s22,s12,root 
   real(kind=kpr) :: r33 
   real(kind=kpr) :: lag,x,R(5),F(2),SmallestPositiveRoot
   real(kind=kpr),dimension(0:6) :: c
   real(kind=kpr),dimension(0:2) :: bcoeff
   real(kind=kpr),dimension(:,:),allocatable :: M
   real(kind=kpr),dimension(:),allocatable :: RealRoots,d
   complex(kind=kpr),dimension(:),allocatable :: solution

   !ts3 is the reference ts

   !Differentiation
   d13=Difference(ts1,ts3)
   d23=Difference(ts2,ts3) 

   !Sample Var/Covariances 
   s11=stdev(d13%y)**2
   s22=stdev(d23%y)**2
   s12=cova(d13%y,d23%y)
   
   c=EvaluateCoefficients(s11,s22,s12)
   print*,''
   write(*,*)' + Coefficients of a Sixth-Order Polynomial +'
        do i=0,size(c)-1
           write(*,'(3x,a,1x,i4,3x,g24.12)') ' Order ',i,c(i)
        end do
   write(*,*)' - Coefficients of a Sixth-Order Polynomial -'

   allocate(d(0:size(c)-2),solution(size(c)-2))
   solution=zero
   d=zero
   d=c(1:6)
   d=d/d(5)
   allocate(M(5,5))
   M=zero
   M=FormCompanionMatrix(d(0:4))
   solution=Roots(M)
   call AnalyzeRoots(solution,RealRoots,SmallestPositiveRoot)
   print*,''
   write(*,*)' + Non-Trivial Real Roots +'
     do i=1,size(RealRoots)
         write(*,'(5x,i4,1x,g24.12)') i,RealRoots(i)
     end do
   write(*,*)' - Non-Trivial Real Roots -'
   print*,''
   if(isNan(SmallestPositiveRoot).eqv..true.) then 
     fmin=zero 
   else
     fmin=SmallestPositiveRoot
   endif  
   tch=zero
   print*,''
   write(*,*)' + Smallest Positive Root +'
           write(*,'(3x,a,1x,g24.12)') ' fmin = ',fmin
   write(*,*)' - Smallest Positive Root -'
   !Determine b coefficients 
   bcoeff=EvaluateBCoefficients(s11,s22,s12,fmin)
   r33=-bcoeff(1)/bcoeff(2)
   R=EvaluateRCoefficients(s11,s22,s12,r33,fmin)
   tch(1)=R(1) 
   tch(2)=R(3)
   tch(3)=r33
   tch(4)=R(2) ! Cova(1,2)
   tch(5)=R(4) ! Cova(1,3)
   tch(6)=R(5) ! Cova(2,3)
  end function GenTCH 

  function EvaluateRCoefficients(s11,s22,s12,r33,fmin) result(R)

   implicit none 
   real(kind=kpr),intent(in) :: s11,s22,s12,r33
   real(kind=kpr),intent(in) :: fmin
   real(kind=kpr),dimension(5):: R 
   real(kind=kpr) :: den,detS,r11,r12,r22,r13,r23
   real(kind=kpr) :: a20,a02,a11,a10,a00 

   detS=abs(s11*s22-s12**2) 
   
   a20=two*sqrt(detS)+fmin*s22
   a02=two*sqrt(detS)+fmin*s11
   a11=sqrt(detS)-fmin*s12 
   a10=sqrt(detS)*(two*r33+s12)
   a00=sqrt(detS)*(two*r33**2+(s12+r33)**2)-fmin*r33*detS 

   den=a20*a02-a11**2
   r13=r33-(a10*(a02-a11))/den 
   r23=r33-(a10*(a20-a11))/den 

   r11=s11-r33+two*r13 
   r12=s12-r33+r13+r23   
   r22=s22-r33+two*r23  

   R(1)=r11
   R(2)=r12
   R(3)=r22
   R(4)=r13
   R(5)=r23

  end function EvaluateRCoefficients

  function EvaluateBCoefficients(s11,s22,s12,fmin) result(B)

   implicit none 
   real(kind=kpr),intent(in) :: s11,s22,s12
   real(kind=kpr),intent(in) :: fmin
   real(kind=kpr),dimension(0:2) :: b
   real(kind=kpr) :: detS

   detS=abs(s11*s22-s12**2) 
   b(2)=sqrt(detS)+two*(s11+s22-s12)*fmin+three*sqrt(detS)*fmin**2  
   b(1)=-sqrt(detS)*s12-(2*s12**2+three*detS/two)*fmin+ &
        -sqrt(detS)*(s11+s22)*fmin**2-(detS*fmin**3)/two
   b(0)=sqrt(detS)*s12**2+(s12**2*(s11+s22))*fmin+sqrt(detS)*(s12*fmin)**2     

  end function EvaluateBCoefficients

  function tch(ts1,ts2,ts3) 

   implicit none 

   type(series),intent(in) :: ts1,ts2,ts3
   real(kind=kpr) :: tch(3)
   type(series) :: d12,d13,d23

   real(kind=kpr) :: v12,v13,v23 
   real(kind=kpr),dimension(3) :: L(3)  
   real(kind=kpr),dimension(3,3) :: A(3,3)  

   !Differentiation
   d12=Difference(ts1,ts2)
   d13=Difference(ts1,ts3)
   d23=Difference(ts2,ts3) 
   call WriteTimeSeries(d12,'d12.dat')
   call WriteTimeSeries(d13,'d13.dat')
   call WriteTimeSeries(d23,'d23.dat')

   !Sample Variances 
   v12=stdev(d12%y)**2
   v13=stdev(d13%y)**2
   v23=stdev(d23%y)**2
  
   !Jacobi Matrix
   A=one
   A(1,3)=zero
   A(2,2)=zero 
   A(3,1)=zero  

   !Right-Hand Term
   L(1)=v12
   L(2)=v13
   L(3)=v23 

   ! Solve
   call inverse(A,3)
   tch=matmul(A,L)

  end function tch

  subroutine SelectCommonObs(ts1,ts2,comm1,comm2)

   implicit none 

   type(series),intent(in) :: ts1,ts2
   type(series),intent(out) :: comm1,comm2

   integer :: nel 
   integer :: i,j,ifound
   integer :: i1,i2,ife 

   nel=0
   outer1: do i=1,size(ts1%x)
    inner1: do j=1,size(ts2%x)
     if(ts1%x(i).eq.ts2%x(j))then 
      nel=nel+1
      cycle outer1  
     endif     
     enddo inner1 
   enddo outer1 

   !Handle Formal Errors
   i1=0
   i2=0
   if(allocated(ts1%sy)) i1=1
   if(allocated(ts2%sy)) i2=1
   if(i1.eq.1.and.i2.eq.1)then
    ife=1
   else
    ife=0   
   endif    
   if(ife.eq.0) allocate(comm1%x(nel),comm1%y(nel),& 
                         comm2%x(nel),comm2%y(nel))
   if(ife.eq.1) allocate(comm1%x(nel),comm1%y(nel),comm1%sy(nel),& 
                         comm2%x(nel),comm2%y(nel),comm2%sy(nel))
   nel=0
   outer2: do i=1,size(ts1%x)
    inner2: do j=1,size(ts2%x)
     if(ts1%x(i).eq.ts2%x(j))then 
      nel=nel+1
      comm1%x(nel)=ts1%x(i)
      comm1%y(nel)=ts1%y(i)
      comm2%x(nel)=ts2%x(j)
      comm2%y(nel)=ts2%y(j)
      if(ife.eq.1)then
        comm1%sy(nel)=ts1%sy(i)
        comm2%sy(nel)=ts2%sy(j)
      endif    
      cycle outer2  
     endif     
     enddo inner2 
   enddo outer2 

  end subroutine SelectCommonObs  

  pure function Difference(ts1,ts2) result(diff)

   implicit none 

   type(series),intent(in) :: ts1,ts2
   type(series) :: diff

   integer :: nel 
   integer :: i,j,ifound
   integer :: i1,i2,ife 

   nel=0
   outer1: do i=1,size(ts1%x)
    inner1: do j=1,size(ts2%x)
     if(ts1%x(i).eq.ts2%x(j))then 
      nel=nel+1
      cycle outer1  
     endif     
     enddo inner1 
   enddo outer1 

   !Handle Formal Errors
   i1=0
   i2=0
   if(allocated(ts1%sy)) i1=1
   if(allocated(ts2%sy)) i2=1
   if(i1.eq.1.and.i2.eq.1)then
    ife=1
   else
    ife=0   
   endif    
   if(ife.eq.0) allocate(diff%x(nel),diff%y(nel))
   if(ife.eq.1) allocate(diff%x(nel),diff%y(nel),diff%sy(nel))
   nel=0
   outer2: do i=1,size(ts1%x)
    inner2: do j=1,size(ts2%x)
     if(ts1%x(i).eq.ts2%x(j))then 
      nel=nel+1
      diff%x(nel)=ts1%x(i)
      diff%y(nel)=ts1%y(i)-ts2%y(j)
      if(ife.eq.1)then
       diff%sy(nel)=sqrt(ts1%sy(i)**2+ts2%sy(j)**2)
      endif    
      cycle outer2  
     endif     
     enddo inner2 
  enddo outer2 

  end function Difference  

  subroutine SortVector(a,indx)

  ! SortVector does an indexed heap ascending sort of the
  ! vector stored in a.
  ! The sorting is not actually carried out.  Rather an index array is
  ! created which defines the sorting.  This array may be used to sort
  ! or index the array, or to sort or index related arrays keyed on the
  ! original array.
  ! Once the index array is computed, the sorting can be carried out
  ! implicitly as follows: a(indx(1:n)).

  implicit none
  real ( kind = kpr ) , intent ( in ) , dimension ( : ) :: a
  integer , dimension ( : ) , intent ( inout ) , allocatable ::  indx

  ! Local Variables 
  integer ::   n
  real ( kind = kpr) :: aval
  integer ::  i
  integer ::  indxt
  integer ::  ir
  integer ::  j
  integer ::  l

  if ( allocated ( indx ) .eqv. .true. ) deallocate ( indx ) 

  n = size ( a ) 

  allocate ( indx ( n ) ) 

  if ( allocated ( indx ) .neqv. .true. ) stop ' [SortVector] indx Allocation Error'

  if ( n < 1 ) then
    return
  end if

  do i = 1, n
    indx(i) = i
  end do

  if ( n == 1 ) then
    return
  end if

  l = n / 2 + 1
  ir = n

   lout: do
     if ( 1 < l ) then
       l = l - 1
       indxt = indx(l)
       aval = a(indxt)
     else
       indxt = indx(ir)
       aval = a(indxt)
       indx(ir) = indx(1)
       ir = ir - 1
       if ( ir == 1 ) then
         indx(1) = indxt
         exit
       end if
     end if
     i = l
     j = l + l
     lwhile: do while ( j <= ir )
       if ( j < ir ) then
         if ( a(indx(j)) < a(indx(j+1)) ) then
           j = j + 1
         end if
       end if
       if ( aval < a(indx(j)) ) then
         indx(i) = indx(j)
         i = j
         j = j + j
       else
         j = ir + 1
       end if
     end do lwhile 
     indx(i) = indxt
   end do lout

  end subroutine SortVector

  function SingularSpectrumReconstruction(self,m,nComponents,logfile) result(output)

    implicit none 

    class(series),intent(in) :: self
    integer,intent(in) :: m,nComponents
    character(len=*),intent(in),optional :: logfile
    class(series),allocatable :: output 
    
    !Local Variables 
    type(series) :: temp 
    
    !Handle Output
    call AllocateTimeSeries(ts=temp,nel=size(self.x),kind='xy')
    temp.x=self.x

    !Apply SSA to self without Normalizing the time series
    call SSA(x=self.x,y=self.y,m=m,ncmp=nComponents,yrec=temp.y,normalization=.false.,logfile=logfile)

    !Allocate the output
    allocate(output,source=temp)

    !Deallocate temp
    call DeAllocateTimeSeries(t=temp)

  end function SingularSpectrumReconstruction

  subroutine GapDetector(self,timetag,samplingperiod,igap,thisIsWhereTheGapsAre)

   implicit none 
   
   class(series),intent(in) :: self 
   character(len=*),intent(in) :: timetag
   real(kind=kpr),intent(in) :: samplingperiod
   integer, intent(out) :: igap
   real(kind=kpr),allocatable,intent(out),optional :: thisIsWhereTheGapsAre(:,:)

   !Local Variables
   integer :: ierr,ifound,icounter 
   integer :: i,j
   real(kind=kpr) :: k,tol,threshold
   real(kind=kpr),allocatable :: diff(:)
   real(kind=kpr),allocatable :: BiTime(:,:) 

   tol=1.0e-05_kpr

   !Construct an Evenly Sampled Array
   icounter=0
   k=self%x(1)
   do while (k.ge.self%x(1).and.k.le.self%x(size(self%x)))
    icounter=icounter+1
    k=k+SamplingPeriod
   end do  
   allocate(BiTime(icounter,3),stat=ierr)
   if(ierr/=0)stop' [GapDetector] 2DTime Allocation Error'
   BiTime=zero
   BiTime(:,2)=two
   k=self%x(1)-SamplingPeriod
   do i=1,size(BiTime,dim=1)
    k=k+SamplingPeriod
    BiTime(i,1)=k
   enddo
    
   !Compare BiTime and self%x 
   uneven: do i=1,size(self%x)
     ifound=0
     even: do j=1,size(BiTime,dim=1)
      if(abs(BiTime(j,1)-self%x(i)).le.tol)then   
        if(abs(BiTime(j,1)-self%x(i)).eq.zero)then
           BiTime(j,3)=self%x(i) 
           BiTime(j,2)=zero ! No need to Interp
           ifound=1
        else
           BiTime(j,3)=self%x(i) 
           BiTime(j,2)=one ! Gentle Interpolation 
           ifound=1
        endif
      endif  
      if(ifound.eq.1) cycle uneven 
     enddo even 
   enddo uneven

   !Count the Gaps 
   igap=0
   do i=1,size(BiTime,dim=1)
      if(BiTime(i,2).eq.one.or.BiTime(i,2).eq.two)then 
        igap=igap+1
      endif  
   enddo 

   !Handle Optional Arguments
   select case (present(thisIsWhereTheGapsAre))
    case (.true.) 
      allocate(thisIsWhereTheGapsAre,source=BiTime)
   endselect
   deallocate(BiTime)

   !! Form Differences
   !allocate(diff(size(self%x)-1),stat=ierr)
   !if(ierr/=0)stop' [DetectGaps] diff Allocation Error '
   !diff=zero
   !do i=1,size(diff) 
    !diff(i)=self%x(i+1)-self%x(i)
   !end do

   !! Detect Gaps 
   !threshold=0.1_kpr
   !igap=0
   !do i=1,size(diff)
    !if(abs(diff(i)-samplingperiod).gt.threshold)then 
      !igap=igap+1
      !print*,diff(i)
    !endif  
   !enddo  

  end subroutine GapDetector 

  subroutine ConvertTimeTag(self,from,to)

   implicit none 

   class(series),intent(inout) :: self
   character(len=*),intent(in) :: from,to 

   ! Local Variables
   integer ::  i,n,ial
   real(kind=kpr) :: timetag
   character(len=:),allocatable :: StringSelector
   logical :: lError 

  StringSelector=from//to
  n=size(self%x)
  select case (StringSelector)
   case ('MJDYR','mjdyr')
        do i=1,n
         call mjd2yr(self%x(i),timetag,lError) 
         self%x(i)=timetag
        enddo 
   case ('YRMJD','yrmjd')
        do i=1,n
         call yr2mjd(self%x(i),timetag,lError) 
         self%x(i)=timetag
        enddo 
  end select      

  end subroutine ConvertTimeTag
 
  function FillTheGapsWithTheReconstructedValues(self,Reconstructed,thisIsWhereTheGapsAre) result(output)

   implicit none 
   class(series),intent(in) :: self
   class(series),intent(in) :: Reconstructed 
   real(kind=kpr),intent(in),dimension(:,:) :: thisIsWhereTheGapsAre
   type(series) :: output 

   !Local Variables
   integer ::  i,j,n

   if(.not. allocated(self.sy))stop' [FillTheGapsWithTheReconstructedValues] self.sy not allocated'
   !Allocate output
   n=size(thisIsWhereTheGapsAre,dim=1)
   call AllocateTimeSeries(ts=output,nel=n,kind='xye')
   !Fill it in
   output%x=thisIsWhereTheGapsAre(:,1)
   output%y=zero
   output%sy=zero
   do i=1,n
      if (thisIsWhereTheGapsAre(i,2).eq.two)then  
        do j=1,size(Reconstructed%x)
           if(Reconstructed%x(j).eq.thisIsWhereTheGapsAre(i,1))then 
              output%y(i)=Reconstructed%y(j)
              output%sy(i)=100.0_kpr
           endif   
        enddo   
      else
        do j=1,size(self%x)
           if(self%x(j).eq.thisIsWhereTheGapsAre(i,1))then 
              output%y(i)=self%y(j)
              output%sy(i)=self%sy(j)
           endif   
        enddo   
      endif
   enddo 

  end function FillTheGapsWithTheReconstructedValues

  function FillTheGapsWithZeros(self,thisIsWhereTheGapsAre) result(output)

   implicit none 

   class(series),intent(in) :: self
   real(kind=kpr),intent(in),dimension(:,:) :: thisIsWhereTheGapsAre
   type(series) :: output 

   !Local Variables
   integer ::  i,j,n

   !Allocate output
   n=size(thisIsWhereTheGapsAre,dim=1)
   call AllocateTimeSeries(ts=output,nel=n,kind='xye')
   !Fill it in
   output%x=thisIsWhereTheGapsAre(:,1)
   output%y=zero
   output%sy=zero
   do i=1,n
      if (thisIsWhereTheGapsAre(i,2).eq.zero.or.& 
          thisIsWhereTheGapsAre(i,2).eq.one)then 
        do j=1,size(self%x)
           if(self%x(j).eq.thisIsWhereTheGapsAre(i,1))then 
              output%y(i)=self%y(j)
              output%sy(i)=self%sy(j)
           endif   
        enddo 
      endif
   enddo 

  end function FillTheGapsWithZeros 

  function Sort(self,key) result(output) 

   implicit none

   class(series),intent(inout) :: self
   character(len=*),intent(in) :: key
   class(series),allocatable :: output 


   ! Local Variables
   integer ::  i,n,ial
   integer, dimension(:), allocatable :: indTime

  ! Sort in an increasing the first column
  ! observations and we should apply the
  ! sorting to the other elements on the
  ! remaining columns.

  n=size(self%x)
  allocate(indTime(n),stat=ial)
  if(ial.ne.0)stop' [sortThysType] indTime allocation Error '
  indTime=0
  select case (key)
   case ('x','1')
    call SortVector(self%x,indTime)
   case ('y','2')
    call SortVector(self%y,indTime)
   case ('s','3')
    call SortVector(self%sy,indTime)
  end select  
  ! Substitution
  self%x=self%x(indTime)
  self%y=self%y(indTime)
  self%sy=self%sy(indTime)

  allocate(output,source=self)

 end function Sort 
 
 subroutine writeSREFPositionTimeSeriesHeader(self,iUnit)

   implicit none 
   class(positionTimeSeries) :: self 
   integer,intent(in) :: iUnit

   ! Local Variables
   integer :: dat(8) 
   integer :: nElements
   integer(kind=4) :: rslt
   character(len=:),allocatable :: cuteString
   character(len=:),allocatable :: friendly,iso_fmt
   character(len=max_hostnam_length+1) :: HostName

   friendly='%W, %L %d, %Y %H:%m:%s %N'
   iso_fmt='%Y-%M-%DT%h:%m:%s.%x%z'
   cuteString='~oO◜.◝oO~~oO◜.◝Oo~~oO◜.◝Oo~~oO◜.◝Oo~~oO◜.◝Oo~~oO◜.◝Oo~'

   nElements=size(self.e.x)
   call date_and_time(values=dat)
   write(unit=iUnit,fmt='("#",1x,a)') cuteString
   write(unit=iUnit,fmt='("#",1x,a)') ' File created on  ' 
   write(unit=iUnit,fmt='("#",2x,a)')fmtdate(dat,iso_fmt)
   write(unit=iUnit,fmt='("#",2x,a)')fmtdate(dat,friendly)
   rslt=hostnam(HostName)
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' The HOST is ', trim(HostName) 
   write(unit=iUnit,fmt='("#",1x,a,1x,f9.4)')   ' Site Latitude             = ' , self.lat 
   write(unit=iUnit,fmt='("#",1x,a,1x,f9.4)')   ' Site Longitude            = ' , self.lon 
   write(unit=iUnit,fmt='("#",1x,a,2x,a)')   ' Technique                 = ' , self.technique 
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' First Column              = ' , ' Time [decy] ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Second Column             = ' , ' East ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Third Column              = ' , ' North '  
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Fourth Column             = ' , ' Radial ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,i10)') ' Number of Observations    = ' , nElements
   write(unit=iUnit,fmt='("#",1x,a)') cuteString

 end subroutine writeSREFPositionTimeSeriesHeader


 subroutine writePositionTimeSeriesHeader(self,iUnit)

   implicit none 
   class(positionTimeSeries) :: self 
   integer,intent(in) :: iUnit

   ! Local Variables
   integer :: dat(8) 
   integer :: nElements
   integer(kind=4) :: rslt
   character(len=:),allocatable :: cuteString
   character(len=:),allocatable :: friendly,iso_fmt
   character(len=max_hostnam_length+1) :: HostName

   friendly='%W, %L %d, %Y %H:%m:%s %N'
   iso_fmt='%Y-%M-%DT%h:%m:%s.%x%z'
   cuteString='~oO◜.◝oO~~oO◜.◝Oo~~oO◜.◝Oo~~oO◜.◝Oo~~oO◜.◝Oo~~oO◜.◝Oo~'

   nElements=size(self.e.x)
   call date_and_time(values=dat)
   write(unit=iUnit,fmt='("#",1x,a)') cuteString
   write(unit=iUnit,fmt='("#",1x,a)') ' File created on  ' 
   write(unit=iUnit,fmt='("#",2x,a)')fmtdate(dat,iso_fmt)
   write(unit=iUnit,fmt='("#",2x,a)')fmtdate(dat,friendly)
   rslt=hostnam(HostName)
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' The HOST is ', trim(HostName) 
   write(unit=iUnit,fmt='("#",1x,a,1x,f9.4)')   ' Site Latitude             = ' , self.lat 
   write(unit=iUnit,fmt='("#",1x,a,1x,f9.4)')   ' Site Longitude            = ' , self.lon 
   write(unit=iUnit,fmt='("#",1x,a,2x,a)')   ' Technique                 = ' , self.technique 
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' First Column              = ' , ' Time [decy] ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Second Column             = ' , ' East ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Third Column              = ' , ' North '  
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Fourth Column             = ' , ' Radial ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Fifth Column              = ' , ' sigma East ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Sixth Column              = ' , ' sigma North '  
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Seventh Column            = ' , ' sigma Radial ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Eighth Column             = ' , ' c(E,N) ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Ninth Column              = ' , ' c(E,U) '  
   write(unit=iUnit,fmt='("#",1x,a,1x,a)')   ' Tenth Column              = ' , ' c(U,N) ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,i10)') ' Number of Observations    = ' , nElements
   write(unit=iUnit,fmt='("#",1x,a)') cuteString

 end subroutine writePositionTimeSeriesHeader

 subroutine WriteTimeSeriesHeader(iUnit,unit,nElements)

   implicit none 
   integer,intent(in) :: iUnit
   character(len=*),intent(in) :: unit 
   integer,intent(in) :: nElements

   ! Local Variables
   integer :: dat(8) 
   integer(kind=4) :: rslt
   character(len=:),allocatable :: cuteString
   character(len=:),allocatable :: friendly,iso_fmt
   character(len=max_hostnam_length+1) :: HostName

   friendly='%W, %L %d, %Y %H:%m:%s %N'
   iso_fmt='%Y-%M-%DT%h:%m:%s.%x%z'
   cuteString='~oO◜.◝oO~~oO◜.◝Oo~~oO◜.◝Oo~~oO◜.◝Oo~~oO◜.◝Oo~~oO◜.◝Oo~'

   call date_and_time(values=dat)
   write(unit=iUnit,fmt='("#",1x,a,a)') waves, waves 
   write(unit=iUnit,fmt='("#",1x,a)') ' File created on  ' 
   write(unit=iUnit,fmt='("#",2x,a)')fmtdate(dat,iso_fmt)
   write(unit=iUnit,fmt='("#",2x,a)')fmtdate(dat,friendly)
   rslt=hostnam(HostName)
   write(unit=iUnit,fmt='("#",1x,a,1x,a)') ' The HOST is ', trim(HostName) 
   !write(unit=iUnit,fmt='("#",1x,a,1x,a)') ' Units       = ' , trim(unit) 
   !write(unit=iUnit,fmt='("#",1x,a,1x,a)') ' First Column    = ' , ' Time [decy] ' 
   !write(unit=iUnit,fmt='("#",1x,a,1x,a)') ' Second Column   = ' , ' x1 ' 
   !write(unit=iUnit,fmt='("#",1x,a,1x,a)') ' Third Column    = ' , ' x2 '  
   !write(unit=iUnit,fmt='("#",1x,a,1x,a)') ' Fourth Column   = ' , ' x3 ' 
   !write(unit=iUnit,fmt='("#",1x,a,1x,a)') ' Fifth Column    = ' , ' Time [mjd] ' 
   write(unit=iUnit,fmt='("#",1x,a,1x,i10)') ' Number of Lines = ' , nElements
   write(unit=iUnit,fmt='("#",1x,a,a)') waves, waves 

 end subroutine WriteTimeSeriesHeader  

 end module HandleTimeSeries
