/**
     @file NonLinearity.cpp
     @author Greg Paxton
     @date 4/08/2009

	$Date:$
    $Revision:$

    @copyright
     (©) Copyright 2009 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 Corrects the nonlinearity in the signals
     
*/

#include <cmath>
#include <iostream>
#include <vector>
#include <stdexcept>
#include "SOFIE_namespace.h"
#include "GATS_Utilities.hpp"
#include "NonLinearity.h"
#include "ChannelSignals.h"

using namespace std;
using namespace GATS_Utilities;

NonLinearity::NonLinearity(ConfigFile& cf)
{
	cout << "Before getting values" << endl;
	Cnl   = cf.GetRealValues("NonLinearityCorrection","Cnl");
	GAcal = cf.GetReal("NonLinearityCorrection","GAcal");
	cout << "After getting values" << endl;
}

//-----------------------------------------------------------------------------------------
//  Object method that sets up the nonlinearity correction for all channels.
//      1. Extracts a channel from the input MultiSignal object
//      2. Passes each channel into the correction function
//
//  Input:
//     inSignal........Input signals as multiple signals (All channels)
//
//  Output:
//     outSignal.......The corrected signals as multiple signals
//
//  Source: Greg Paxton, GATS
//-----------------------------------------------------------------------------------------
MultiSignal NonLinearity::correctAllSignals(MultiSignal &inSignal) 
{
    
    valarray<double> timeArray;
    vector<ChannelSignals> signalArray;

    Journal auditJournal;
      
    string processName = "Nonlinearity Correction";
    ChannelSignals inChannel;
    ChannelSignals outChannel;
    string comments;

    inSignal.getTime(timeArray);    

    auditJournal = inSignal.getAudit();

    for(unsigned int indx = 0; indx < inSignal.numberOfSignals(); indx++)
    {
        try {
                
            inChannel = inSignal.getChannel(indx);        
            outChannel = correction(inChannel);
            signalArray.push_back(outChannel);   

            } catch (exception &ex) {
            
            cerr << "Exception thrown in NonLinearity :" << endl;
            cerr << "Channel Gas ID: " << inChannel.getID() << endl;
            cerr << "Exception Message: " << ex.what() << endl;
            cerr << "Additional comments: " << comments << endl;
            cerr << "Current Audit Trail: ";
            inSignal.getAudit().print(cerr);
            cerr << endl;
            
            throw;
            }
    }

    auditJournal.append(processName);
    
    MultiSignal outSignal(inSignal.getEventNumber(), inSignal.get150kmTime(), 
                            inSignal.getModeFlag(), timeArray, signalArray, auditJournal) ;
                
    return outSignal;
}

//-----------------------------------------------------------------------------------------
//  Object method that sets up the nonlinearity correction for one channel.
//		1. It retrieves the channel signals
//		2. Retrieves the attenuation for the strong and weak bands
//		3. Passes each band's signals into a smoothing function
//		4. Corrects the signal nonlinearity by using a ratio of the corrected response
//			and smoothed signal to reduce noise added in via the nonlinearity correction.
//
//  Input:
//     in..............Input signal as channels (Strong and Weak bands)
//
//  Output:
//     ChannelSignals..The corrected signals as channels
//
//  Source: Greg Paxton, GATS
//-----------------------------------------------------------------------------------------
ChannelSignals NonLinearity::correction(ChannelSignals &in) 
{
    int gasId = in.getID();
    double strong_atten = in.getAttenuation(SOFIE::Strong);
    double weak_atten   = in.getAttenuation(SOFIE::Weak);
	double smoothedResponse = 0.0;
    int strongBands[] = {1, 3, 6, 7,  9, 11, 13, 16};
    int weakBands[]   = {2, 4, 5, 8, 10, 12, 14, 15};

	cout << "Performing NonLinearity correction for channel " << gasId << endl;

    std::valarray<double> strong = in.getSignal(SOFIE::Strong);
    std::valarray<double> weak = in.getSignal(SOFIE::Weak);
    std::valarray<double> diff = in.getSignal(SOFIE::Diff);

	std::valarray<double> strongSmoothed = cosBell(strong, 4);
	std::valarray<double> weakSmoothed = cosBell(weak, 4);

    for (unsigned int i=0; i < strong.size(); i++)
    {
        //strong[i] = sofie_nonlinearity(strongBands[gasId], 1, strong[i], strong_atten);
        smoothedResponse = sofie_nonlinearity(strongBands[gasId], 1, strongSmoothed[i], strong_atten);
		strong[i] = strong[i] * (smoothedResponse/strongSmoothed[i]);
    }
    for (unsigned int i=0; i < weak.size(); i++)
    {
        //weak[i] = sofie_nonlinearity(weakBands[gasId], 1, weak[i], weak_atten);
        smoothedResponse = sofie_nonlinearity(weakBands[gasId], 1, weakSmoothed[i], weak_atten);
		weak[i] = weak[i] * (smoothedResponse/weakSmoothed[i]);
    }

    return ChannelSignals(gasId, strong_atten, weak_atten, strong, weak, diff);

}

