//Postgauss appr. spectral weight of 2D Hubbard model in PARA and AF phases. 
//18.5.21

#include "stdafx.h"
#include<iostream>
#include <math.h>
#include <string>
#include <vector>
#include <complex>
#include <array>
#include <memory>
#include <fstream>
#include <string_view>
#include <ctime>
#include <iomanip>
#include <ppl.h>
#include "List.h"
#include "mkl_lapacke.h"
#include "mkl_service.h"

using namespace std;
using namespace concurrency;

//parameters and constants
constexpr int Nx = 4, N = Nx*Nx;
/*Ttab{0.0287141,0.0574282,0.0861423,0.114856,0.117728,0.120599,0.123471,0.126342,0.129213,0.132085,0.134956,0.137828,0.140699}
,0.403333,0.403333,0.403333,0.403333,0.403333,0.403333,0.403333,0.403333}

{0.2552,0.246695,0.216641,0.143289,0.13102,0.116904,0.100163,0.079145,0.048582,3.33333*10^-7,3.33333*10^-7,3.33333*10^-7,3.33333*10^-7}*/
constexpr double T = 1.;
constexpr double tp = -0.;
constexpr double U = 4.0;
constexpr int Nw = 1;
/**/
constexpr double dop = .1;
constexpr double den = 1. - dop;
constexpr double mu = .1;
constexpr double mur = .1;// mu - U*den / 2.;
constexpr double pg = .000000001;//for AF
constexpr double eta = .03;   //disorder or regulator
using Complex = std::complex<double>;
using ComplexVector = List<Complex, N>;
using ComplexVectorAF = List<Complex, N/2>;

double A[Nw];
constexpr Complex I(0., 1.);
constexpr double pi = 3.1415926535897932385;
constexpr double dk = 2.*pi/Nx;

//functions declarations in functions1.h
MKL_Complex16* StdToMkl(Complex* complexPtr);
Complex ixp(Complex y);
int mx(int w); int mxAF(int w); int my(int w);
int cx(int w); int cy(int w);
int vsum(int a, int b);int vsum(int a, int b, int c);
double wb(int n);
Complex ss1a(Complex a, Complex b, Complex c, Complex z, double mur);
Complex ss1b(Complex a, Complex b, Complex c, Complex z, double mur);
Complex ss2a(Complex a, Complex b, Complex c, Complex d, Complex z, double mur);
Complex ss2b(Complex a, Complex b, Complex c, Complex d, Complex z, double mur);
Complex SS1(Complex a, Complex b, Complex c, Complex z, double mur);
Complex SS2(Complex a, Complex b, Complex c, Complex d, Complex z, double mur);
double swP(Complex z, int kx, int ky, double mur);
double swAF(Complex z, int Kx, int Ky, double pg, double mur);

void print(string_view s);

int main()
{
	print("starting...");
	double w[Nw] = { .1};
	Complex z[Nw];
	for (int i = 0; i < Nw; i++) 
	{
        z[i] = pi*T*I;
		cout << "z=" << z[i] << endl;
	}
	cout << "mu=" << mu << "mur=" << mur << ";   U=" << U << ";   T=" << T << endl;
	
	parallel_for(int(0), Nw, [&](int i)
	//for (int i = 0; i < Nw; i++)
	{
		//A[i] = swAF(z[i], 1, 3, mu, pg);
	//}
	});
	cout << "A_AF={" << swAF(z[0], 0, 0, pg, mur) << "};"<< endl;
	cout << "A_P={" << swP(z[0], 0, 0, mur) << "};" << endl;
	print("end");
	getchar();
	return 0;
}
#include "functions1.h"