/*
    This file is part of Mitsuba, a physically based rendering system.

    Copyright (c) 2007-2014 by Wenzel Jakob and others.

    Mitsuba is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License Version 3
    as published by the Free Software Foundation.

    Mitsuba is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#if !defined(__SUN_H)
#define __SUN_H

#include <mitsuba/mitsuba.h>

#define EARTH_MEAN_RADIUS 6371.01	// In km
#define ASTRONOMICAL_UNIT 149597890	// In km

MTS_NAMESPACE_BEGIN

struct DateTimeRecord {
	int year;
	int month;
	int day;
	Float hour;
	Float minute;
	Float second;

	std::string toString() const {
		std::ostringstream oss;
		oss << "DateTimeRecord[year = " << year
			<< ", month= " << month
			<< ", day = " << day
			<< ", hour = " << hour
			<< ", minute = " << minute
			<< ", second = " << second << "]";
		return oss.str();
	}
};

struct LocationRecord {
	Float longitude;
	Float latitude;
	Float timezone;

	std::string toString() const {
		std::ostringstream oss;
		oss << "LocationRecord[latitude = " << latitude
			<< ", longitude = " << longitude
			<< ", timezone = " << timezone << "]";
		return oss.str();
	}
};

struct SphericalCoordinates {
	Float elevation;
	Float azimuth;

	inline SphericalCoordinates() { }

	inline SphericalCoordinates(Float elevation, Float azimuth)
		: elevation(elevation), azimuth(azimuth) { }

	inline SphericalCoordinates(Stream *stream) {
		elevation = stream->readFloat();
		azimuth = stream->readFloat();
	}

	void serialize(Stream *stream) const {
		stream->writeFloat(elevation);
		stream->writeFloat(azimuth);
	}

	std::string toString() const {
		std::ostringstream oss;
		oss << "SphericalCoordinates[elevation = " << radToDeg(elevation)
			<< ", azimuth = " << radToDeg(azimuth) << "]";
		return oss.str();
	}
};

Vector toSphere(const SphericalCoordinates coords) {
	Float sinTheta, cosTheta, sinPhi, cosPhi;

	math::sincos(coords.elevation, &sinTheta, &cosTheta);
	math::sincos(coords.azimuth, &sinPhi, &cosPhi);

	return Vector(sinPhi*sinTheta, cosTheta, -cosPhi*sinTheta);
}

SphericalCoordinates fromSphere(const Vector &d) {
	Float azimuth = std::atan2(d.x, -d.z);
	Float elevation = math::safe_acos(d.y);
	if (azimuth < 0)
		azimuth += 2*M_PI;
	return SphericalCoordinates(elevation, azimuth);
}

/**
 * \brief Compute the elevation and azimuth of the sun as seen by an observer
 * at \c location at the date and time specified in \c dateTime.
 *
 * Based on "Computing the Solar Vector" by Manuel Blanco-Muriel,
 * Diego C. Alarcon-Padilla, Teodoro Lopez-Moratalla, and Martin Lara-Coira,
 * in "Solar energy", vol 27, number 5, 2001 by Pergamon Press.
 */
