/**
	@class AllenCurve
	@brief This class provides an interface to the Allen Solar Limb Darkening function.

	@datecreated September 20, 2006

	@version 1.0

	@author
		- Lance Deaver

	@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.

	@see		AllenCurve.h
	@bug		None known

 */


#include "AllenCurve.h"
#include "GATS_Exception.h"
#include "GATS_Utilities.hpp"

#include <iostream>
#include <cmath>

using namespace std;
using namespace GATS_Utilities;

AllenCurve::AllenCurve()
{
}

AllenCurve::~AllenCurve()
{
}
/**
 	Routine returns the solar limb darkening curve calculated according
 	the the "Allen" formulation.  The limb darkening curve is a fractional
  	solar intensity (1 at sun center 0 off the Sun edge) versus angular
  	displacement from Sun top Edge as viewed from the Earth.

 	@author
 		- Mark Hervig
 	@param[in] angle	The angle from the top edge of the Solar disk (arcminutes). If this value is
 			< 0 this will indicate the angle is above the top edge.  A value that is > 0
 			indicates an angle that is lower or Earthward from the top edge so from Earth
 			a value of 16 arcminutes would roughly be solar center.
 	@param[in]  wavelength The requested wavelength for the Solar intensity (microns).
 			This value must be within  0.26 and 10.5 microns.
	@param[in]	distance The distance to the Sun (km).

 	@retval  double The normalized solar intensity (compared to Sun center).

 	@exception	ValueOutOfRange  Input distance or wavelength are not within required limits.
 	
 	@note	This routine has been modified from the Original model used from HALORET.  This modification
 			is an additional U and V coefficient at .26 microns.  The original code only went to
 			.31 microns which did not model the SOFIE Band 1 channel.
*/
double AllenCurve::getSolarIntensity(const double angle, const double wavelength, const double distance)
{
	//- setup some variables

	const double radsun	= 6.9595e5;	//km
	const double pi		= M_PI;
	const double convert	= pi/10800;

	static double outcon	= 0.0;
	static double wavnow	= -1.0 ;   //always force a first time calculation
	static double urun		= 0.95;
	static double vrun		= -0.20;
	static double constant	= 0.05;
	static double acon		= 0.25;
//	static bool   MODi		= true;


	const 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};
	const 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};
	const 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};

	const int iallen = 15;


	if(distance < 10.f) {
		std::string errMsg ="AllenCurve::getSolarIntensity: Input distance:" + ConvertToString(distance) + " must be greater 10 km.";
		THROW_GATS_EXCEPTION(ValueOutOfRange, errMsg);
	}
	//- do the calculation

//	bool gotWave	= false;
	double intensity = 0.0;
	double sigma  = asin(1000.0/distance)/3;
	double sigsqr = 0.5/(sigma*sigma);
	double ratio  = distance/radsun;
	double edge   = asin(1.0/ratio);

	if (wavnow != wavelength)
	{
		if ((wavelength > wave[iallen-1]) || (wavelength < wave[0]))
		{
			std::string errMsg ="AllenCurve::getSolarIntensity: Input wavelength:" + ConvertToString(wavelength)	 + " must be within "
				+ ConvertToString(wave[0]) + " and " + ConvertToString(wave[iallen-1]) + " microns.";
			THROW_GATS_EXCEPTION(ValueOutOfRange, errMsg );
		}

		for (int i = 1; i < iallen; i++)
		{
			if (wave[i] >= wavelength)
			{
				urun		= (U[i] * (wavelength - wave[i-1]) + U[i-1] * (wave[i] - wavelength)) / (wave[i] - wave[i-1]);
				vrun		= (V[i] * (wavelength - wave[i-1]) + V[i-1] * (wave[i] - wavelength)) / (wave[i] - wave[i-1]);
				constant	= 1 - urun;
				acon		= constant - vrun;
				wavnow		= wavelength;
				outcon = pi * urun / (edge * acon);
//				MODi		= true;
//				gotWave		= true;
				break;
			}
		}

//		if (!gotWave)
//		{
//			cerr << "wavelength out of bounds in getSolarIntensity: " << wavelength << endl;
//			exit(1);
//		}
	}

//	if (MODi)
//	{
//		MODi = false;
//		outcon = pi * urun / (edge * acon);
//	}

	double theta = abs(angle) * convert;
//addition added to reference angles relative to top edge;  ET/LD

	theta = (angle <0.f) ?  theta+edge : theta-edge;
// end of additions;
	if (theta <= edge)
	{
		double sinthsq = pow((ratio * sin(theta)),2);
		intensity  = constant + urun * sqrt(1.0 - sinthsq) - vrun * sinthsq;
	}
	else
	{
		double delta   = theta - edge;
		double coefficient   = delta * (delta * sigsqr + outcon);
		 if (coefficient <= 50.0)
		 {
		 	intensity = acon * exp(-coefficient);
		 }
	}

	//- done

   return (intensity == intensity) ? intensity: 0.0 ;
 }