/**
    @file CalcDrift.cpp
    @author Brian Magill
    @creation date 1/3/2007

	$Date:$
    $Revision:$

    @copyright
     (©) 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.
        
    @brief Calculates the drift signal for a given channel
*/
#include <iostream>
#include <stdexcept>
#include "CalcDrift.h"
#include "DetermDriftCorr.h"
#include "DetermDiffDrift.h"

#include "SOFIE_namespace.h"

using namespace std;

void CalcDrift::operator()(valarray<double> const &time, ChannelSignals const &in, 
                        DriftParam & param, std::string &comments) const 
{

//  Set parameters for strong, weak, and difference signals to values that won't alter the signals

unsigned long vecSize = 3;
vector<double> slope(vecSize);
vector<double> stdev(vecSize);
vector<double> sigRef(vecSize);
vector<int> quality(vecSize);
double delta_time;
double Va;
double Ka;

valarray<double> signal;
DetermDriftCorr detDrift;
DetermDiffDrift diffDrift;
bool correctionValid =  true;

if (end_time - start_time  < minInterval)
    {
    correctionValid =  false;
    end_time = minInterval + start_time;
    }
    
unsigned long index  = start_index(time);
unsigned long length = end_index(time) - index + 1;

valarray<double> time_intv = time[slice(index, length, 1)];

in.getSignal(SOFIE::Strong, signal);
valarray<double> strong_intv = signal[slice(index, length, 1)];
in.getSignal(SOFIE::Weak, signal);                               
valarray<double> weak_intv = signal[slice(index, length, 1)];
in.getSignal(SOFIE::Diff, signal);
valarray<double> diff_intv = signal[slice(index, length, 1)];

if (correctionValid)
    {
    
    detDrift = DetermDriftCorr(time_intv, strong_intv);
    slope[SOFIE::Strong]  = detDrift.getSlope();
    stdev[SOFIE::Strong]  = detDrift.getSlopeError();
    delta_time = detDrift.getTimeRef() - ref_time;
    Ka = detDrift.getSlope();
    Va = detDrift.getVref();
    sigRef[SOFIE::Strong] = Va*(1. + Ka*delta_time);
// David's Code Does this for Band 16
//  std::cout << " band " << in.getID() << std::endl;
    if(in.getID() == SOFIE::NO) {
//std::cout << " band 16 " << std::endl; exit(23);
      slope[SOFIE::Strong] = 0.0;
      stdev[SOFIE::Strong] = 0.0;
      sigRef[SOFIE::Strong] = Va;
    }

    detDrift = DetermDriftCorr(time_intv, weak_intv);
    slope[SOFIE::Weak]  = detDrift.getSlope();
    stdev[SOFIE::Weak]  = detDrift.getSlopeError();
    delta_time = detDrift.getTimeRef() - ref_time;
    Ka = detDrift.getSlope();
    Va = detDrift.getVref();
    sigRef[SOFIE::Weak] = Va*(1. + Ka*delta_time);
    
    diffDrift = DetermDiffDrift(time_intv, diff_intv);
    slope[SOFIE::Diff]  = diffDrift.getSlope();
    stdev[SOFIE::Diff]  = diffDrift.getSlopeError();
    delta_time = detDrift.getTimeRef() - ref_time;    
    Ka = diffDrift.getSlope();
    Va = diffDrift.getVref();
    sigRef[SOFIE::Diff] = Va - Ka*delta_time;
    
    quality = vector<int> (vecSize, 0);
    }
else
    {
    
    detDrift = DetermDriftCorr(time_intv, strong_intv);
    sigRef[SOFIE::Strong] = detDrift.getVref();

    detDrift = DetermDriftCorr(time_intv, weak_intv);
    sigRef[SOFIE::Weak] = detDrift.getVref();

    diffDrift = DetermDiffDrift(time_intv, diff_intv);
    sigRef[SOFIE::Diff] = diffDrift.getVref();

    slope   = vector<double> (vecSize, 0.0);
    stdev   = vector<double> (vecSize, 0.0);
    quality = vector<int> (vecSize, 1);
    
    }

double interval = end_time - start_time; 
double ref_timeA = ref_time;
param = DriftParam(in.getID(), slope, stdev, sigRef, quality, ref_timeA, interval);
  
comments = "";

};

unsigned long CalcDrift::start_index(valarray<double> const &time) const
{
    unsigned long index = 0;
    bool foundIndex = false;
    
    for(unsigned long i = 0; i < time.size(); i++)
        if (time[i] >= start_time)
        {
        index = i;
	foundIndex = true;
        break;
        }
    if (!foundIndex)
	throw runtime_error("CalcDrift::start_index  : Could not find index.  All times less than beginning of calculation range");

return index;
}

unsigned long CalcDrift::end_index(valarray<double> const &time) const
{
    unsigned long index = 0;
    bool foundIndex = false;

    for(unsigned long i = 0; i < time.size(); i++)
        if (time[i] >= end_time)
        {
        index = i;
	foundIndex = true;	
        break;
        }

if (!foundIndex)
	throw runtime_error("CalcDrift::end_index  : Could not find index. All times less than end of calculation range");
		 
return index;
}