#include "pbox.h"

Block NewLeft(const Block& b1, const Block& b2, const Vector& bas)
    {
    Block res;
    res.H11 = bas(1) * bas(1) * b1.H11 + bas(2) * bas(2) * 2.0 
    			- 2 * bas(1) * bas(2) * b1.L_inner;
    res.L_inner = bas(2);

    return res;
    }

Block NewRight(const Block& b3, const Block& b4, const Vector& bas)
    {
    Block res;
    res.H11 = bas(1) * bas(1) * b4.H11 + bas(2) * bas(2) * 2.0 
    			- 2 * bas(1) * bas(2) * b4.L_inner;
    res.L_inner = bas(2);

    return res;
    }

Real System::GetGroundState(WaveFunction& p)
    {
    Matrix H(4,4), evecs(4,4);
    Vector evals(4);
    H = 0.0;
    H(1,1) = b1.H11;
    H(2,2) = b2.H11;
    H(3,3) = b3.H11;
    H(4,4) = b4.H11;
    H(1,2) = H(2,1) = -b1.L_inner * b2.L_inner;
    H(2,3) = H(3,2) = -b2.L_inner * b3.L_inner;
    H(3,4) = H(4,3) = -b3.L_inner * b4.L_inner;
    EigenValues(H,evals,evecs);
    p.v = evecs.Column(1);
    if(p.v.sumels() < 0.0) p.v *= -1.0;
    Real energy = evals(1);
    return energy;
    }