package dmrgvisu;

import info.monitorenter.gui.chart.traces.Trace2DReplacing;

import org.netlib.lapack.DSYEV;
import org.netlib.util.*;

import javax.swing.JTextArea;

public class DMRG
{
	public int length;
	public int size;
	public int sweeps;
	public WaveFunction psi;
	public block siteblock;
	public block [] allblocks;
	public double energy;
	public double wf;
	public int i;
	public int sweep;
	public int EnergyRead;
	public int WFRead;
	public double y;
	public int lock;
	public int potentialChooser;
	public double h;
	public double OldY;
	public double OldWf;
	public double OldEn;
	public int FullHalf;
	public double MinEnergy;
	/** Energieausgabe */
	public JTextArea EnergyText;

	/** Sweepausgabe */
	public JTextArea SweepText;
	
	public Trace2DReplacing exWFTrace;
	public Trace2DReplacing exEnTrace;

    public DMRG() 
	{
		this.FullHalf = 0;
		this.MinEnergy = 10000;
        this.length = 10000;			// Anzahl der Stützstellen; Größe des Arrays
		this.size = 10;					// Halbe Größe des Systems (System=[-size:size])
		this.sweeps = 100000;			// Maximalanzahl der Sweeps
		this.psi = new WaveFunction();	// Wellenfunktion des Systems
		double energy = 0.0;			// Energie am Ort i
		this.i = 1;						// Derzeitiger Rechenschritt
		this.sweep = 0;					// Derzeitiger Sweep
		this.EnergyRead = -1;			// -1: 1. Durchlauf; 0: Energie wurde noch nicht gelesen; 1: Energie wurde gelesen.
		this.WFRead = -1;				// -1: 1. Durchlauf; 0: Wellenfunktion wurde noch nicht gelesen; 1: Wellenfunktion wurde gelesen.
		this.y = 0;						// Derzeitiger Ausgabeort (zumeist ungleich i)
		this.wf = 0;					// Normierter Anteil der Wellenfunktion am Ort i
		this.lock = 0;					// 1: Berechnung wird ausgeführt; 0: Berechnung wird nicht ausgeführt
		this.potentialChooser = 0;		// Auswahl des Potentials (0: Kein Potential; 1: x^2-Potential; 2: x^2+10*x^4-Potential; 3: -x^2+0.02*x^4-Potential; 4: Kastenpotential)
		this.h = 2.0*this.size/this.length;		// Diskretisierungsschrittgröße
		this.siteblock = new block(this);		// Erzeugung des Standard-Siteblocks
		this.allblocks = new block[length];		// Erzeugung des Block-Arrays
		for (int i=0; i<this.length; i++)		// Befüllung des Block-Arrays mit Standard-Blöcken
			this.allblocks[i] = new block(this);
		allblocks[1] = siteblock;		// Befüllung des ersten (nicht nullten) Blocks mit einem Block mit Potential
		this.OldY = 0;
		this.OldWf = 0.0;
		this.OldEn = 0.0;
    }
	
    public DMRG(int Samplings, int SystemSize, int Sweeps, int Potential) 
	{
		this.FullHalf = 0;
		this.MinEnergy = Samplings;
		//System.out.println(Potential);
	        this.length = Samplings;			// Anzahl der Stützstellen; Größe des Arrays
		this.size = SystemSize;				// Halbe Größe des Systems (System=[-size:size])
		this.sweeps = Sweeps;			// Maximalanzahl der Sweeps
		this.psi = new WaveFunction();	// Wellenfunktion des Systems
		double energy = 0.0;			// Energie am Ort i
		this.i = 0;						// Derzeitiger Rechenschritt
		this.sweep = 0;					// Derzeitiger Sweep
		this.EnergyRead = -1;			// -1: 1. Durchlauf; 0: Energie wurde noch nicht gelesen; 1: Energie wurde gelesen.
		this.WFRead = -1;				// -1: 1. Durchlauf; 0: Wellenfunktion wurde noch nicht gelesen; 1: Wellenfunktion wurde gelesen.
		this.y = 0;						// Derzeitiger Ausgabeort (zumeist ungleich i)
		this.wf = 0;					// Normierter Anteil der Wellenfunktion am Ort i
		this.lock = 0;					// 1: Berechnung wird ausgeführt; 0: Berechnung wird nicht ausgeführt
		this.potentialChooser = Potential;	// Auswahl des Potentials (0: Kein Potential; 1: x^2-Potential; 2: x^2+10*x^4-Potential; 3: -x^2+0.02*x^4-Potential; 4: Kastenpotential)
		this.h = 2.0*this.size/this.length;		// Diskretisierungsschrittgröße
		this.siteblock = new block(this);		// Erzeugung des Standard-Siteblocks
		this.allblocks = new block[length];		// Erzeugung des Block-Arrays
		for (int j=0; j<this.length; j++)		// Befüllung des Block-Arrays mit Standard-Blöcken
		{
			this.allblocks[j] = new block(this);
			this.allblocks[j].H11 += this.V(j);	// Macht keinen Unterschied...
		}
//		allblocks[i] = siteblock;		// Befüllung des nullten Blocks mit einem Block mit Potential
//		allblocks[0].H11 += this.V(0);
		this.OldY = 0;
		this.OldWf = 0.0;
		this.OldEn = 0.0;
    }
	
