#define THIS_IS_MAIN

#include "pbox.h"

int exlen = 0;
double exacten(int length);
double exact(int length, int i);

int main()
    {
    Block siteblock;
    cerr << "Input length, number of iterations: ";
    int i, length, nsweeps;
    cin >> length >> nsweeps;
    cout << "length, nsweeps = " << length SP nsweeps << endl;
    cout << "Exact energy = " << exacten(length) << endl;
    exlen = length;
    Array1<Block> allblocks(length);
    WaveFunction psi;
    Real energy;

// Warmup sweep
    allblocks[1] = siteblock;
    for(i = 1; i < length/2; i++)
	{
	Block rightblock = allblocks(i).Reflect();

	System S(allblocks(i),siteblock,siteblock,rightblock);

	energy = S.GetGroundState(psi);
	cout << i+1 SP psi.v(2) SP energy SP 0 << endl;

	DensityMatrix rho(psi,Left);
	Vector basis = rho.NewBasis();

	allblocks[i+1] = NewLeft(allblocks(i),siteblock,basis);
	}

// Finite System sweeps
    for(int swp = 1; swp <= nsweeps; swp++)
	{
// We assume reflection symmetry:
	allblocks[length/2 + 2] = allblocks(length/2 - 1).Reflect();
	cout << endl;

	for(i = length/2+2; i > 3; i--)
	    {
	    System S(allblocks(i-3),siteblock,siteblock,allblocks(i));

	    energy = S.GetGroundState(psi);
	    cout << i-1 SP psi.v(3) SP energy SP swp - 0.5 << endl;

	    DensityMatrix rho(psi,Right);
	    Vector basis = rho.NewBasis();

	    allblocks[i-1] = NewRight(siteblock,allblocks(i),basis);
	    }

	cout << endl << 1 SP psi.v(1) SP energy SP swp << endl;
	for(i = 1; i < length/2-1; i++)
	    {
	    System S(allblocks(i),siteblock,siteblock,allblocks(i+3));

	    energy = S.GetGroundState(psi);
	    cout << i+1 SP psi.v(2) SP energy SP swp << endl;

	    DensityMatrix rho(psi,Left);
	    Vector basis = rho.NewBasis();

	    allblocks[i+1] = NewLeft(allblocks(i),siteblock,basis);
	    }
	}
    return 0;
    }

static const Real Pi = M_PI;

double exacten(int length)
    {
    double q = Pi / (length + 1.0);
    return 2 * (1.0 - cos(q));
    }

double exact(int length, int i)
    {
    double q = Pi / (length + 1.0);
    double norm = 0.0;
    for(int j=1; j <= length; j++)
	norm += sin(j * q) * sin(j * q);
    norm = 1.0 / sqrt(norm);
    return sin(i * q) * norm;
    }