// pbox.h

#include <iostream.h>

#ifndef _pbox_h
#define _pbox_h
#include "matrix.h"
#include <math.h>

class Block
    {
public:
    Real H11;
    Real L_inner;

    Block()			// Default: construct a one-site block
	: H11(2.0), L_inner(1.0) { }
    Block Reflect() const { return *this; }

    ARRAY1H_DEFS(Block)
    };

inline ostream& operator << (ostream &s, const Block &l) 
    { return s << l.H11 << " " << l.L_inner; }

Block NewLeft(const Block& b1, const Block& b2,const Vector&);
Block NewRight(const Block& b3, const Block& b4,const Vector&);

class WaveFunction;

class System
    {
public:
    const Block& b1;
    const Block& b2;
    const Block& b3;
    const Block& b4;

    System(const Block& bb1,const Block&bb2, 
    		const Block&bb3, const Block&bb4)
	: b1(bb1), b2(bb2), b3(bb3), b4(bb4) { }

    Real GetGroundState(WaveFunction& p);
    };

class WaveFunction
    {
public:
    Vector v;
    WaveFunction() : v(4) {}
    };

enum LR {Left, Right};

class DensityMatrix
    {
public:
    Real a,b;

    DensityMatrix(const WaveFunction& psi,LR lr) 
	{
	if(lr == Left)
	    { a = psi.v(1); b = psi.v(2); }
	else
	    { a = psi.v(4); b = psi.v(3); }
	}
    Vector NewBasis()
	{
	Vector res(2);
	if(a == 0.0 && b == 0.0)
	    a = 1.0;
	Real norm = sqrt(a*a+b*b);
	res(1) = a / norm;
	res(2) = b / norm;
	return res;
	}
    };

#ifdef THIS_IS_MAIN

ARRAY1CC_DEFS(Block)

#endif

#endif