//-----------------------------------------------------------------------------------------
//  Object method smooths the input signal using a cosine bell function
//
//  Input:
//     in..............Input signal as valarray doubles
//     halfWidth.......The half width of the cosine bell function. (Note:  The outside
//						points will always be zero weighted.)
//
//  Output:
//     smoothed...normal exit: The smoothed data
//                error  exit: The input array with no smoothing applied
//
//  Source: Greg Paxton, GATS
//-----------------------------------------------------------------------------------------
std::valarray<double> NonLinearity::cosBell(std::valarray<double> &in, int halfWidth) 
{
	unsigned int points = 2 * halfWidth + 1;
	double sumSi = 0.0;
	const double PI = 3.14159265;
	std::valarray<double> SiNorm(points);
	std::valarray<double> smoothed(in.size());

    if(points > in.size())
    {
		cerr << "In NonLinearity::cosBell the number of points to smooth was larger than the array passed in!" << endl;
		cerr << "\t\tPassing back the input array without applying the smoothing!" << endl;
		return in;
    }

	// Determine the sum of Si for use in normalizing the
	// weighting function
	for(int i = (-halfWidth); i < ((int)points - halfWidth); i++)
	{
		sumSi += 1 + cos((PI/halfWidth)*i);
	}

	// Determine the normalized cosBell weighting for each point
	for(int i = (-halfWidth); i < ((int)points - halfWidth); i++)
	{
		SiNorm[i + halfWidth] = (1 + cos((PI/halfWidth)*i))/sumSi;
	}

	// Smooth the data using the cosBell weighting function
    for(unsigned int i = halfWidth; i < (in.size() - halfWidth); i++)
    {
		double sum = 0.0;
		for(int j = (-halfWidth); j < ((int)points - halfWidth); j++)
		{
			sum += SiNorm[j + halfWidth]  * in[i+j];
		}
		smoothed[i] = sum;
    }

	// Simpy use the first smoothed point for all points prior
	for (int i = 0; i < halfWidth; i++)
	{
		smoothed[i] = smoothed[halfWidth];
	}

	// Simpy use the last smoothed point for all points post
	for (unsigned int i = (in.size() - halfWidth); i < in.size(); i++)
	{
		smoothed[i] = smoothed[((in.size() - halfWidth) - 1)];
	}

	return smoothed;
}

//-----------------------------------------------------------------------------------------
//  Object method that deals with the nonlinear response of SOFIE radiometer measurements.
//
//  There are 2 options:
//
//    1) Correct (remove) nonlinearity in the input radiometer signal.
//
//    2) Make the input linear signal be nonlinear.
//
//  Input:
//     band.......band number (integer from 1-16)
//     resp.......the input response (counts), can range from 0 - 32768
//     GA.........balance attenuator setting (gain) for the measurement (can range from 0-1)
//     option.....if 1 then correct nonlinearity in resp, if 2 then make resp nonlinear
//
//  Output:
//     resp_out...option 1: response corrected for nonlinearity
//                option 2: response with nonlinearity applied 
//
//  Source: Mark Hervig, GATS
//-----------------------------------------------------------------------------------------
double NonLinearity::sofie_nonlinearity(int band, int option, double resp, double GA)
{

// check the input for problems

assert(band >= 1 && band <= 16);
/*
   if (band < 1 || band > 16) {
     cout << "band number out of range in sofie_nonlinearity.cpp:  " << band << endl;
     return resp;
   }
*/
//assert( resp >= 0. && resp <= 32768.0);
//   if (resp < 0 || resp > 32768.0) {
//     cout << "resp out of range in sofie_nonlinearity.cpp:  " << resp << endl;
//     return resp;
//   }

assert( GA >= 0. && GA <= 1.);
//   if (GA < 0.0 || GA > 1.0) {
//     cout << "ballance attenuator setting out of range in sofie_nonlinearity.cpp:  " << GA << endl;
//     return resp;
//   }

assert(option >= 0 || option <= 2);
//   if (option < 1 || option > 2) {
//     cout << "option is out of range in sofie_nonlinearity.cpp:  " << option << endl;
//     return resp;
//   }

// Cnl are the nonlinearity constants derived in ground calibration,
// see sofiedata.org for further explanation.  The balance attenuator
// setting (electronic gain, GA) for these calibration results was GAcal = 0.83

/* Removed the following lines so that the data values can be read from cfg
   static const double GAcal = 0.83;
   static const double Cnl[16] =   // the nonlinearity constants,  1/counts
                {
                    0.,  // band 1
                    0.,
                    0.,
                    0.,
                    1.79e-6, // band 5
                    1.56e-6,
                    9.58e-6,
                    8.55e-6,
                    0.80e-6, // band 9
                    1.60e-6,
                    1.74e-6,
                    2.56e-6,
                    5.01e-6, // band 13
                    3.15e-6,
                    1.93e-6,
                    2.20e-6 // band 16
                };
*/

// correct nonlinearity or apply nonlinearity

   double resp_out = 0.0;  // the output response

   if (option == 0) {
     resp_out = resp;  // no correction
   }

   if (option == 1) {
     resp_out = resp / (1.0 -  Cnl[band-1] * resp * GAcal / GA);  // correct nonlinearity
   }

   if (option == 2) {
     resp_out = resp / (1.0 +  Cnl[band-1] * resp * GAcal / GA);  // apply nonlinearity
   }

   return resp_out;
}