	public static block NewLeft(block b1, block b2, double []bas, double P, double h)
	{
		block res = new block();
		res.H11 = bas[0]*bas[0]*b1.H11 + bas[1]*bas[1]*(b2.H11+P) - 2.0*bas[0]*bas[1]*b1.L_inner;
		res.L_inner = bas[1]/Math.pow(h, 2.0);
		return res;
	}
	
	public static block NewRight(block b3, block b4, double []bas, double P, double h)
	{
		block res = new block();
		res.H11 = bas[0]*bas[0]*b4.H11 + bas[1]*bas[1]*(b3.H11+P) - 2.0*bas[0]*bas[1]*b4.L_inner;
		res.L_inner = bas[1]/Math.pow(h, 2.0);
		return res;
	}

	public void DoNextCalc()
	{
		if ((this.sweep-1)/2 == this.sweeps)
			return;
		if (this.sweep == 0)
		{
			if (this.i < this.length-4)
			{
				block rightblock = this.allblocks[this.i];
				DMRGsystem S = new DMRGsystem(this.allblocks[this.i], this.siteblock, this.siteblock, rightblock);
				this.energy = S.GetGroundState(this.psi, this);
				this.y = 2.0*(((double)this.i+1.0)/(this.length-1.0)-0.5)*this.size;
				this.wf = this.psi.V[1];
				DensityMatrix rho = new DensityMatrix(this.psi, -1);
				double [] basis = rho.NewBasis();
				this.allblocks[this.i+1] = DMRG.NewLeft(this.allblocks[this.i], siteblock, basis, this.V(i+1), this.h);
				this.i++;
			}
			else
			{
				this.sweep = 1;
				this.SweepText.setText("Current sweep: " + ((this.sweep+1)/2));
				this.i = this.length-1;
			}
		}
		else
		{
			if (this.sweep%2 == 1)
			{
				if (this.i > 3)
				{
					DMRGsystem S = new DMRGsystem(this.allblocks[this.i-3], this.siteblock, this.siteblock, this.allblocks[this.i]);
					this.energy = S.GetGroundState(this.psi, this);
					this.y = 2.0*(((double)this.i-1.0)/(this.length-1.0)-0.5)*this.size;
					this.wf = this.psi.V[2];
					DensityMatrix rho = new DensityMatrix(this.psi, 1);
					double [] basis = rho.NewBasis();
					this.allblocks[this.i-1] = DMRG.NewRight(this.siteblock, this.allblocks[this.i], basis, this.V(i-1), this.h);
					this.i--;
				}
				else
				{
					if (this.i == -1)
					{
					//	this.y = 2.0*(((double)1)/(this.length-1.0)-0.5)*this.size;
					//	this.wf = this.psi.V[0];
						this.sweep++;
						this.SweepText.setText("Current sweep: " + ((this.sweep+1)/2));
						this.i = 0;
					}
					else
					{
						this.i = -1;
					}
				}
			}
			else
			{
				if (this.i < this.length-4)
				{
					DMRGsystem S = new DMRGsystem(this.allblocks[this.i], this.siteblock, this.siteblock, this.allblocks[this.i+3]);
					this.energy = S.GetGroundState(this.psi, this);
					this.y = 2.0*(((double)this.i+1.0)/(this.length-1.0)-0.5)*this.size;
					this.wf = this.psi.V[1];
					DensityMatrix rho = new DensityMatrix(this.psi, -1);
					double [] basis = rho.NewBasis();
					this.allblocks[this.i+1] = DMRG.NewLeft(this.allblocks[this.i], this.siteblock, basis, this.V(i+1), this.h);
					this.i++;
				}
				else
				{
					this.i = (this.FullHalf==0?this.length-1:this.length/2+2);
					this.sweep++;
					this.SweepText.setText("Current sweep: " + ((this.sweep+1)/2));
				}
			}
		}
	}
	
