/*************************************************************************
 *  Copyright (C) 2015 by Andrius Štikonas <andrius@stikonas.eu>         *
 *                                                                       *
 *  This program is free software; you can redistribute it and/or modify *
 *  it under the terms of the GNU General Public License as published by *
 *  the Free Software Foundation; either version 3 of the License, or    *
 *  (at your option) any later version.                                  *
 *                                                                       *
 *  This program 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/>.*
 *************************************************************************/

// Calculation of mutual information in the setup of http://arxiv.org/abs/1503.08161

#include <cmath>
#include <complex>
#include <iostream>

static double alpha, beta;
std::complex<double> crossRatio (std::complex<double>, std::complex<double>, std::complex<double>, std::complex<double>);
double Fitzpatrick(std::complex<double>, std::complex<double>, double);

int main()
{
	// Parameters that can be changed:
	double tPlus = 0;  // time on the left boundary
	double tMinus = 0; // time on the right boundary
	alpha = 0.4; // encodes the conformal dimention of local operator h_Psi
	double y = 3; // endpoint of interval A
	double L = 15; // length of interval A
	double epsilon = 0.01; // smearing parameter
	beta = 10; // inverse temperature
// 	double tOmega = /*5.4418*/; // thermal state is perturbed by operator inserted at time -tOmega
	double c = 600; // central charge. Must be large in our approximation
	// End of parameters
	for(unsigned int i = 0; i < 1000; ++i)
	{
		double tOmega = i*2*L/1000;

		// Operator insertion points: Left boundary
		std::complex<double> x1 (0, -epsilon), x4 (0, epsilon), x1bar, x4bar;
		x1bar = conj(x1);
		x4bar = conj(x4);
		double x2    = y - tOmega - tMinus;
		double x2bar = y + tOmega + tMinus;
		double x3    = L + x2;
		double x3bar = L + x2bar;

		// Operator insertion points: Right boundary
		std::complex<double> x6    (y - tPlus - tOmega, beta/2);
		std::complex<double> x6bar (y + tPlus + tOmega, -beta/2);
		std::complex<double> x5    = L + x6;
		std::complex<double> x5bar = L + x6bar;

		// Cross-ratios for S_A
		std::complex<double> zA    = crossRatio(x1, x2, x3, x4);
		std::complex<double> zAbar = crossRatio(x1bar, x2bar, x3bar, x4bar);

		// Cross-ratios for S_B
		std::complex<double> zB    = crossRatio(x1, x5, x6, x4);
		std::complex<double> zBbar = crossRatio(x1bar, x5bar, x6bar, x4bar);

		// Cross-ratios for S_{A union B}
		std::complex<double> z2    = crossRatio(x1, x2, x6, x4);
		std::complex<double> z2bar = crossRatio(x1bar, x2bar, x6bar, x4bar);
		std::complex<double> z5    = crossRatio(x1, x5, x3, x4);
		std::complex<double> z5bar = crossRatio(x1bar, x5bar, x3bar, x4bar);

		// Now we calculate entanglement entropies using Fitzpatrick, Kaplan, Walters formula.
		double S_A = c/6 * log(Fitzpatrick(zA, zAbar, y < tMinus + tOmega && tMinus + tOmega < y + L ? 2*M_PI : 0));
		double S_B = c/6 * log(Fitzpatrick(zB, zBbar, 0));
		double S_union = c/6 * log(Fitzpatrick(z2, z2bar, y < tMinus + tOmega ? 2*M_PI : 0) * Fitzpatrick(z5, z5bar, y +L < tMinus + tOmega ? -2*M_PI : 0));
		double S_union_analytic = c/3 * log(beta/M_PI/epsilon * sin(M_PI*alpha)/alpha) + c/6* log(sinh(M_PI/beta*(tMinus + tOmega - y)) * cosh(M_PI/beta*(tPlus + tOmega - y)) * sinh(M_PI/beta*(tMinus + tOmega - y - L)) * cosh(M_PI/beta*(tPlus + tOmega - y - L)));
		double S_thermal = 2*c/3 * log(sinh(M_PI*L/beta)/cosh(M_PI/beta*(tMinus-tPlus)));
		double I = S_A + S_B - S_union + S_thermal;
		std::cout << tOmega << "\t" << I << std::endl;
}
	return 0;
}

std::complex<double> crossRatio (std::complex<double> x1, std::complex<double> x2, std::complex<double> x3, std::complex<double> x4)
{
	double pb = M_PI/beta;
	return sinh(pb*(x1-x2))*sinh(pb*(x3-x4))/sinh(pb*(x1-x3))/sinh(pb*(x2-x4));
}

double Fitzpatrick(std::complex<double> z, std::complex<double> zbar, double phase)
{
	// phase is necessary to take into account nontrivial monodromy
	std::complex<double> one = 1; // one as a complex number
	std::complex<double> i(0,1); // imaginary unit
	double alphaExpr = 0.5-alpha/2;
	std::complex<double> exponent1 = exp(phase*i*alphaExpr);
	std::complex<double> exponent2 = exp(phase*i*alpha);
	return real(exponent1 * pow(z, alphaExpr) * pow(zbar, alphaExpr)
			* (one - exponent2 * pow(z, alpha)) * (one - pow(zbar, alpha))
			/ ( alpha*alpha * (one-z) * (one-zbar) ) );
}