SphericalCoordinates computeSunCoordinates(const DateTimeRecord &dateTime, const LocationRecord &location) {
	// Main variables
	double elapsedJulianDays, decHours;
	double eclipticLongitude, eclipticObliquity;
	double rightAscension, declination;
	double elevation, azimuth;

	// Auxiliary variables
	double dY;
	double dX;

	/* Calculate difference in days between the current Julian Day
	   and JD 2451545.0, which is noon 1 January 2000 Universal Time */
	{
		// Calculate time of the day in UT decimal hours
		decHours = dateTime.hour - location.timezone +
			(dateTime.minute + dateTime.second / 60.0 ) / 60.0;

		// Calculate current Julian Day
		int liAux1 = (dateTime.month-14) / 12;
		int liAux2 = (1461*(dateTime.year + 4800 + liAux1)) / 4
			+ (367 * (dateTime.month - 2 - 12 * liAux1)) / 12
			- (3 * ((dateTime.year + 4900 + liAux1) / 100)) / 4
			+ dateTime.day - 32075;
		double dJulianDate = (double) liAux2 - 0.5 + decHours / 24.0;

		// Calculate difference between current Julian Day and JD 2451545.0
		elapsedJulianDays = dJulianDate - 2451545.0;
	}

	/* Calculate ecliptic coordinates (ecliptic longitude and obliquity of the
	   ecliptic in radians but without limiting the angle to be less than 2*Pi
	   (i.e., the result may be greater than 2*Pi) */
	{
		double omega = 2.1429 - 0.0010394594 * elapsedJulianDays;
		double meanLongitude = 4.8950630 + 0.017202791698 * elapsedJulianDays; // Radians
		double anomaly = 6.2400600 + 0.0172019699 * elapsedJulianDays;

		eclipticLongitude = meanLongitude + 0.03341607 * std::sin(anomaly)
			+ 0.00034894 * std::sin(2*anomaly) - 0.0001134
			- 0.0000203 * std::sin(omega);

		eclipticObliquity = 0.4090928 - 6.2140e-9 * elapsedJulianDays
			+ 0.0000396 * std::cos(omega);
	}

	/* Calculate celestial coordinates ( right ascension and declination ) in radians
	   but without limiting the angle to be less than 2*Pi (i.e., the result may be
	   greater than 2*Pi) */
	{
		double sinEclipticLongitude = std::sin(eclipticLongitude);
		dY = std::cos(eclipticObliquity) * sinEclipticLongitude;
		dX = std::cos(eclipticLongitude);
		rightAscension = std::atan2(dY, dX);
		if (rightAscension < 0.0)
			rightAscension += 2*M_PI;
		declination = std::asin(std::sin(eclipticObliquity) * sinEclipticLongitude);
	}

	// Calculate local coordinates (azimuth and zenith angle) in degrees
	{
		double greenwichMeanSiderealTime = 6.6974243242
			+ 0.0657098283 * elapsedJulianDays + decHours;

		double localMeanSiderealTime = degToRad((Float) ((greenwichMeanSiderealTime * 15
			+ location.longitude)));

		double latitudeInRadians = degToRad(location.latitude);
		double cosLatitude = std::cos(latitudeInRadians);
		double sinLatitude = std::sin(latitudeInRadians);

		double hourAngle = localMeanSiderealTime - rightAscension;
		double cosHourAngle = std::cos(hourAngle);

		elevation = std::acos(cosLatitude * cosHourAngle
			* std::cos(declination) + std::sin(declination) * sinLatitude);

		dY = -std::sin(hourAngle);
		dX = std::tan(declination) * cosLatitude - sinLatitude * cosHourAngle;

		azimuth = std::atan2(dY, dX);
		if (azimuth < 0.0)
			azimuth += 2*M_PI;

		// Parallax Correction
		elevation += (EARTH_MEAN_RADIUS / ASTRONOMICAL_UNIT) * std::sin(elevation);
	}

	return SphericalCoordinates((Float) elevation, (Float) azimuth);
}

SphericalCoordinates computeSunCoordinates(const Vector& sunDir, const Transform &worldToLuminaire) {
	return fromSphere(normalize(worldToLuminaire(sunDir)));
}

SphericalCoordinates computeSunCoordinates(const Properties &props) {
	/* configure position of sun */
	if (props.hasProperty("sunDirection")) {
		if (props.hasProperty("latitude") || props.hasProperty("longitude")
			|| props.hasProperty("timezone") || props.hasProperty("day")
			|| props.hasProperty("time"))
			SLog(EError, "Both the 'sunDirection' parameter and time/location "
					"information were provided -- only one of them can be specified at a time!");

		return computeSunCoordinates(
			props.getVector("sunDirection"),
			props.getTransform("toWorld", Transform()).inverse());
	} else {
		LocationRecord location;
		DateTimeRecord dateTime;

		location.latitude  = props.getFloat("latitude", 35.6894f);
		location.longitude = props.getFloat("longitude", 139.6917f);
		location.timezone  = props.getFloat("timezone", 9);
		dateTime.year      = props.getInteger("year", 2010);
		dateTime.day       = props.getInteger("day", 10);
		dateTime.month     = props.getInteger("month", 7);
		dateTime.hour      = props.getFloat("hour", 15.0f);
		dateTime.minute    = props.getFloat("minute", 0.0f);
		dateTime.second    = props.getFloat("second", 0.0f);

		SphericalCoordinates coords = computeSunCoordinates(dateTime, location);

		SLog(EDebug, "Computed sun position for %s and %s: %s",
			location.toString().c_str(), dateTime.toString().c_str(),
			coords.toString().c_str());

		return coords;
	}
}

/* The following is from the implementation of "A Practical Analytic Model for
   Daylight" by A.J. Preetham, Peter Shirley, and Brian Smits */

/* All data lifted from MI. Units are either [] or cm^-1. refer when in doubt MI */

// k_o Spectrum table from pg 127, MI.
Float k_oWavelengths[64] = {
	300, 305, 310, 315, 320, 325, 330, 335, 340, 345,
	350, 355, 445, 450, 455, 460, 465, 470, 475, 480,
	485, 490, 495, 500, 505, 510, 515, 520, 525, 530,
	535, 540, 545, 550, 555, 560, 565, 570, 575, 580,
	585, 590, 595, 600, 605, 610, 620, 630, 640, 650,
	660, 670, 680, 690, 700, 710, 720, 730, 740, 750,
	760, 770, 780, 790
};

Float k_oAmplitudes[65] = {
	10.0, 4.8, 2.7, 1.35, .8, .380, .160, .075, .04, .019, .007,
	.0, .003, .003, .004, .006, .008, .009, .012, .014, .017,
	.021, .025, .03, .035, .04, .045, .048, .057, .063, .07,
	.075, .08, .085, .095, .103, .110, .12, .122, .12, .118,
	.115, .12, .125, .130, .12, .105, .09, .079, .067, .057,
	.048, .036, .028, .023, .018, .014, .011, .010, .009,
	.007, .004, .0, .0
};