	public double[] GetNextWF()
	{
		double [] Point = new double[2];
		Point[0] = -1;
		Point[1] = -1;
		if (((this.WFRead == -1 && this.EnergyRead == -1) || (this.WFRead == 1 && this.EnergyRead == 1)))
		{	
			if (this.lock == 0)
			{
				this.lock = 1;
				this.DoNextCalc();
				this.OldY = this.y;
				this.OldWf = this.wf;
				this.OldEn = this.energy;
				this.WFRead = 0;
				this.EnergyRead = 0;
				this.lock = 0;
			}
		}
		if (this.WFRead == 0)
		{
			Point[0] = this.y;
			Point[1] = this.wf;
			this.WFRead = 1;
		}
		return Point;
	}
	
	public double[] GetNextEnergy()
	{
		double [] Point = new double[2];
		Point[0] = -1;
		Point[1] = -1;
		if (this.EnergyRead == 0)
		{
			if (this.lock == 0)
			{
				this.lock = 1;
				Point[0] = this.y;
				Point[1] = this.energy;
				if (this.energy < this.MinEnergy && this.sweep > 0)
					this.EnergyText.setText("Minimal energy: " + (this.MinEnergy=this.energy));
				this.EnergyRead = 1;
				this.lock = 0;
			}
		}
		return Point;
	}
	
	public void CreateExactSolution()
	{
		double norming = 1;
		if (this.potentialChooser == 0)
			norming = Math.sqrt(this.length/2.0-1.0/4.0*(this.length+1.0)/Math.PI*Math.sin(2.0*this.length*Math.PI/(this.length+1.0)));	// Normalisierungskonstante der exakten Wellenfunktion für den Fall ohne Potential
		if (this.potentialChooser == 1)
			norming = Math.sqrt(Math.sqrt(Math.PI)/((this.h*(this.length+1))/(this.length)));
		for (int j=0; j<this.length;j++)
		{
			if (this.potentialChooser == 0)
			{	
				this.exWFTrace.addPoint(2.0*(((double)j)/(this.length-1.0)-0.5)*this.size, Math.sin(j*Math.PI/(this.length-1))/norming);
				this.exEnTrace.addPoint(2.0*(((double)j)/(this.length-1.0)-0.5)*this.size, (2*(1.0-Math.cos(((double)(Math.PI/(this.length+1.0))))))/Math.pow(this.h, 2.0));
			}
			if (this.potentialChooser == 1)
			{
				this.exWFTrace.addPoint(2.0*(((double)j)/(this.length-1.0)-0.5)*this.size, Math.exp(-0.5*Math.pow(this.h, 2.0)*Math.pow(j-(this.length-1)/2.0, 2.0))/norming);
				this.exEnTrace.addPoint(2.0*(((double)j)/(this.length-1.0)-0.5)*this.size, 1.0);
			}
		}
		//System.out.println((2*(1.0-Math.cos(((double)(Math.PI/(this.length+1.0))))))/Math.pow(this.h, 2.0));
	}

	public void SetExWFTrace(Trace2DReplacing Trace)
	{
		this.exWFTrace = Trace;
	}

	public Trace2DReplacing GetExWFTrace()
	{
		return this.exWFTrace;
	}
	
	public void SetExEnTrace(Trace2DReplacing Trace)
	{
		this.exEnTrace = Trace;
	}
	
	public Trace2DReplacing GetExEnTrace()
	{
		return this.exEnTrace;
	}
	
	public double V(int j)
	{
		double V=0.0;
		if (this.potentialChooser == 1)			// Harmonisch (derzeit Doppelmulde)
			V = Math.pow(2.0*(((double)j)/(this.length-1)-0.5)*this.size,2.0);
		if (this.potentialChooser == 2)			// Anharmonisch
			V = Math.pow(2.0*(((double)j)/(this.length-1)-0.5)*this.size,2.0)+10.0*Math.pow(2.0*(((double)j)/(this.length-1)-0.5)*this.size,4.0);
		if (this.potentialChooser == 3)			// Doppelmulden
			V = -Math.pow(2.0*(((double)j)/(this.length-1)-0.5)*this.size,2.0)+0.02*Math.pow(2.0*(((double)j)/(this.length-1)-0.5)*this.size,4.0);
		if (this.potentialChooser == 4)			// Rechteck
			if (j == this.length/2.0 || j == this.length/2.0+1)
				V = -0.5;
		return V;

	}
}
