// Copyright (c) 2019, Christoph Peters
// All rights reserved.
// 
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//     * Redistributions of source code must retain the above copyright
//       notice, this list of conditions and the following disclaimer.
//     * Redistributions in binary form must reproduce the above copyright
//       notice, this list of conditions and the following disclaimer in the
//       documentation and/or other materials provided with the distribution.
//     * Neither the name of the Karlsruhe Institute of Technology nor the
//       names of its contributors may be used to endorse or promote products
//       derived from this software without specific prior written permission.
// 
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY
// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


#pragma once
#include <cmath>
#include "complex_algebra.hpp"


void trigonometricToExponentialMomentsReal3(float_complex pOutExponentialMoment[3],const float pTrigonometricMoment[3]){
	float zerothMomentPhase=3.14159265f*pTrigonometricMoment[0]-1.57079633f;
	pOutExponentialMoment[0]=float_complex(cosf(zerothMomentPhase),sinf(zerothMomentPhase));
	pOutExponentialMoment[0]=0.0795774715f*pOutExponentialMoment[0];
	pOutExponentialMoment[1]=pTrigonometricMoment[1]*float_complex(0.0f,6.28318531f)*pOutExponentialMoment[0];
	pOutExponentialMoment[2]=pTrigonometricMoment[2]*float_complex(0.0f,6.28318531f)*pOutExponentialMoment[0]+pTrigonometricMoment[1]*float_complex(0.0f,3.14159265f)*pOutExponentialMoment[1];
	pOutExponentialMoment[0]=2.0f*pOutExponentialMoment[0];
}

void trigonometricToExponentialMoments2(float_complex pOutExponentialMoment[2],const float_complex pTrigonometricMoment[2]){
	float zerothMomentPhase=3.14159265f*pTrigonometricMoment[0].x-1.57079633f;
	pOutExponentialMoment[0]=float_complex(cosf(zerothMomentPhase),sinf(zerothMomentPhase));
	pOutExponentialMoment[0]=0.0795774715f*pOutExponentialMoment[0];
	pOutExponentialMoment[1]=float_complex(0.0f,6.28318531f)*pOutExponentialMoment[0]*pTrigonometricMoment[1];
	pOutExponentialMoment[0]=2.0f*pOutExponentialMoment[0];
}


void levinsonsAlgorithm2(float_complex pOutSolution[2],const float_complex pFirstColumn[2]){
	pOutSolution[0]=float_complex(1.0f/(pFirstColumn[0].x),0.0f);
	float_complex dotProduct;
	float_complex flippedSolution1;
	float factor;
	dotProduct=pOutSolution[0].x*pFirstColumn[1];
	factor=1.0f/(1.0f-absSqr(dotProduct));
	flippedSolution1=float_complex(pOutSolution[0].x,0.0f);
	pOutSolution[0]=float_complex(factor*pOutSolution[0].x,0.0f);
	pOutSolution[1]=factor*(-flippedSolution1.x*dotProduct);
}

void levinsonsAlgorithm3(float_complex pOutSolution[3],const float_complex pFirstColumn[3]){
	pOutSolution[0]=float_complex(1.0f/(pFirstColumn[0].x),0.0f);
	float_complex dotProduct;
	float_complex flippedSolution1;
	float_complex flippedSolution2;
	float factor;
	dotProduct=pOutSolution[0].x*pFirstColumn[1];
	factor=1.0f/(1.0f-absSqr(dotProduct));
	flippedSolution1=float_complex(pOutSolution[0].x,0.0f);
	pOutSolution[0]=float_complex(factor*pOutSolution[0].x,0.0f);
	pOutSolution[1]=factor*(-flippedSolution1.x*dotProduct);
	dotProduct=pOutSolution[0].x*pFirstColumn[2]+pOutSolution[1]*pFirstColumn[1];
	factor=1.0f/(1.0f-absSqr(dotProduct));
	flippedSolution1=conjugate(pOutSolution[1]);
	flippedSolution2=float_complex(pOutSolution[0].x,0.0f);
	pOutSolution[0]=float_complex(factor*pOutSolution[0].x,0.0f);
	pOutSolution[1]=factor*(pOutSolution[1]-flippedSolution1*dotProduct);
	pOutSolution[2]=factor*(-flippedSolution2.x*dotProduct);
}


void computeAutocorrelation2(float_complex pOutAutocorrelation[2],const float_complex pSignal[2]){
	pOutAutocorrelation[0]=pSignal[0]*conjugate(pSignal[0])+pSignal[1]*conjugate(pSignal[1]);
	pOutAutocorrelation[1]=pSignal[0]*conjugate(pSignal[1]);
}

void computeAutocorrelation3(float_complex pOutAutocorrelation[3],const float_complex pSignal[3]){
	pOutAutocorrelation[0]=pSignal[0]*conjugate(pSignal[0])+pSignal[1]*conjugate(pSignal[1])+pSignal[2]*conjugate(pSignal[2]);
	pOutAutocorrelation[1]=pSignal[0]*conjugate(pSignal[1])+pSignal[1]*conjugate(pSignal[2]);
	pOutAutocorrelation[2]=pSignal[0]*conjugate(pSignal[2]);
}