// k_g Spectrum table from pg 130, MI.
Float k_gWavelengths[4] = {
	759, 760, 770, 771
};

Float k_gAmplitudes[4] = {
	0, 3.0, 0.210, 0
};

// k_wa Spectrum table from pg 130, MI.
Float k_waWavelengths[13] = {
	689, 690, 700, 710, 720,
	730, 740, 750, 760, 770,
	780, 790, 800
};

Float k_waAmplitudes[13] = {
	0, 0.160e-1, 0.240e-1, 0.125e-1,
	0.100e+1, 0.870, 0.610e-1, 0.100e-2,
	0.100e-4, 0.100e-4, 0.600e-3,
	0.175e-1, 0.360e-1
};

/* Wavelengths corresponding to the table below */
Float solWavelengths[38] = {
	380, 390, 400, 410, 420, 430, 440, 450,
	460, 470, 480, 490, 500, 510, 520, 530,
	540, 550, 560, 570, 580, 590, 600, 610,
	620, 630, 640, 650, 660, 670, 680, 690,
	700, 710, 720, 730, 740, 750
};

/* Solar amplitude in watts / (m^2 * nm * sr) */
Float solAmplitudes[38] = {
	16559.0, 16233.7, 21127.5, 25888.2, 25829.1,
	24232.3, 26760.5, 29658.3, 30545.4, 30057.5,
	30663.7, 28830.4, 28712.1, 27825.0, 27100.6,
	27233.6, 26361.3, 25503.8, 25060.2, 25311.6,
	25355.9, 25134.2, 24631.5, 24173.2, 23685.3,
	23212.1, 22827.7, 22339.8, 21970.2, 21526.7,
	21097.9, 20728.3, 20240.4, 19870.8, 19427.2,
	19072.4, 18628.9, 18259.2
};

Spectrum computeSunRadiance(Float theta, Float turbidity) {
    InterpolatedSpectrum k_oCurve(k_oWavelengths, k_oAmplitudes, 64);
   	InterpolatedSpectrum k_gCurve(k_gWavelengths, k_gAmplitudes, 4);
    InterpolatedSpectrum k_waCurve(k_waWavelengths, k_waAmplitudes, 13);
    InterpolatedSpectrum solCurve(solWavelengths, solAmplitudes, 38);
    Float data[91], wavelengths[91];  // (800 - 350) / 5  + 1

    Float beta = 0.04608365822050f * turbidity - 0.04586025928522f;

	// Relative Optical Mass
    Float m = 1.0f / (std::cos(theta) + 0.15f *
		std::pow(93.885f - theta/M_PI*180.0f, (Float) -1.253f));

    Float lambda;
	int i = 0;
    for(i = 0, lambda = 350; i < 91; i++, lambda += 5) {
		// Rayleigh Scattering
		// Results agree with the graph (pg 115, MI) */
		Float tauR = math::fastexp(-m * 0.008735f * std::pow(lambda/1000.0f, (Float) -4.08));

		// Aerosol (water + dust) attenuation
		// beta - amount of aerosols present
		// alpha - ratio of small to large particle sizes. (0:4,usually 1.3)
		// Results agree with the graph (pg 121, MI)
		const Float alpha = 1.3f;
		Float tauA = math::fastexp(-m * beta * std::pow(lambda/1000.0f, -alpha));  // lambda should be in um

		// Attenuation due to ozone absorption
		// lOzone - amount of ozone in cm(NTP)
		// Results agree with the graph (pg 128, MI)
		const Float lOzone = .35f;
		Float tauO = math::fastexp(-m * k_oCurve.eval(lambda) * lOzone);

		// Attenuation due to mixed gases absorption
		// Results agree with the graph (pg 131, MI)
		Float tauG = math::fastexp(-1.41f * k_gCurve.eval(lambda) * m / std::pow(1 + 118.93f
			* k_gCurve.eval(lambda) * m, (Float) 0.45f));

		// Attenuation due to water vapor absorbtion
		// w - precipitable water vapor in centimeters (standard = 2)
		// Results agree with the graph (pg 132, MI)
		const Float w = 2.0;
		Float tauWA = math::fastexp(-0.2385f * k_waCurve.eval(lambda) * w * m /
				std::pow(1 + 20.07f * k_waCurve.eval(lambda) * w * m, (Float) 0.45f));

		data[i] = solCurve.eval(lambda) * tauR * tauA * tauO * tauG * tauWA;
		wavelengths[i] = lambda;
    }

    InterpolatedSpectrum interpolated(wavelengths, data, 91);
	Spectrum discretized;
	discretized.fromContinuousSpectrum(interpolated);
	discretized.clampNegative();

	return discretized;
}


MTS_NAMESPACE_END

#endif /* __SUN_H */
