/**
    @file ApplyPointing.cpp
    @author Brian Magill
    @creation date 4/2/2008

	$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 Applies the pointing signal correction
*/

#include <iostream>
#include "EventVar.h"
#include "ApplyPointing.h"
#include "BoxCarSmooth.h"
#include "PointingCutoff.h"
#include "SOFIE_namespace.h"

using namespace std;

ApplyPointing::ApplyPointing(valarray<double> const & az, valarray<double> const & el, 
				valarray<double> const & extent, PointingDriftScaling const & scale, bool flag, double cutoff)
    {
    lockdownAz.resize(az.size());
    lockdownEl.resize(el.size());
    solarExtent.resize(extent.size());
    lockdownAz = az;
    lockdownEl = el;
    solarExtent = extent;
    calcScaling = scale;
    mode = flag;
    cutoffTime = cutoff;
    }

ChannelSignals ApplyPointing::operator()(valarray<double> const &timeArray, 
                                            ChannelSignals const &in, PointingSigParam const &pointingParam)
    {
    const int GOOD = 0;

    double scaleFactor;
	long cutoffIndex = 0;
    int gasId = in.getID();
	double wave[16] = {0.292, 0.33, 0.867, 1.037, 2.462, 2.618, 2.785, 2.939, 3.064, 3.186, 3.384, 3.479, 4.324, 4.646, 5.006, 5.316};
    int strongBands[] = {1, 3, 6, 7,  9, 11, 13, 16};
    int weakBands[]   = {2, 4, 5, 8, 10, 12, 14, 15};
    valarray<double> strong = in.getSignal(SOFIE::Strong);
    valarray<double> weak = in.getSignal(SOFIE::Weak);
    valarray<double> diff = in.getSignal(SOFIE::Diff);  
    valarray<double> smooth_model;  
    double strong_atten = in.getAttenuation(SOFIE::Strong); 
    double weak_atten   = in.getAttenuation(SOFIE::Weak);

    double A; 
    double B; 
    double C; 
	double u = 0.0;
	double v = 0.0;
    BoxCarSmooth smooth(10);
    PointingCutoff modelCutoff(cutoffTime, mode);
	std::cout << "\tIn ApplyPointing" << std::endl;
      
    if(pointingParam.getQuality(SOFIE::Strong) == GOOD)
        {
		std::cout << "\t\tProcessing Strong with SSF" << std::endl;
        A = pointingParam.getAzimCoef(SOFIE::Strong);
        B = pointingParam.getElevCoef(SOFIE::Strong); 
        C = pointingParam.getOffset(SOFIE::Strong);    
        valarray<double> strongModel =  (A*lockdownAz + B*lockdownEl + C);   
        smooth(strongModel, smooth_model);
        smooth_model = modelCutoff(timeArray, 1.0/smooth_model, cutoffIndex);
        scaleFactor = strongModel[cutoffIndex];
        strong = scaleFactor*strong*smooth_model;
        }
	else
	{
		int band = strongBands[gasId] - 1;
		std::cout << "\t\tProcessing Strong with Allen Curve for band " << band << std::endl;
		bool gotUV = getAllenUVFromWavelength(wave[band], u, v);
		if (gotUV)
		{
			valarray<double> strongModel(lockdownAz.size());
			double avgSolarRadius = getAverageSolarRadius(solarExtent);
			for (unsigned int i = 0; i < lockdownAz.size(); i++)
			{
				if (solarExtent[i] > -1.0e23)
				{
					strongModel[i] = getIntensityFromAllenCurve(lockdownAz[i], lockdownEl[i], avgSolarRadius, u, v);
				}
				else
				{
					strongModel[i] = -1.0e24;
				}
			}
			smooth(strongModel, smooth_model);
			smooth_model = modelCutoff(timeArray, 1.0/smooth_model, cutoffIndex);
			scaleFactor = strongModel[cutoffIndex];
			strong = scaleFactor*strong*smooth_model;
		}
	}

    if(pointingParam.getQuality(SOFIE::Strong) == GOOD)
	{        
		std::cout << "\t\tProcessing Weak with SSF" << std::endl;
        A = pointingParam.getAzimCoef(SOFIE::Weak);
        B = pointingParam.getElevCoef(SOFIE::Weak);  
        C = pointingParam.getOffset(SOFIE::Weak);
        valarray<double> weakModel =  (A*lockdownAz + B*lockdownEl + C);   
        smooth(weakModel, smooth_model);  
        smooth_model = modelCutoff(timeArray, 1.0/smooth_model, cutoffIndex);
        scaleFactor = weakModel[cutoffIndex];
        weak = scaleFactor*weak*smooth_model;
	}
	else
	{
		std::cout << "\t\tProcessing Weak with Allen Curve" << std::endl;
		int band = weakBands[gasId] - 1;
		bool gotUV = getAllenUVFromWavelength(wave[band], u, v);
		if (gotUV)
		{
			valarray<double> weakModel(lockdownAz.size());
			double avgSolarRadius = getAverageSolarRadius(solarExtent);
			for (unsigned int i = 0; i < lockdownAz.size(); i++)
			{
				if (solarExtent[i] > -1.0e23)
				{
					weakModel[i] = getIntensityFromAllenCurve(lockdownAz[i], lockdownEl[i], avgSolarRadius, u, v);
				}
				else
				{
					weakModel[i] = -1.0e24;
				}
			}
			smooth(weakModel, smooth_model);
			smooth_model = modelCutoff(timeArray, 1.0/smooth_model, cutoffIndex);
			scaleFactor = weakModel[cutoffIndex];
			weak = scaleFactor*weak*smooth_model;
		}
	}

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

bool ApplyPointing::getAllenUVFromWavelength(double wavelength, double &u, double &v)
{
	double wave[15] = {0.26, 0.32,0.35,0.37,0.40,0.45,0.5,0.6,0.8,1.0,1.5,2.0,3.0,5.0,10.05};
	double U[15]    = {0.66, 0.86,0.96,1.01,0.93,0.99,0.95,0.86,0.71,0.63,0.56,0.48,0.36,0.23,0.14};
	double V[15]    = {0.31, 0.05,-0.08,-0.15,-0.06,-0.17,-0.20,-0.21,-0.20,-0.19,-0.21,-0.18,-0.13,-0.08,-0.04};

	if ((wavelength > wave[14]) or (wavelength < wave[0]))
	{
		std::cerr << "Wavelength for Allen Curve out of bounds in ApplyPointing::getAllenUVFromWavelength()!" << std::endl; 
		return false;
	}

	for (int i = 0; i < 15; i++)
	{
		if (wave[i] >= wavelength)
		{
			u = (U[i] * (wavelength - wave[i-1]) + U[i-1] * (wave[i] - wavelength)) / (wave[i] - wave[i-1]);
			v = (V[i] * (wavelength - wave[i-1]) + V[i-1] * (wave[i] - wavelength)) / (wave[i] - wave[i-1]);
			return true;
		}
	}
	std::cerr << "Wavelength for Allen Curve out of bounds in ApplyPointing::getAllenUVFromWavelength()!" << std::endl; 
	return false;
}

double ApplyPointing::getAverageSolarRadius(valarray<double> const &extent)
{
	unsigned int count = 0;
	unsigned int samples = 20;
	unsigned int midPoint = 0;
	double sum = 0.0;
	double value = 0.0;

	midPoint = int(extent.size() / 2);

	if (midPoint < samples)
	{
		std::cerr << "Not enough solar extent samples in ApplyPointing::getAverageSolarRadius!" << std::endl;
		return 0.0;
	}

	while (count < samples)
	{
		// Must divide extent by 2 to get radius
		value =  (extent[midPoint - count]/2);
		sum += value;
		count++;
		//std::cout << "\t\t\t\tCurrent radius = " << value << "\tCurrent sum = " << sum << std::endl;
	}
	//std::cout << "\t\t\t\t\tAverage Solar radius = " << (sum/samples) << std::endl;
	return sum/samples;
}

double ApplyPointing::getIntensityFromAllenCurve(double az, double el, double radius, double u, double v)
{
	if (radius == 0.0)
	{
		std::cerr << "Extent is zero.  Cowardly refusing to divide by zero in ApplyPointing::getIntensityFromAllenCurve!" << std::endl;
		return -1.0e24;
	}

	double normAz = az/radius;
	double normEl = 1 - el/radius;
	double normRadiusSquare = (normAz * normAz) + (normEl * normEl);

	if (!isfinite(normAz) || !isfinite(normEl) || !isfinite(normRadiusSquare) || normRadiusSquare > 1.0)
	{
		if(!isfinite(normAz))
		{
			std::cerr << "\t\t\t\t\tNormalized Az = " << normAz << "\tAz = " << az << "\tradius = " << radius << std::endl;
		}
		if(!isfinite(normEl))
		{
			std::cerr << "\t\t\t\t\tNormalized El = " << normEl << "\tEl = " << el << "\tradius = " << radius << std::endl;
		}
		if(!isfinite(normRadiusSquare))
		{
			std::cerr << "\tNormalized Radius^2 = " << normRadiusSquare << "\tEl = " << el << "\tAz = " << az << "\tradius = " << radius << std::endl;
		}
		return -1.0e24;
	}

	return (1 - u + u * sqrt(1 - normRadiusSquare) - v * (normRadiusSquare));
}