/**
     @file FOV_Offset.cpp
     @author Lance Deaver

	$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 FOV Offsets between band 1-16 and band 3.
     
*/

#include "FOV_Offset.h"

#include "F77Functions.h"
#include "SigEarthLocParam.h"
#include "SOFIE_namespace.h"
#include "GATS_Utilities.hpp"


using namespace GATS_Utilities;

FOV_Offset::FOV_Offset(Event& L0, ConfigFile& cf)
{
	offsets = std::vector<double>(32,0.0);
    if(cf.ValidEntry("FOV_Offsets", "FOV_Offsets") ) {
		std::vector<float> vals = cf.GetRealValues("FOV_Offsets", "FOV_Offsets");
        assert(vals.size() == 2* 16);
		std::copy(vals.begin(), vals.end(), offsets.begin()) ;
    } 
//	std::cout << " string " << L0.getSRSSFlag() << std::endl;
	mode = ('s' == L0.getSRSSFlag() )  ? 0 : 1 ;
}

//-----------------------------------------------------------------------------------------
//  Object method that sets up the FOV Offset correction 

//  Input:
//     inSignal........Input signals as multiple signals (All channels)
//
//  Output:
//     outSignal.......The corrected signals as multiple signals
//
//  Source: Lance Deaver, GATS
//-----------------------------------------------------------------------------------------
MultiSignal FOV_Offset::correct(MultiSignal &inSignal, SigEarthLocParam &earthInfo ) 
{
	std::cout << "Inside FOV_Correction" << std::endl;


    valarray<double> timeArray;
    vector<ChannelSignals> signalArray;
    Journal auditJournal= inSignal.getAudit();
    string comments;  
    string processName = "FOV Offset Correction";
    ChannelSignals Channel;


    try {
		inSignal.getTime(timeArray);
	    std::vector<double> ViewAngles = earthInfo.getViewAngle();
/*
		std::vector<double> alts = earthInfo.getAltitudes() ; 
for(int i = 0 ; i< alts.size(); ++i) 
	std::cout << ViewAngles[i]*180./M_PI << "  " << alts[i] << std::endl;
exit(23);
*/
		std::vector<double> voff(ViewAngles.size() );
		std::valarray<double> S(ViewAngles.size() );
		int N = (int)ViewAngles.size();

for(int i = 0; i<= SOFIE::NO ; ++i) {
		Channel = inSignal.getChannel(i);
		std::valarray<double> sigs =  Channel.getSignal(SOFIE::Strong) ;
        std::valarray<double> sigw =  Channel.getSignal(SOFIE::Weak) ;
		assert (sigs.size() == ViewAngles.size() && sigw.size() == ViewAngles.size()  );

		int istrong = (4*i) + mode ;
    	int iweak = istrong+2;
	
	    if (i == SOFIE::NO || i == SOFIE::H2O )  swap(istrong,iweak);

//std::cout << "Offset Strong is " << offsets[istrong] << std::endl;

		std::transform(ViewAngles.begin(), ViewAngles.end(), voff.begin(),
     		std::bind2nd( std::minus<double>(),offsets.at(istrong) *  M_PI/(180*60)   ) ); 
		::linerd_( &voff[0],  &sigs[0], &ViewAngles[0], &S[0], &N, &N );
		sigs = S;

		std::transform(ViewAngles.begin(), ViewAngles.end(), voff.begin(),
     		std::bind2nd( std::minus<double>(),offsets.at(iweak) *  M_PI/(180*60)   ) );
		::linerd_( &voff[0], &sigw[0],  &ViewAngles[0], &S[0], &N, &N   );
		sigw = S;

 		signalArray.push_back( ChannelSignals(Channel.getID(),
                 Channel.getAttenuation(SOFIE::Strong),
                 Channel.getAttenuation(SOFIE::Weak),
            sigs,//      Channel.getSignal(SOFIE::Strong),
            sigw, //      Channel.getSignal(SOFIE::Weak),
                  Channel.getSignal(SOFIE::Diff)
        ));
}
/*
std::valarray<double> S = Channel.getSignal(SOFIE::Strong) ;  
for(int i =0; i< tpa_alts.size() ; ++i) 
	std::cout <<  ViewAngles[i] << "  " <<  tpa_alts[i] <<  std::endl;
std::cout << "&" << std::endl << "#next " << std::endl;

for(int i =0; i< tpa_alts.size() ; ++i)
    std::cout <<  S[i] << "  " <<  tpa_alts[i] <<  std::endl;
std::cout << "&" << std::endl ;
*/

/*
	
        Channel = inSignal.getChannel(SOFIE::NO);   //NO channel
        valarray<double> strong = Channel.getSignal(SOFIE::Strong);   // get the strong band
		valarray<double> diffsig = Channel.getSignal(SOFIE::Diff); 

		int npts=(int)strong.size() ;  // number of points
        int iStart=0;  // index of starting location


 if( inSignal.getModeFlag() ) {  // sunset
        double *t_ptr = std::find_if( &timeArray[0], &timeArray[timeArray.size()],
            std::bind2nd(std::greater_equal<double>(), BalanceStopTime  ) ) ;
		iStart = std::distance(&timeArray[0], t_ptr);
        npts = timeArray.size() - iStart ;

 } else {   // sunrise
		double *t_ptr = std::find_if( &timeArray[0], &timeArray[timeArray.size()],
			std::bind2nd(std::greater<double>(), BalanceStartTime  ) ) ;
        npts = std::distance(&timeArray[0], t_ptr);
        iStart=0;
}
*/


// std::cout << *t_ptr << "  " << EventStopTime << "  " << timeArray[npts-1] << std::endl;
//std::cout << " npts " << npts << " size " << strong.size() << std::endl;
//exit(23);
//		int npts = strong.size();
/*
		double maxBump= -0.003 , weight=0.0 ; 
        double minFitz=140., maxFitz=190., startBA = BalanceStartTime,
			endBA=BalanceStopTime; 
		char mode= (char)inSignal.getModeFlag(), band=(char)16;
	
		int chisFlag, bumpFlag;
	
		valarray<double> residual = strong , maxUncertianty(npts), 
	    	 minUncertianty(npts), Aparameters(9), Sigmaparameters(9) ;
*/
	
// OK now call David's correction Code

//cout << " before davids code "   << npts << endl;

/*
	::sofieNO8_natural (
		&strong[iStart], 
		&TP_Alts[iStart],
		&timeArray[iStart],
		&mode,
		&residual[iStart],
//		&band,
		&npts,
		&resChi[0],
		&maxChi,
		&chisFlag,
		&Aparameters[0],
		&Sigmaparameters[0],
		&maxUncertianty[0],
		&minUncertianty[0],
		&bumpFlag,
		&maxBump,
		&weight,
      &minFitz,
      &maxFitz,
      &startBA,
      &endBA );

if(chisFlag) {   // Bad fit

 cerr << " **NO_Thermal Correction Not Performed**" << endl 
     << "Chi Squared Value " << ConvertToString( *resChi) << " Exceeds Limit: " << ConvertToString(maxChi) ; 
inSignal.getAudit().print(cerr);
cerr << endl;

resChi[0] *=  -1;

return inSignal; 
}

//std::cout << " before get sig " << std::endl;
valarray<double> r1 =  diffsig ; 
    ::sofieNO8_natural (
        &diffsig[iStart],
        &TP_Alts[iStart],
        &timeArray[iStart],
        &mode,
        &r1[iStart],
//      &band,
        &npts,
        &resChi[1],
        &maxChi,
        &chisFlag,
        &Aparameters[0],
        &Sigmaparameters[0],
        &maxUncertianty[0],
        &minUncertianty[0],
        &bumpFlag,
        &maxBump,
        &weight,
      &minFitz,
      &maxFitz,
      &startBA,
      &endBA );
if(chisFlag) {   // Bad fit

resChi[1] *= -1.;
r1 = diffsig;
}

//std::cout << "after getsig " << std::endl;





cout << "This is the number " <<  signalArray.size()  <<  " " << Channel.getID() << endl;

        signalArray.push_back( ChannelSignals(Channel.getID(), 
	             Channel.getAttenuation(SOFIE::Strong),
                 Channel.getAttenuation(SOFIE::Weak),
	             residual, 
                  Channel.getSignal(SOFIE::Weak),
	              r1 //  Channel.getSignal(SOFIE::Diff)  
				)
	              );
		
//cout << "This is the nubmer " <<  signalArray.size()  << endl;	exit(23);


	for(int i = SOFIE::NO+1 ; i< inSignal.numberOfSignals() ; ++i) {
	    signalArray.push_back( inSignal.getChannel(i) );
	}
		

for(int i =0; i< residual.size(); ++i) 
    cout << ConvertToStringPrec(strong[i]) << "  " << ConvertToStringPrec(TP_Alts[i])  << endl;
cout << "&" << endl;

cout << "# " << ConvertToStringPrec( *resChi) << std::endl;
for(int i =0; i< residual.size(); ++i) 
    cout << ConvertToStringPrec(residual[i]) << "  " << ConvertToStringPrec(TP_Alts[i])  << endl;

exit(23);
*/    

    } catch (exception &ex) {
            
            cerr << "Exception thrown in FOV_Offsets :" << endl;
            cerr << "Channel Gas ID: " << Channel.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);
    
     cout << "finished FOV_Offset correction" << endl;
     return MultiSignal(inSignal.getEventNumber(), inSignal.get150kmTime(), 
                             inSignal.getModeFlag(), timeArray, signalArray, auditJournal) ;

}