void computeCorrelation2(float_complex pOutCorrelation[2],const float_complex pLHS[2],const float_complex pRHS[2]){
	pOutCorrelation[0]=pLHS[0]*pRHS[0]+pLHS[1]*pRHS[1];
	pOutCorrelation[1]=pLHS[1]*pRHS[0];
}


void computeImaginaryCorrelation3(float pOutCorrelation[3],const float_complex pLHS[3],const float_complex pRHS[3]){
	pOutCorrelation[0]=pLHS[0].x*pRHS[0].y+pLHS[0].y*pRHS[0].x+pLHS[1].x*pRHS[1].y+pLHS[1].y*pRHS[1].x+pLHS[2].x*pRHS[2].y+pLHS[2].y*pRHS[2].x;
	pOutCorrelation[1]=pLHS[1].x*pRHS[0].y+pLHS[1].y*pRHS[0].x+pLHS[2].x*pRHS[1].y+pLHS[2].y*pRHS[1].x;
	pOutCorrelation[2]=pLHS[2].x*pRHS[0].y+pLHS[2].y*pRHS[0].x;
}


float_complex evaluateFastHerglotzTransform2(const float_complex circlePoint,const float_complex pExponentialMoment[2],const float_complex pEvaluationPolynomial[2]){
	float_complex conjCirclePoint=conjugate(circlePoint);
	float polynomial1=pEvaluationPolynomial[0].x;
	float_complex polynomial0=pEvaluationPolynomial[1]+polynomial1*conjCirclePoint;
	float_complex dotProduct=polynomial1*pExponentialMoment[1];
	return pExponentialMoment[0]+2.0f*(dotProduct)/(polynomial0);
}

float_complex evaluateFastHerglotzTransform3(const float_complex circlePoint,const float_complex pExponentialMoment[3],const float_complex pEvaluationPolynomial[3]){
	float_complex conjCirclePoint=conjugate(circlePoint);
	float polynomial2=pEvaluationPolynomial[0].x;
	float_complex polynomial1=pEvaluationPolynomial[1]+polynomial2*conjCirclePoint;
	float_complex polynomial0=pEvaluationPolynomial[2]+conjCirclePoint*polynomial1;
	float_complex dotProduct=polynomial1*pExponentialMoment[1]+polynomial2*pExponentialMoment[2];
	return pExponentialMoment[0]+2.0f*(dotProduct)/(polynomial0);
}


float evaluateFourierSeries2(const float_complex circlePoint,const float_complex pFourierCoefficient[2]){
	// Minimize sequential dependencies by computing powers of circlePoint
	// by multiplying previous powers by powers of two
	float_complex circlePointPower1=circlePoint;
	float result=pFourierCoefficient[1].x*circlePointPower1.x-pFourierCoefficient[1].y*circlePointPower1.y;
	return 2.0f*result+pFourierCoefficient[0].x;
}


float evaluateFourierSeriesReal3(const float_complex circlePoint,const float pFourierCoefficient[3]){
	// Minimize sequential dependencies by computing powers of circlePoint
	// by multiplying previous powers by powers of two
	float_complex circlePointPower1=circlePoint;
	float_complex circlePointPower2=circlePointPower1*circlePointPower1;
	float result=pFourierCoefficient[1]*circlePointPower1.x+pFourierCoefficient[2]*circlePointPower2.x;
	return 2.0f*result+pFourierCoefficient[0];
}


void prepareReflectanceSpectrum2(float_complex pOutExponentialMoment[2],float_complex pOutEvaluationPolynomial[2],const float_complex pTrigonometricMoment[2]){
	trigonometricToExponentialMoments2(pOutExponentialMoment,pTrigonometricMoment);
	levinsonsAlgorithm2(pOutEvaluationPolynomial,pOutExponentialMoment);
	pOutEvaluationPolynomial[0]=6.28318531f*pOutEvaluationPolynomial[0];
	pOutEvaluationPolynomial[1]=6.28318531f*pOutEvaluationPolynomial[1];
}

void prepareReflectanceSpectrumReal3(float_complex pOutExponentialMoment[3],float_complex pOutEvaluationPolynomial[3],const float pTrigonometricMoment[3]){
	trigonometricToExponentialMomentsReal3(pOutExponentialMoment,pTrigonometricMoment);
	levinsonsAlgorithm3(pOutEvaluationPolynomial,pOutExponentialMoment);
	pOutEvaluationPolynomial[0]=6.28318531f*pOutEvaluationPolynomial[0];
	pOutEvaluationPolynomial[1]=6.28318531f*pOutEvaluationPolynomial[1];
	pOutEvaluationPolynomial[2]=6.28318531f*pOutEvaluationPolynomial[2];
}


float evaluateReflectanceSpectrum2(const float phase,const float_complex pExponentialMoment[2],const float_complex pEvaluationPolynomial[2]){
	float_complex circlePoint;
	circlePoint=float_complex(cosf(phase),sinf(phase));
	float_complex herglotzTransform;
	herglotzTransform=evaluateFastHerglotzTransform2(circlePoint,pExponentialMoment,pEvaluationPolynomial);
	return fast_atan(herglotzTransform.y/herglotzTransform.x)*0.318309886f+0.5f;
}

float evaluateReflectanceSpectrum3(const float phase,const float_complex pExponentialMoment[3],const float_complex pEvaluationPolynomial[3]){
	float_complex circlePoint;
	circlePoint=float_complex(cosf(phase),sinf(phase));
	float_complex herglotzTransform;
	herglotzTransform=evaluateFastHerglotzTransform3(circlePoint,pExponentialMoment,pEvaluationPolynomial);
	return fast_atan(herglotzTransform.y/herglotzTransform.x)*0.318309886f+0.5f;
}


void prepareReflectanceSpectrumLagrange2(float_complex pOutLagrangeMultiplier[2],const float_complex pTrigonometricMoment[2]){
	float_complex pExponentialMoment[2];
	trigonometricToExponentialMoments2(pExponentialMoment,pTrigonometricMoment);
	float_complex pEvaluationPolynomial[2];
	levinsonsAlgorithm2(pEvaluationPolynomial,pExponentialMoment);
	pEvaluationPolynomial[0]=6.28318531f*pEvaluationPolynomial[0];
	pEvaluationPolynomial[1]=6.28318531f*pEvaluationPolynomial[1];
	float_complex pAutocorrelation[2];
	computeAutocorrelation2(pAutocorrelation,pEvaluationPolynomial);
	pExponentialMoment[0]=0.5f*pExponentialMoment[0];
	computeCorrelation2(pOutLagrangeMultiplier,pAutocorrelation,pExponentialMoment);
	float normalizationFactor=1.0f/(3.14159265f*pEvaluationPolynomial[0].x);
	pOutLagrangeMultiplier[0]=normalizationFactor*(float_complex(-0.0f,-1.0f)*pOutLagrangeMultiplier[0]);
	pOutLagrangeMultiplier[1]=normalizationFactor*(float_complex(-0.0f,-1.0f)*pOutLagrangeMultiplier[1]);
	pOutLagrangeMultiplier[0]=float_complex(pOutLagrangeMultiplier[0].x,0.0f);
}

void prepareReflectanceSpectrumLagrangeReal3(float pOutLagrangeMultiplier[3],const float pTrigonometricMoment[3]){
	float_complex pExponentialMoment[3];
	trigonometricToExponentialMomentsReal3(pExponentialMoment,pTrigonometricMoment);
	float_complex pEvaluationPolynomial[3];
	levinsonsAlgorithm3(pEvaluationPolynomial,pExponentialMoment);
	pEvaluationPolynomial[0]=6.28318531f*pEvaluationPolynomial[0];
	pEvaluationPolynomial[1]=6.28318531f*pEvaluationPolynomial[1];
	pEvaluationPolynomial[2]=6.28318531f*pEvaluationPolynomial[2];
	float_complex pAutocorrelation[3];
	computeAutocorrelation3(pAutocorrelation,pEvaluationPolynomial);
	pExponentialMoment[0]=0.5f*pExponentialMoment[0];
	computeImaginaryCorrelation3(pOutLagrangeMultiplier,pAutocorrelation,pExponentialMoment);
	float normalizationFactor=1.0f/(3.14159265f*pEvaluationPolynomial[0].x);
	pOutLagrangeMultiplier[0]=normalizationFactor*pOutLagrangeMultiplier[0];
	pOutLagrangeMultiplier[1]=normalizationFactor*pOutLagrangeMultiplier[1];
	pOutLagrangeMultiplier[2]=normalizationFactor*pOutLagrangeMultiplier[2];
}


float evaluateReflectanceSpectrumLagrange2(const float phase,const float_complex pLagrangeMultiplier[2]){
	float_complex conjCirclePoint;
	conjCirclePoint=float_complex(cosf(-phase),sinf(-phase));
	float lagrangeSeries;
	lagrangeSeries=evaluateFourierSeries2(conjCirclePoint,pLagrangeMultiplier);
	return fast_atan(lagrangeSeries)*0.318309886f+0.5f;
}

float evaluateReflectanceSpectrumLagrangeReal3(const float phase,const float pLagrangeMultiplier[3]){
	float_complex conjCirclePoint;
	conjCirclePoint=float_complex(cosf(-phase),sinf(-phase));
	float lagrangeSeries;
	lagrangeSeries=evaluateFourierSeriesReal3(conjCirclePoint,pLagrangeMultiplier);
	return fast_atan(lagrangeSeries)*0.318309886f+0.5f;
}


