TDVPData.m

classdef TDVPData
	% The 'TDVPData' class is used to load, process and analyze data
	% generated by VMPS/TDVP code
    %
	
	% Defines all properties of objects from class 'TDVPData'.
    % They are accessed by objName.property
    properties
		tCalc;		% t of the simulation
		dt;			% timesteps
		t;			% t of observables
		tS;			% t of star extraction
		omega;		% frequency axis for star picture
		para;
		tresults;
		version;
		nChains;
		chainLabel; % Name of the Chains
		sysLabel;	% Name of System States
		% for SBM:
		alpha;		% if Spin-Boson
		s;			% if Spin-Boson
		spin;		% Spin Observable
		% for DPMES
		CTShift;	% string
		% for all
		lastIdx;    % last Observable Idx
		occC;		% chain occupation
		occCd;		% chain occupation, diabatic projection
		occCa;		% chain occupation, adiabatic projection
		occCc;		% chain occupation, |TT><LE+| coherence projection
		occS;		% star occupation
		xC;			% chain displacement
		xCd;		% chain displacement, diabatic projection
		xCa;		% chain displacement, adiabatic projection
		xCc;		% chain displacement, |TT><LE+| coherence projection
		xC2;		% chain <x^2>
		xC2d;		% chain <x^2>, diabatic projection
		xC2a;		% chain <x^2>, adiabatic projection
		xS;			% star displacement
		jC;			% chain current
		rho;		% rduced density matrix
		rhoOscRes;	% residual of populations after exponential fit
		stateProj;	% state projection amplitude
		sysState;	% state of system
		Heff;		% effective potential for system
        mps;        % MPS and Vmat of sites 1&2
        Vmat;
		treeMPS;	% The TreeMPS
		LegLabel;	%
		Comment;	% 
		folder;		% for easy identification
	end
	
	methods
			
		
		function obj = TDVPData(fname,varargin)
			% obj = TDVPData(para,tresults)
			%
			% Load TDVP dataset into a TDVPData object
			if nargin > 1
				% TDVPData(para,tresults)
				% fname = para, varargin{1} = tresults
				if isstruct(fname) && isstruct(varargin{1})
					temp.para = fname;
					temp.tresults = varargin{1};
				else
					fprintf('please program me\n');
				end
			elseif nargin == 0
				return;		% empty initialization
			elseif nargin == 1 && ischar(fname)
				try
					temp = load(fname);
					% Deserialize if needed
					Vars = fields(temp);
					for ii = 1:size(Vars)
						x = temp.(Vars{ii});
						if isa(x,'uint8')
							temp.(Vars{ii}) = hlp_deserialize(x);
	% 						eval(sprintf('%s = hlp_deserialize(%s);',Vars(ii).name,Vars(ii).name));
						end
					end
				catch err
					fprintf('load failed for: %s\n',fname);
					return;
				end
			elseif iscell(fname)
				% TDVPData({file1,file2,file3,...}
				% to be loaded into TDVPData array!
				L = length(fname);
				obj(L,1) = TDVPData();
				for ii = 1:L
					obj(ii) = TDVPData(fname{ii});
				end
				return;
			end
			
			% assign parameters
			if isfield(temp, 'para') && isfield(temp,'tresults')
				obj.para = temp.para;
				obj.tresults = temp.tresults;
			else
				return
			end
			if isfield(obj.para.tdvp,'version')
				obj.version = str2double(obj.para.tdvp.version(2:end));
			else
				obj.version = 0;
			end
			
			if isfield(obj.tresults,'lastIdx')
				obj.lastIdx = obj.tresults.lastIdx;
			else
				obj.lastIdx = size(obj.tresults.nx,1);
			end
			
			if isfield(obj.tresults,'spin')
				obj.spin = zeros(obj.lastIdx,3,'single');				% t x 3
				obj.spin(:,1) = obj.tresults.spin.sx(1:obj.lastIdx);
				obj.spin(:,2) = obj.tresults.spin.sy(1:obj.lastIdx);
				obj.spin(:,3) = obj.tresults.spin.sz(1:obj.lastIdx);
			end
			
			if obj.version >= 50 && isfield(obj.tresults,'n')
				obj.occC = obj.tresults.n;
			elseif isfield(obj.tresults,'nx')
				obj.occC = real(obj.tresults.nx);
			end
			if obj.version >= 72 && isfield(obj.tresults,'na')
				obj.occCa = obj.tresults.na;
				if isempty(obj.occC)
					obj.occC = squeeze(sum(real(obj.occCa),3));
				end
			end
			if obj.version >= 72 && isfield(obj.tresults,'nd')
				obj.occCd = obj.tresults.nd;
				if isempty(obj.occC)
					obj.occC = squeeze(sum(real(obj.occCd),3));
				end
			end
			if obj.version >= 72 && isfield(obj.tresults,'nc')
				obj.occCc = obj.tresults.nc;
			end
			
			if obj.version >= 66 && isfield(obj.tresults,'x')
				if ndims(obj.tresults.x) == 3
					obj.xC = obj.tresults.x;
				elseif ~isempty(strfind(obj.para.tdvp.Observables,'.x.'))		% support for version <=72
					obj.xCd = obj.tresults.x;
					if isempty(obj.xC)
						obj.xC = squeeze(sum(real(obj.xCd),3));
					end
				elseif ~isempty(strfind(obj.para.tdvp.Observables,'.x2.'))		% support for version <=72
					obj.xCa = obj.tresults.x;
					if isempty(obj.xC)
						obj.xC = squeeze(sum(real(obj.xCd),3));
					end
				end
			end
			if obj.version >= 72 && isfield(obj.tresults,'xa')
				obj.xCa = obj.tresults.xa;
				if isempty(obj.xC)
					obj.xC = squeeze(sum(real(obj.xCa),3));
				end
			end
			if obj.version >= 72 && isfield(obj.tresults,'xd')
				obj.xCd = obj.tresults.xd;
				if isempty(obj.xC)
					obj.xC = squeeze(sum(real(obj.xCd),3));
				end
			end
			if obj.version >= 72 && isfield(obj.tresults,'xc')
				obj.xCc = obj.tresults.xc;
			end
			if obj.version >= 72 && isfield(obj.tresults,'x2')
				obj.xC2 = obj.tresults.x2;
			end
			if obj.version >= 72 && isfield(obj.tresults,'x2a')
				obj.xC2a = obj.tresults.x2a;
				if isempty(obj.xC2)
					obj.xC2 = squeeze(sum(real(obj.xC2a),3));
				end
			end
			if obj.version >= 72 && isfield(obj.tresults,'x2d')
				obj.xC2d = obj.tresults.x2d;
				if isempty(obj.xC2)
					obj.xC2 = squeeze(sum(real(obj.xC2d),3));
				end
			end
			
			if obj.version >= 42 && isfield(obj.tresults,'t')
				obj.t = obj.tresults.t;
				obj.dt    = max(diff(obj.tresults.t));				% important for fourier smoothing
			else
				obj.t = obj.para.tdvp.t;
				obj.dt    = max(diff(obj.t));				% important for fourier smoothing
			end
			obj.tCalc = obj.para.tdvp.t;
			
			
			if isfield(obj.para,'nChains')
				obj.nChains = obj.para.nChains;
			else
				obj.nChains = 1;
			end
			
			if isfield(obj.tresults,'rho')
				obj.rho = obj.tresults.rho;
			else
				obj.rho = [];
			end
			obj.rhoOscRes = [];
			
			if isfield(obj.tresults,'system') && isfield(obj.tresults.system,'state')
				obj.sysState = obj.tresults.system.state;
			else
				obj.sysState = [];
			end
			
			
			% save star observables
			if isfield(obj.tresults,'star')
				if isfield(obj.tresults.star,'n')
					obj.occS  = obj.tresults.star.n;
					if obj.version <= 42 && size(obj.occS,1) == length(obj.tresults.star.omega)
						obj.occS = obj.occS.';
					end
				else
					obj.occS = 0;
				end
				if isfield(obj.tresults.star,'x')
					obj.xS    = obj.tresults.star.x;
				else
					obj.xS = 0;
				end
				obj.tS    = obj.tresults.star.t;
				if size(obj.tresults.star.omega,1) ~= 1
					obj.omega = obj.tresults.star.omega;		% w x nChains
				else
					obj.omega = obj.tresults.star.omega.';		% w x 1
				end
			end
			
			% extract simulation parameters SBM
			if ~isempty(strfind(obj.para.model,'SpinBoson'))
				if isa(obj.para.chain,'struct')
					if isfield(obj.para,'alpha')
						obj.alpha = obj.para.alpha;
					end
					if isfield(obj.para,'s')
						obj.s = obj.para.s;
					end

				elseif isa(obj.para.chain,'cell')
					obj.alpha = cellfun(@(x) x.alpha, obj.para.chain);
					obj.s     = cellfun(@(x) x.s    , obj.para.chain);
				end
			end
			
			if isfield(obj.tresults,'stateProjection')
				obj.stateProj = obj.tresults.stateProjection;
			end
			
			if isfield(obj.tresults,'Heff')
				obj.Heff = obj.tresults.Heff;
			end
			
			if isfield(obj.tresults,'mps')
                obj.mps = obj.tresults.mps;
                obj.tresults.mps = [];
                obj.Vmat = obj.tresults.Vmat;
                obj.tresults.Vmat = [];
			end
			
			if isfield(temp,'treeMPS')
				obj.treeMPS = temp.treeMPS;
			end
			
			obj.chainLabel = arrayfun(@(x) sprintf('$%d$',x),1:obj.nChains,'UniformOutput',false);
			obj.sysLabel   = arrayfun(@(x) sprintf('$%d$',x),1:2,'UniformOutput',false);			% TODO: needs fixing
			if ~isempty(strfind(obj.para.model,'DPMES'))
				% Need fully enclosed Latex expressions to enable embedding into other expressions.
				obj.sysLabel = {'$\mathrm{TT}$','$\mathrm{LE}^+$','$\mathrm{LE}^-$','$\mathrm{CT}^+$','$\mathrm{CT}^-$'};
				% any particular vars for DPMES?
				temp = regexp(obj.para.folder,'CT([-.0-9]*)','tokens');
				if ~isempty(temp)
					obj.CTShift = temp{1}{1};
				else
					temp = regexp(obj.para.folder,'Delta([-.0-9]*)State','tokens');
					if ~isempty(temp)
						obj.CTShift = temp{1}{1};
					else
						obj.CTShift = 0;	% give up
					end
				end
				if strcmp(obj.para.model,'DPMES4-5C')
					obj.sysLabel   = {'$\mathrm{TT}$','$\mathrm{LE}^+$','$\mathrm{CT}^+$','$\mathrm{CT}^-$'};
					obj.chainLabel = {	'$A_{1,1}$',	'$A_{1,2}$',	'$B_{1}$',	'$A_2$',	'$B_{2}$'	};
				end
				if strcmp(obj.para.model,'DPMES5-7C')
					obj.chainLabel = {	'$A_{1,1}$',	'$A_{1,2}$',	'$A_2$',	'$B_{1}$',	'$B_{2,1}$',  '$B_{2,2}$',  '$B_{2,3}$'	};
				end
				if strcmp(obj.para.model,'DPMESclust7-1')
					obj.chainLabel = {	'$B_{1,1}$',	'$B_{1,2}$',	'$A_{1,1}$',	'$A_{1,2}$',	'$B_{2,1}$',  '$B_{2,2}$',  '$A_2$'	};
				end
				if strcmp(obj.para.model,'DPMES-Tree4')
					obj.chainLabel = {	'$B_{1,1}$',	'$B_{1,2}$',	'$A_{1,1}$',	'$A_{1,2}$',	'$A_2$',	'$B_{2,1}$',  '$B_{2,2}$'	};
				end
				if strcmp(obj.para.model,'DPMES-Tree6')
					%	'B21',	'B22',	'B11',	'B12',	'B13',	'B14',  'A2',	'A11',	'A12'
					obj.chainLabel = {	'$B_{2,1}$',	'$B_{2,2}$',	'$B_{1,1}$',	'$B_{1,2}$', '$B_{1,3}$',	'$B_{1,4}$',	'$A_2$',	'$A_{1,1}$',  '$A_{1,2}$'	};
				end
			end
			obj.LegLabel = '';
			obj.folder = obj.para.folder;
		end
		
		function obj = setLegLabel(obj,entry)
			if iscell(entry)
				for i = 1:length(obj)
					obj(i) = obj(i).setLegLabel(entry{i});
				end
				return;
			end
			if ~isa(entry,'string')
				obj.LegLabel = num2str(entry);
			else
				obj.LegLabel = varName;
			end
		end
		
		function obj = setComment(obj, str)
			obj.Comment = str;
		end
		
		function obj = setHeffTo(obj,type)
			%% obj = setHeffTo(obj,type)
			% type: {'pes','tes'}	to select between TES and PES in the plotting routines
			%
			switch type
				case 'pes'
					obj.Heff = obj.tresults.pes;
				case 'tes'
					obj.Heff = obj.tresults.Heff;
			end
		end
		
		function out = gettRhoiiSystem(obj)
			% DEPRECATED: Use getData('rhoii') instead
			% returns the diagonal of the system reduced density matrix
			out = 0;
			if isfield(obj.tresults,'PPCWavefunction')
				out = obj.tresults.PPCWavefunction;
			elseif isfield(obj.tresults,'rho')
				out = cell2mat(arrayfun(@(i) diag(squeeze(obj.tresults.rho(i,:,:))),1:size(obj.tresults.rho,1),'UniformOutput',false)).';
			end
		end
		
		function out = getData(obj,type,varargin)
			%% out = getData(obj,type,varargin)
			%	returns certain data or observables:
			%
			%	getData('rhoii')				t x dk system RDM populations
			%	getData('rhoii','state',n)		t x 1  system RDM population of state n
			%
			full = 0;		% truncate to simulated time - tresults.lastIdx
			idxOffset = 0;
			out = 0;
						
			% Parse additional inputs
			p = inputParser;
			addParameter(p,'state',0,@isnumeric);		% used in heff
			addParameter(p,'full' ,0,@islogical);
			p.addOptional('factor',0,@isnumeric);
			parse(p,varargin{:});
			
			switch lower(type)
				case 'calctime'
					out = obj.para.tdvp.calcTime;
					idxOffset = 1;
				case 'calctime-d'
					out = diff(obj.para.tdvp.calcTime);
					idxOffset = 2;
				case 'rhoii'
					if isfield(obj.tresults,'PPCWavefunction')
						rhoii = obj.tresults.PPCWavefunction;
					elseif isfield(obj.tresults,'rho')
						rhoii = cell2mat(arrayfun(@(i) diag(squeeze(obj.tresults.rho(i,:,:))),1:obj.lastIdx,'UniformOutput',false)).';
					else
						warning('TDVPData:getData:rhoii','rho does not exist');
						rhoii = 0;
					end
					out = rhoii;
					h.state = p.Results.state;
					if h.state ~= 0 && h.state <= size(out,2)
						% return only a single state
						out = out(:,h.state);
					end
					full = 1;
				case 'rhoii-osc-res'
					% oscillating residuals
					plotFit = 1;						% 1: plot the fitted functions into separate figure to evaluate fitness.
					rhoii = abs(obj.getData('rhoii'));		% t x dk(1)
					out = zeros(size(rhoii));
					% fit and remove exponential components to obtain the residuals
					if plotFit
						f = gcf;
						figure; hold all; ax = gca;
						plot(obj.t(1:obj.lastIdx),rhoii);
						ax.ColorOrderIndex = 1;					% reset to have fit in same color as original data
					end
					for m = 1:size(rhoii,2)
						fun{1} = @(x,xdata) x(1)+x(2).*exp(x(3).*xdata);					% exponential model with shift
						fun{2} = @(x,xdata) x(1)+x(2).*exp(x(3).*xdata)+x(4).*exp(x(5).*xdata);					% exponential model with shift
						a0{1} = double([rhoii(end,m),rhoii(1,m)-rhoii(end,m),-0.01]);								% initial guesses
						if rhoii(1,m) == 0
							a0{1}(3) = mean(diff(rhoii(100:200,m)))/obj.t(2);			% rise if starting from 0
						end
						a0{2} = double([rhoii(end,m),rhoii(1,m)-rhoii(end,m),-0.01,rhoii(1,m)-rhoii(end,m),+0.01]);								% initial guesses
						for k = 1:length(fun)
							[x{k},r(k),res{k}] = lsqcurvefit(fun{k},a0{k},double(obj.t(1:obj.lastIdx).'), double(rhoii(:,m)));
						end
						[~,ind] = min(r);		% pick function with minimal residual norm;
						if plotFit, plot(obj.t(1:obj.lastIdx),fun{ind}(x{ind},obj.t(1:obj.lastIdx))); end
						out(:,m) = res{ind};		% save the residual
					end
					if plotFit,	figure(f); end
					full = 1;
				case 'rhoii-osc-res-med'
					% oscillating residuals
					rhoii = abs(obj.getData('rhoii'));		% t x dk(1)
					% fit by smoothing and remove exponential components to obtain the residuals
					out = TDVPData.movAvgRes(rhoii,170/obj.t(2));		% 170 or 350
					full = 1;
				case 'rhoii-osc-res-mean'
					% oscillating residuals by removing DC via mean
					out = TDVPData.meanRes(abs(obj.getData('rhoii')));
					full = 1;
				case 'rhoij'
					% out{1}: t x #off-diagonal elements
					% out{2}: [from,to] x #off-diagonal elements
					%			carries information to identify objects in out{1}
					if ~isfield(obj.tresults,'rho') || isempty(obj.tresults.rho)
						error('TDVPData:getData:rhoij','rho does not exist')
					end
					out = {};
					r = obj.tresults.rho(1:obj.lastIdx,:,:);
					d = size(r);
					out{1} = zeros(d(1),d(2)*(d(2)-1)/2-1);				% calculates number of off-diagonals
					out{2} = zeros(2,size(out,2));
					iC = 1;
					for jj = 1:d(2)
						for ii = jj:d(2)
							if ii ~= jj
								out{1}(:,iC) = r(:,ii,jj);				% get offdiag
								out{2}(:,iC) = [jj,ii];
								iC = iC+1;
							end
						end
					end
					full = 1;
				case 'vne'
					% von Neumann Entropy = -tr(r ln r)
					rdm = obj.tresults.rho;						% t x i x j
					out = zeros(size(rdm,1),1);					% t x 1
					rdm = permute(rdm,[2,3,1]);					% i x j x t
					for ii = 1:length(out)
						out(ii) = real(-sum(sum(rdm(:,:,ii).*log(rdm(:,:,ii).'))));
					end
                case 'sys-env-x'
					out = obj.getSysEnvObs('x');			% returns t x NC cell array; using mps(:,2), Vmat(:,2)
					full = 1;
                case 'sys-env-n'
					out = obj.getSysEnvObs('n');			% returns t x NC cell array; using mps(:,2), Vmat(:,2)
                    full = 1;
				case 'sz'
					out = obj.tresults.spin.sz;
				case 'spin'
				case 'spin-ttm'
					out = struct();
					[sx,sy,sz] = spinop('Z');
					finalT = obj.t(end)*p.Results.factor;					% how far to extrapolate in time, TODO: make argument
					n = round(finalT/obj.dt);
					rhoT = zeros(length(obj.tresults.TTM.T)*4,1);
					Esigma = zeros(n,3);
					T = reshape(obj.tresults.TTM.T, 4,[]);			% creates [T(1) T(2) T(3) ...]
					for i = 1:n
						if i == 1
							rho = [1,0,0,0]';
						else
							rho = T*rhoT;
						end
						rhoT = [rho; rhoT(1:end-4)];					% prepend new vector rho(i)
						rho = reshape(rho,[2,2]);						% reshape rho(i) for observables
						Esigma(i,1) = trace(sx*rho);
						Esigma(i,2) = trace(sy*rho);
						Esigma(i,3) = trace(sz*rho);
					end
					out.sx = real(Esigma(:,1));
					out.sy = real(Esigma(:,2));
					out.sz = real(Esigma(:,3));
					out.t  = 0:obj.dt:finalT;
					return;
				case 'heff-e'
					% calculate energy from heff & sysState
					d = size(obj.Heff);																% t x D x dk x D x dk
					obj.Heff = permute(reshape(obj.Heff, [d(1),d(2)*d(3), d(4)*d(5)]), [2,3,1]);	% (D * dk) x (D * dk) x t
					obj.sysState = reshape(permute(obj.sysState, [3,2,1]), [d(2)*d(3),d(1)]);		% (D * dk) x t
					out = arrayfun(@(i) (obj.sysState(:,i)'*obj.Heff(:,:,i)*obj.sysState(:,i))/(obj.sysState(:,i)'*obj.sysState(:,i)),1:obj.lastIdx);
					out = real(out);
				case 'heff'
					assert(~isempty(obj.Heff),'Heff was not extracted in simulation');
					h.state = p.Results.state;
					[V,D] = arrayfun(@(i) eig(squeeze(obj.Heff(i,h.state,:,h.state,:)),'vector'),[1:obj.lastIdx]','UniformOutput',false);
					D = real(cell2mat(D'))';
					[D,I] = sort(D,2);																% sort eigenvalues ascending, get I to sort V
					V = arrayfun(@(i) V{i}(:,I(i,:)),[1:obj.lastIdx]','UniformOutput',false);		% reorder eigenvectors accordingly
					out = {};
					out{1} = D;																		% t x dk_eig
					out{2} = V;																		% t x dk x dk_eig
				case 'heff-pop'
					h.state = p.Results.state;
					temp = obj.getData('heff','state',h.state);
					D = temp{1}; V = temp{2};
					% V{i}: dk x D_eig
					Vtemp = cell2mat(V);												% creates (dk*t) x D_eig
					Vtemp = reshape(Vtemp,[size(V{1},1),length(V),size(V{1},1)]);		% dk x t x D_eig
					Vtemp = permute(Vtemp, [2,3,1]);									% t x D_eig x dk

					% obj.sysState: t x dk x D
					tempState = permute(obj.sysState(1:obj.lastIdx,:,h.state),[1,3,2]);				% t x 1 x dk
					pop = sum(bsxfun(@times,Vtemp ,tempState),3);
					pop = pop.*conj(pop);															% t x dk_eig x 1
					out = {};
					out{1} = D;																		% t x dk_eig
					out{2} = pop;																	% t x dk_eig
				case 'heff-pop-diab'
					h.state = p.Results.state;
					temp = obj.getData('heff','state',h.state);
					D = temp{1}; V = temp{2};
					% V{i}: dk x D_eig
					Vtemp = cell2mat(V);												% creates (dk*t) x D_eig
					Vtemp = reshape(Vtemp,[size(V{1},1),length(V),size(V{1},1)]);		% dk x t x D_eig
					Vtemp = permute(Vtemp, [2,3,1]);									% t x D_eig x dk

					% obj.sysState: t x dk x D
					tempState = permute(obj.sysState(1:obj.lastIdx,:,h.state),[1,3,2]);				% t x 1 x dk
					pop = sum(bsxfun(@times,Vtemp ,tempState),3);
					pop = pop.*conj(pop);															% t x dk_eig x 1

					pop = bsxfun(@times,(Vtemp.*conj(Vtemp)),pop);									% t x dk_eig x dk
					pop = permute(pop,[1,3,2]);
					out = {};
					out{1} = D;								% t x dk_eig
					out{2} = pop;							% t x dk x dk_eig
				case 'heff-swap'
					% swaps Bond states according to matrix defined in 'state':
					% h.state has to be a matrix defining: [tIdx, n1, n2], which pair of D states n1, n2 have to be swapped.
					h.state = p.Results.state;
					if ~ismatrix(h.state)
						error('for this, h.state has to be a matrix defining: [tIdx, n1, n2], which D states have to be swapped');
					end
					obj2 = obj;		% copy obj.
					D = size(obj.sysState,3);
					for ii = 0:size(h.state,1)
						if ii == 0
							start = 1;
							stop = h.state(1,1)-1;
							thisPerm = eye(D);										% start off with identity matrix
						else
							start = stop+1;
							if ii < size(h.state,1)
								stop  = h.state(ii+1,1)-1;												% last index within this cut, should always be < obj.lastIdx
							else
								stop = obj.lastIdx;
							end
							newPerm = eye(D);
							newPerm([h.state(ii,2),h.state(ii,3)],:) = newPerm([h.state(ii,3),h.state(ii,2)],:);
							thisPerm = thisPerm*newPerm;
						end
					
						% apply permutation
						obj2.sysState(start:stop,:,:) = contracttensors(obj.sysState(start:stop,:,:),3,3, thisPerm,2,1);
						
						temp = contracttensors(obj.Heff(start:stop,:,:,:,:),5,2, thisPerm,2,1);			% t x dk x D x dk x D,
						temp = contracttensors(temp, 5,3, thisPerm,2,1);								% t x dk x dk x D x D,
						obj2.Heff(start:stop,:,:,:,:) = permute(temp,[1,4,2,5,3]);						% t x D x dk x D x dk
					end
					out = obj2;
				case 'heff-current'
					% best on an already swapped Heff.
					% calculates the transition rates / currents from one diabatic state to another
					% is equivalent to evaluating all the elements which can make up the current operator j = i[H,N]
					% out: t x D' x dk' x D x dk
					rho_SE = bsxfun(@times, permute(conj(obj.sysState), [1 3 2]) , permute(obj.sysState, [1,4,5,3,2]));  % t x D' x dk' x 1 x 1 .* t x 1 x 1 x D x dk  = t x D' x dk' x D x dk
					
					out = -2 * imag(obj.Heff(1:obj.lastIdx,:,:,:,:) .* conj(rho_SE(1:obj.lastIdx,:,:,:,:)));		% is 2*Im( Tr(H rho*) )
					full = 1;
				case 'heff-current-squared'
					% Possibly relevant for superexchange. Uses H^2 
					% calculates the transition rates / currents from one diabatic state to another
					% is equivalent to evaluating all the elements which can make up the current operator j = i[H,N]
					% out: t x D' x dk' x D x dk
					rho_SE = bsxfun(@times, permute(conj(obj.sysState), [1 3 2]) , permute(obj.sysState, [1,4,5,3,2]));  % t x D' x dk' x 1 x 1 .* t x 1 x 1 x D x dk  = t x D' x dk' x D x dk
					
					heffsquared = squeeze(sum(sum(bsxfun(@times, obj.Heff, permute(obj.Heff, [1,6,7,2,3,4,5])),4),5)); 
					out = -2 * imag(heffsquared(1:obj.lastIdx,:,:,:,:) .* conj(rho_SE(1:obj.lastIdx,:,:,:,:)));		% is 2*Im( Tr(H rho*) )
					full = 1;
				case 'heff-current-v2'
					% calculates the transition rates / currents from one diabatic state to another
					% is equivalent to evaluating all the elements which can make up the current operator j = i[H,N]
					rho_SE = bsxfun(@times, permute(conj(obj.sysState), [1 3 2]) , permute(obj.sysState, [1,4,5,3,2]));  % t x D' x dk' x 1 x 1 .* t x 1 x 1 x D x dk  = t x D' x dk' x D x dk
					rho_SE = rho_SE(1:obj.lastIdx,:,:,:,:);
					
					d = size(rho_SE);
					Nop = zeros(d(5)); Nop(2,2) = 1;												% select which state to look at.
					HN = contracttensors(obj.Heff(1:obj.lastIdx,:,:,:,:),5,5,Nop,2,1);				% multiply into last dk
					
					out = -2 * imag(HN .* conj(rho_SE));		% is 2*Im( Tr(H rho*) )
					full = 1;
				case 'heff-diab-state'
					% Transform Heff with the Eigensystem of the LE+ staets only.
					% This eliminates the coupling terms for different environment states but same system states, assuming that the trafo is similar for all system states.
					% only cross-terms between different electronic states are still present
					% Find and Sort EigenSystem of the LE+ Heff
					% DEPRECATED: USE heff-full-diab-partdiag
					assert(~isempty(obj.Heff),'Heff was not extracted in simulation');
					d = size(obj.Heff);																% t x D' x dk' x D x dk
					obj.Heff = reshape(obj.Heff,d(1),d(2)*d(3),d(4)*d(5));							% t x D' * dk' x D * dk
					% heff is 'blockdiagonal' with dk running faster than D
					sub_idx = (1:d(2)) + (p.Results.state-1)*d(2);
					[V,D] = arrayfun(@(i) eig(squeeze(obj.Heff(i,sub_idx,sub_idx)),'vector'),[1:obj.lastIdx]','UniformOutput',false);
					D = real(cell2mat(D'))';
					[D,I] = sort(D,2);																% sort eigenvalues ascending, get I to sort V
					V = arrayfun(@(i) kron(eye(d(2)),V{i}(:,I(i,:))),[1:obj.lastIdx]','UniformOutput',false);		% reorder eigenvectors accordingly and expand to D*dk
					
					% Transform Heff to remove coupling between same system states
					heff = arrayfun(@(i) V{i}'*squeeze(obj.Heff(i,:,:))*V{i}, [1:obj.lastIdx]','UniformOutput',false);
					D = arrayfun(@(i) diag(squeeze(heff{i})),1:obj.lastIdx,'UniformOutput',false);
					D = real(cell2mat(D))';
					
					% VC for V-contribution/ diabatic Character 
					VC = cell2mat(V);																% creates (D*dk*t) x D*dk_eig
					VC = reshape(VC,[size(V{1},1),length(V),size(V{1},1)]);							% D*dk x t x D*dk_eig
					VC = permute(VC, [3,1,2]);														% D*dk_eig x D*dk x t
					
					d = size(VC);
					VC = reshape(VC, [d(1),sqrt(d(2))*[1,1], d(3)]);								% D*dk_eig x D x dk x t
					VC = squeeze(sum(VC.*conj(VC),2));												% D*dk_eig x dk x t
					% Now: VC has percentage of dk contribution on each D*dk_eig and t
					
					out = {};
					out{1} = D;																		% t x D*dk_eig
					out{2} = V;																		% t x D*dk x D*dk_eig, cell
					out{3} = VC;																	% D*dk_eig x dk x t
				case 'heff-full-nondiag'
					%% Do not Diagonalise Heff
					% return heff-full diagonal entries without diagonalisation
					% do not sort, since this allows using plain lines with correct colors
					assert(~isempty(obj.Heff),'Heff was not extracted in simulation');
					d = size(obj.Heff);																% t x D' x dk' x D x dk
					heff = reshape(obj.Heff,d(1),d(2)*d(3),d(4)*d(5));								% t x D' * dk' x D * dk
					D = arrayfun(@(i) diag(squeeze(heff(i,:,:)))',[1:obj.lastIdx]','UniformOutput',false);
					D = real(cell2mat(D));															% t x D*dk
% 					[D,I] = sort(D,2);																% sort eigenvalues ascending, get I to sort V
					
					V = arrayfun(@(i) eye(d(2)*d(3)),[1:obj.lastIdx]','UniformOutput',false);		% identity since no basis trafo! compatible with other functions!
% 					V = arrayfun(@(i) V{i}(:,I(i,:)),[1:obj.lastIdx]','UniformOutput',false);		% Reorder to energy sorting
					
					%% Make output
					out = {};
					out{1} = D;																		% t x D*dk
					out{2} = V;																		% t x D*dk x D*dk_eig, cell
				case 'heff-full-diab-nondiag'
					% Get diabatic character of surfaces from Heff trafo
					% VC: diabatic character normalised to 1
					out = obj.getData('heff-full-nondiag');
					D = out{1};																		% t x D*dk_eig
					V = out{2};																		% t x D*dk x D*dk_eig, cell
					
					% VC for V-contribution/ diabatic Character 
					VC = diabaticCharacter(V);														% D*dk_eig x dk x t
					% Now: VC has percentage of dk contribution on each D*dk_eig and t
					
					out = {};
					out{1} = D;																		% t x D*dk_eig
					out{2} = V;																		% t x D*dk x D*dk_eig, cell
					out{3} = VC;																	% D*dk_eig x dk x t
				case 'heff-full-pop-nondiag'
					%% return heff-full, state, population, and diabatic population without diagonalisation
					out = obj.getData('heff-full-nondiag');
					D = out{1};																		% t x D*dk_eig
					V = out{2};																		% t x D*dk x D*dk_eig, cell
					
					Vtemp = permute(cell2mat(permute(V,[2,3,1])),[3,1,2]);							% t x D*dk x D*dk_eig
					
					% obj.sysState: t x dk x D
					d = size(obj.sysState); d(1) = obj.lastIdx;
					tempState = permute(obj.sysState(1:obj.lastIdx,:,:),[1,3,2]);					% t x D x dk
					tempState = reshape(tempState, [d(1),d(2)*d(3)]);								% t x D*dk
					tempState = bsxfun(@times,conj(Vtemp),tempState);								% t x D*dk x D*dk_eig
					
					pop = squeeze(sum(tempState,2));
					pop = pop.*conj(pop);															% t x D*dk_eig
					% pop is sysState rotated into adiab Basis then taken probability -> should be correct
					
					% select popDiab2:
					Vpop = Vtemp.*conj(Vtemp);														% t x D*dk x D*dk_eig
					Vpop = reshape(Vpop,d(1), d(2), d(3), []);										% t x D x dk x D*dk_eig
					Vpop = squeeze(sum(Vpop, 2));													% t x dk x D*dk_eig
					popDiab = bsxfun(@times, Vpop, permute(pop, [1,3,2]));							% t x dk x D*dk_eig
					
					%% Make output
					out = {};
					out{1} = D;																		% t x D*dk
					out{2} = pop;																	% t x D*dk
					out{3} = popDiab;																% t x dk x D*dk_eig
					out{4} = Vtemp;																	% t x D*dk x D*dk_eig
				case 'heff-full-partdiag'
					% Transform Heff with the Eigensystem of the LE+ staets only.
					% This eliminates the coupling terms for different environment states but same system states, assuming that the trafo is similar for all system states.
					% only cross-terms between different electronic states are still present
					% Find and Sort EigenSystem of the LE+ Heff
					assert(~isempty(obj.Heff),'Heff was not extracted in simulation');
					d = size(obj.Heff);																% t x D' x dk' x D x dk
					obj.Heff = reshape(obj.Heff,d(1),d(2)*d(3),d(4)*d(5));							% t x D' * dk' x D * dk
					% heff is 'blockdiagonal' with dk running faster than D
					sub_idx = (1:d(2)) + (p.Results.state-1)*d(2);
					[V,D] = arrayfun(@(i) eig(squeeze(obj.Heff(i,sub_idx,sub_idx)),'vector'),[1:obj.lastIdx]','UniformOutput',false);
					D = real(cell2mat(D'))';
					[D,I] = sort(D,2);																% sort eigenvalues ascending, get I to sort V
					V = arrayfun(@(i) kron(eye(d(2)),V{i}(:,I(i,:))),[1:obj.lastIdx]','UniformOutput',false);		% reorder eigenvectors accordingly and expand to D*dk
					
					% Transform Heff to remove coupling between same system states
					heff = arrayfun(@(i) V{i}'*squeeze(obj.Heff(i,:,:))*V{i}, [1:obj.lastIdx]','UniformOutput',false);
					D = arrayfun(@(i) diag(squeeze(heff{i})),1:obj.lastIdx,'UniformOutput',false);
					D = real(cell2mat(D))';
					
					%% Make output
					out = {};
					out{1} = D;																		% t x D*dk
					out{2} = V;																		% t x D*dk x D*dk_eig, cell
				case 'heff-full-diab-partdiag'
					% Get diabatic character of surfaces from Heff trafo
					% VC: diabatic character normalised to 1
					out = obj.getData('heff-full-partdiag','state',p.Results.state);
					D = out{1};																		% t x D*dk_eig
					V = out{2};																		% t x D*dk x D*dk_eig, cell
					
					% VC for V-contribution/ diabatic Character 
					VC = diabaticCharacter(V);														% D*dk_eig x dk x t
					% Now: VC has percentage of dk contribution on each D*dk_eig and t
					
					out = {};
					out{1} = D;																		% t x D*dk_eig
					out{2} = V;																		% t x D*dk x D*dk_eig, cell
					out{3} = VC;																	% D*dk_eig x dk x t
				case 'heff-full-pop-partdiag'
					%% return heff-full, state, population, and diabatic population without diagonalisation
					out = obj.getData('heff-full-partdiag','state',p.Results.state);
					D = out{1};																		% t x D*dk_eig
					V = out{2};																		% t x D*dk x D*dk_eig, cell
					
					Vtemp = permute(cell2mat(permute(V,[2,3,1])),[3,1,2]);							% t x D*dk x D*dk_eig
					
					% obj.sysState: t x dk x D
					d = size(obj.sysState); d(1) = obj.lastIdx;
					tempState = permute(obj.sysState(1:obj.lastIdx,:,:),[1,3,2]);					% t x D x dk
					tempState = reshape(tempState, [d(1),d(2)*d(3)]);								% t x D*dk
					tempState = bsxfun(@times,conj(Vtemp),tempState);								% t x D*dk x D*dk_eig
					
					pop = squeeze(sum(tempState,2));
					pop = pop.*conj(pop);															% t x D*dk_eig
					% pop is sysState rotated into adiab Basis then taken probability -> should be correct
					
					% select popDiab2:
					Vpop = Vtemp.*conj(Vtemp);														% t x D*dk x D*dk_eig
					Vpop = reshape(Vpop,d(1), d(2), d(3), []);										% t x D x dk x D*dk_eig
					Vpop = squeeze(sum(Vpop, 2));													% t x dk x D*dk_eig
					popDiab = bsxfun(@times, Vpop, permute(pop, [1,3,2]));							% t x dk x D*dk_eig
					
					%% Make output
					out = {};
					out{1} = D;																		% t x D*dk
					out{2} = pop;																	% t x D*dk
					out{3} = popDiab;																% t x dk x D*dk_eig
					out{4} = Vtemp;																	% t x D*dk x D*dk_eig
				case 'heff-full-diag'
					% Diagonalise entire Heff
					% Find and Sort EigenSystem of the full Heff
					assert(~isempty(obj.Heff),'Heff was not extracted in simulation');
					d = size(obj.Heff);																% t x D' x dk' x D x dk
					obj.Heff = reshape(obj.Heff,d(1),d(2)*d(3),d(4)*d(5));							% t x D' * dk' x D * dk
					[V,D] = arrayfun(@(i) eig(squeeze(obj.Heff(i,:,:)),'vector'),[1:obj.lastIdx]','UniformOutput',false);
					D = real(cell2mat(D'))';
					[D,I] = sort(D,2);																% sort eigenvalues ascending, get I to sort V
					V = arrayfun(@(i) V{i}(:,I(i,:)),[1:obj.lastIdx]','UniformOutput',false);		% reorder eigenvectors accordingly
					out = {};
					out{1} = D;																		% t x D*dk_eig
					out{2} = V;																		% t x D*dk x D*dk_eig, cell
				case 'heff-full-diab'
					% 'heff-full-diab-diag'
					% VC: diabatic character normalised to 1
					out = obj.getData('heff-full-diag');
					D = out{1};																		% t x D*dk_eig
					V = out{2};																		% t x D*dk x D*dk_eig, cell
					
					% VC for V-contribution/ diabatic Character 
					VC = diabaticCharacter(V);														% D*dk_eig x dk x t
					% Now: VC has percentage of dk contribution on each D*dk_eig and t
					
					out = {};
					out{1} = D;																		% t x D*dk_eig
					out{2} = V;																		% t x D*dk x D*dk_eig, cell
					out{3} = VC;																	% D*dk_eig x dk x t
				case 'heff-full-pop'
					%% Calculates the Diabatic and Adiabatic population on the HEff surfaces
					% 'heff-full-pop-diag'
					%	Vtemp: trafo from diabat to adiabats
					temp = obj.getData('heff-full-diag');
					D = temp{1}; 																	% t x D*dk_eig
					V = temp{2};																	% t cell x D*dk x D*dk_eig
					
					Vtemp = permute(cell2mat(permute(V,[2,3,1])),[3,1,2]);							% t x D*dk x D*dk_eig
					
					% obj.sysState: t x dk x D
					d = size(obj.sysState); d(1) = obj.lastIdx;
					tempState = permute(obj.sysState(1:obj.lastIdx,:,:),[1,3,2]);					% t x D x dk
					tempState = reshape(tempState, [d(1),d(2)*d(3)]);								% t x D*dk
					tempState = bsxfun(@times,conj(Vtemp),tempState);								% t x D*dk x D*dk_eig
					
					pop = squeeze(sum(tempState,2));
					pop = pop.*conj(pop);															% t x D*dk_eig
					% pop is sysState rotated into adiab Basis then taken probability -> should be correct
					
					% select popDiab2:
					Vpop = Vtemp.*conj(Vtemp);														% t x D*dk x D*dk_eig
					Vpop = reshape(Vpop,d(1), d(2), d(3), []);										% t x D x dk x D*dk_eig
					Vpop = squeeze(sum(Vpop, 2));													% t x dk x D*dk_eig
					popDiab = bsxfun(@times, Vpop, permute(pop, [1,3,2]));							% t x dk x D*dk_eig
					
					out = {};
					out{1} = D;																		% t x D*dk_eig
					out{2} = pop;																	% t x D*dk_eig
					out{3} = popDiab;																% t x dk x D*dk_eig
					out{4} = Vtemp;																	% t x D*dk x D*dk_eig
					
					return;
					
					%% Different approaches to popDiab:
					% Summary:
					% 1 == 4: diabatic population is correct, but distribution over adiabats wrong
					% 2 == 5: adiabatic pop is correct, but division into diabats is wrong!
					% 3:	mixture between 1 and 2 but still wrong
					% can measure similarity with norm(popDiab(:)-popDiab4(:))
					
					% popDiab1: summing over D amplitudes before squaring for populations
					popDiab1 = reshape(tempState, d(1), d(3), d(2), []);							% t x D x dk x D*dk_eig
					popDiab1 = squeeze(sum(popDiab1,2));											% t x dk x D*dk_eig
					popDiab1 = popDiab1.*conj(popDiab1);
					
					% divide pop into popDiab by using mapping information from V
					% this yields results which do not respect inter-adiabat-coherences, 
					% thus are wrong when summing over D*dk_eig to recover diabat populations
					Vpop = Vtemp.*conj(Vtemp);														% t x D*dk x D*dk_eig
					Vpop = reshape(Vpop,d(1), d(2), d(3), []);										% t x D x dk x D*dk_eig
					Vpop = squeeze(sum(Vpop, 2));													% t x dk x D*dk_eig
					popDiab2 = bsxfun(@times, Vpop, permute(pop, [1,3,2]));							% t x dk x D*dk_eig
					
					% popDiab3: take pop and multiply with distribution of popDiab1
					temp = bsxfun(@rdivide, popDiab, sum(popDiab,2));	% normalise
					popDiab3 = bsxfun(@times, temp, permute(pop, [1,3,2]));							% t x dk x D*dk_eig
					
					% popDiab4: tempState^2 then sum over D; very similar to 1
					popDiab4 = reshape(tempState, d(1), d(3), d(2), []);							% t x D x dk x D*dk_eig
					popDiab4 = popDiab4.*conj(popDiab4);
					popDiab4 = squeeze(sum(popDiab4,2));											% t x dk x D*dk_eig
					
					% popDiab5: map from adiab into dk D and trace over D; identical to 2
					popDiab5 = sum(tempState,2);													% t x 1 x D*dk_eig
					popDiab5 = bsxfun(@times, Vtemp, popDiab5);										% t x D*dk x D*dk_eig
					popDiab5 = reshape(popDiab5, d(1), d(3), d(2), []);								% t x D x dk x D*dk_eig
					popDiab5 = squeeze(sum(popDiab5.*conj(popDiab5),2));							% t x dk x D*dk_eig
					
% 					out{3} = popDiab2;	% select between 1-5					
				case 'heff-full-main-bonds'
					% take the dominantly populated surface and derive its breakdown into bond states
					% this allows estimation of participating environments
					out = obj.getData('heff-full-pop');
					D		= out{1};																	% t x D*dk_eig
					pop		= out{2};																	% t x D*dk_eig
					popDiab = out{3};																	% t x dk x D*dk_eig
					Vtemp	= out{4};																	% t x D*dk x D*dk_eig
					
					% find which surface is dominantly populated
					[maxPop,Imax] = max(pop,[],2);
					% take corresponding mapping
					d = size(Vtemp);
					dk = size(popDiab,2);
					D = d(2)/dk;
					
					VtempNew = zeros(d(1),d(2));
					for ii = 1:d(1)
						VtempNew(ii,:) = Vtemp(ii,:,Imax(ii));											% t x D*dk
					end
					
					legLab = cell(d(2),1);
					for ii = 1:d(2)
						[mm,nn] = ind2sub([D,dk],ii);
						legLab{ii} = sprintf('$|$%s,$%d\\rangle$',obj.sysLabel{nn},mm);
					end
					
					out = {};
					out{1} = VtempNew;																	% t x D*dk
					out{2} = legLab;																	% D*dk
				case 'pes-from-tes'
					% get the adiabatic states from the TES
					% use these to get the corresponding PES
					
					%% First make sure TES is selected
					tes = obj.setHeffTo('tes');
					tesout = tes.getData('heff-full-diab');
					D = tesout{1};																	% t x D*dk_eig
					V = tesout{2};																	% t x D*dk x D*dk_eig, cell in timesteps
					VC = tesout{3};																	% D*dk_eig x dk x t
						% V contains map from diabatic System-Environment basis (SE) into adiabatic basis
					%% Get PES into heff
					pes = obj.setHeffTo('pes');
					d = size(pes.Heff);																% t x D' x dk' x D x dk
					pes.Heff = reshape(pes.Heff,d(1),d(2)*d(3),d(4)*d(5));							% t x D' * dk' x D * dk
					
					out = {};
					out{1} = arrayfun(@(ii) V{ii}'*squeeze(pes.Heff(ii,:,:))*V{ii}, [1:pes.lastIdx]','UniformOutput',false);	% t cell: D*dk_eig' x D*dk_eig
					% take only diagonal of potential:
					out{1} = cell2mat(cellfun(@(x) real(diag(x)'),out{1}, 'UniformOutput',false));	% t x D*dk_eig
					out{2} = VC;																	% D*dk_eig x dk x t
				case 'tes-from-pes'
					% get the adiabatic states from the PES
					% use these to get the corresponding TES
					
					%% First make sure TES is selected
					pes = obj.setHeffTo('pes');
					pesout = pes.getData('heff-full-diab');
					D = pesout{1};																	% t x D*dk_eig
					V = pesout{2};																	% t x D*dk x D*dk_eig, cell in timesteps
					VC = pesout{3};																	% D*dk_eig x dk x t
						% V contains map from diabatic System-Environment basis (SE) into adiabatic basis
					%% Get TES into heff
					tes = obj.setHeffTo('tes');
					d = size(tes.Heff);																% t x D' x dk' x D x dk
					tes.Heff = reshape(tes.Heff,d(1),d(2)*d(3),d(4)*d(5));							% t x D' * dk' x D * dk
					
					out = {};
					out{1} = arrayfun(@(ii) V{ii}'*squeeze(tes.Heff(ii,:,:))*V{ii}, [1:tes.lastIdx]','UniformOutput',false);	% t cell: D*dk_eig' x D*dk_eig
					% take only diagonal of potential:
					out{1} = cell2mat(cellfun(@(x) real(diag(x)'),out{1}, 'UniformOutput',false));	% t x D*dk_eig
					out{2} = VC;																	% D*dk_eig x dk x t
				case 'rho-se'
					% get rho from sysState which contains correlations with the environment.
					rho_SE = bsxfun(@times, permute(conj(obj.sysState), [1 3 2]) , permute(obj.sysState, [1,4,5,3,2]));  % t x D' x dk' x 1 x 1 .* t x 1 x 1 x D x dk  = t x D' x dk' x D x dk
					
					% make labels
					d = size(rho_SE);
					legLab = cell(prod(d(2:end)),1);
					for ii = 1:prod(d(2:end))
						[l,m,n,o] = ind2sub(d(2:end),ii);
						legLab{ii} = sprintf('$|%d\\rangle|%d\\rangle \\to |%d\\rangle|%d\\rangle$',o,n,m,l);
					end
					rho_SE = reshape(rho_SE,d(1),[]);
					out = {};
					out{1} = rho_SE(1:obj.lastIdx,:);				% t x D'*dk'*D*dk
					out{2} = legLab;								% D'*dk'*D*dk
				case 'is-heff-selfadjoint'
					% 'out' contains the total absolute error of Heff-Heff' as measure of self-adjointness
					d = size(obj.Heff);
					heffTemp = conj(reshape(permute(obj.Heff, [1,4,5,2,3]),d(1),[]));		% form hermitian conjugate; t x rest
					heff = reshape(obj.Heff, d(1),[]);
					out = sum(abs(heff - heffTemp),2);
					full = 1;
				case 'is-rho-se-selfadjoint'
					% 'out' contains the total absolute error of rhoSE-rhoSE' as measure of self-adjointness
					rho_SE = bsxfun(@times, permute(conj(obj.sysState), [1 3 2]) , permute(obj.sysState, [1,4,5,3,2]));  % t x D' x dk' x 1 x 1 .* t x 1 x 1 x D x dk  = t x D' x dk' x D x dk
					
					d = size(rho_SE);
					rho_SETemp = conj(reshape(permute(rho_SE, [1,4,5,2,3]),d(1),[]));		% form hermitian conjugate; t x rest
					rho_SE = reshape(rho_SE, d(1),[]);
					out = sum(abs(rho_SE - rho_SETemp),2);
					full = 1;
				case 'is-heff-current-antisymmetric'
					% 'out' contains the total absolute error of cur+cur.' as measure of anti-symmetry
					cur = obj.getData('heff-current');
					
					d = size(cur);
					curTemp = reshape(permute(cur, [1,4,5,2,3]),d(1),[]);		% form transpose; t x rest
					cur = reshape(cur, d(1),[]);
					out = sum(abs(cur + curTemp),2);
					full = 1;
				case 'is-heff-current-total-zero'
					% 'out' contains the sum over the upper right diagonal of cur
					% not working, not finished!
					cur = obj.getData('heff-current');
					
					d = size(cur);
					cur = reshape(cur,d(1),d(2)*d(3),[]);				% t x D*dk x D*dk
			end
			
			if ~full && ~iscell(out) && ~isobject(out)
				out = out(1:obj.tresults.lastIdx-idxOffset);	% truncate if wanted
			end
			
			function VC = diabaticCharacter(V)
				%% This calculates the diabatic character from trafo matrices V
				%	Performs sum in D while keeping dk
				%
				%	Input:
				%		V:  tx1 cell  of  D*dk x D*dk_eig trafos
				%
				%	Output:
				%		VC: D*dk_eig x dk x t
				
				% VC for V-contribution/ diabatic Character 
				VC = cell2mat(permute(V,[2,3,1]));												% creates D*dk x D*dk_eig x t
				VC = permute(VC, [2,1,3]);														% D*dk_eig x D*dk x t

				d = size(VC);
				VC = reshape(VC, [d(1),sqrt(d(2))*[1,1], d(3)]);								% D*dk_eig x D x dk x t
				VC = squeeze(sum(VC.*conj(VC),2));												% D*dk_eig x dk x t
				% Now: VC has percentage of dk contribution on each D*dk_eig and t
			end
		end
		
		function [h,pl,obj] = plot(obj,type,varargin)
			%% Simple 1D plots for following Observables:
			% type:
			%	'sz'
			%	'spin'
			%	'calctime'
			%	'calctime-d'
			%	'calctime-d-sec'
			%	'rhoii'
			%	'rho-current'
			%	'rho-dpmes'
			%	'rhoij-real'
			%	'rhoij-imag'
			%	'rhoij-abs'
			%	'hshi'
			%	'stateproj'
			%	'chain-x-t-avg'
			%	'linabs'
			%	'rhoii-ft'
			%	'rhoii-osc-res'
			%	'rhoii-osc-res-med'
			
			% single options:
			%	'-resetcolororder'	:	plot into existing axes given as argument. Before plotting reset the Color order to initial state.
			%
			% additional options ,Name,Value, pair:
			%	'-grid',[m n]		:	Use with object array: plots into specified grid m x n; can specify figure
			%   '-chain',m			:   for which chain to plot certain observables (mostly RC plots)
			%   '-state',m			:	for which adiabatic state to plot heff-pop
			%	'-ylim',[m n]		:	pre-define ylim for heff-pop-diab to obtain correct area sizes
			
			% initialise modifiers
			% plot(...,ax_handle)	Plots into the axes/figure given by handle
			if length(obj) > 1
				plotToGrid = 0;
				gridh = [];					% grid handle, will contain figure and axes handles
				% allow additional option grid
				if nargin-2 > 0
					m    = find(strcmp(varargin,'-grid'));
					fInd = cellfun('isclass',varargin,'matlab.ui.Figure');
					axInd = cellfun('isclass',varargin,'matlab.graphics.axis.Axes');			% override!
					
					if ~isempty(m)						% plot into grid
						if length(varargin) < m+1 || isempty(varargin{m+1})
							[gridSize(1),gridSize(2)] = TDVPData.bestGrid(length(obj));
						else
							gridSize = varargin{m+1};
						end
						plotToGrid = 1;
						if sum(fInd) == 1
							gridh = TDVPData.plotGrid(gridSize(1),gridSize(2),varargin{fInd});
						else
							gridh = TDVPData.plotGrid(gridSize(1),gridSize(2));
						end
						if any(axInd)
							varargin{axInd} = [];														% cannot plot multiple plots into single axes
						end
						if any(fInd)
							varargin{fInd}  = [];
						end
						[varargin{[m,m+1]}] = deal([]);
					end
				end
				% deal with object array
				pl = gobjects(length(obj),1);
				for ii = 1:length(obj)
					if plotToGrid
						[hTemp,plTemp,obj(ii)] = obj(ii).plot(type,varargin{:},gridh.ax(ii));
					else
						[hTemp,plTemp,obj(ii)] = obj(ii).plot(type,varargin{:});
					end
					if isempty(obj(ii).LegLabel)
						t = title(obj(ii).folder);
					else
						t = title(obj(ii).LegLabel);
					end
					if plotToGrid
						t.Units  = 'normalized';
						t.Position = [0.5,0.9];
					end
					plTemp = reshape(plTemp,1,[]);
					if size(pl,2) ~= size(plTemp,2)
						pl(end,size(plTemp,2)) = gobjects(1,1);
					end
					pl(ii,1:length(plTemp)) = plTemp;
					h(ii) = hTemp;
				end
				return;
			end
			h.xlbl = '$t$';
			h.ylbl = [];
			h.leglbl = {};
			h.xdata = [];
			h.ydata = [];
			h.ylim = [];
			h.xlim = [];
			h.pl = [];					% plot handle
			pl = [];					% old plot handle;
			h.chain = [];				% chain number needed for (a)diabatic n and x plots
			h.state = [];				% state number needed for heff plots
			h.color = {};				% cell array with color for each line, if specified
			
			ts = 1;						% scale time axis
			eScale = 0;					% energy units used for DFT: 0 = none, 1 = eV, 2 = cm (wavenumber)
			unicolor = 0;
			resetColorOrder = 0;		% start from index 1
			plotOpt = {};				% additional plotOptions from varargin
			DFTplot = 0;				% plotting DFT data
			DFTshift = 0;				% shift the fft results
			h.normalise = 0;				% normalise to maximum peak (for FFT especially)
			h.distribute = 0;			% distribute plot lines along y; only use if h.normalise = 1, otherwise gets messy
			h.rowheight  = 2;			% defaults for grid
			h.rowwidth  = 7.2;			% defaults for grid
			h.patchthickness = 0.2;		% in % of total ylim
			h.visible = 1;				% default for fig visibility
			
			h.f  = [];
			h.ax = [];
			
			for m = 1:nargin-2
				if isobject(varargin{m})
					if isa(varargin{m},'matlab.ui.Figure')
						h.f  = varargin{m};
					elseif isa(varargin{m},'matlab.graphics.axis.Axes')
						h.ax = varargin{m};									% the last axis added will be used!
					end
					continue;
				elseif isempty(varargin{m})
					continue;
				end
				switch lower(varargin{m})
					case '-fsev'
						% fs timescale for H in eV and eV scale in DFT plots
						ts = 0.658;
						eScale = 1;
						h.xlbl = '$t/\mathrm{fs}$';
					case '-cmev'
						% cm scale in DFT plots
						ts = 0.658;
						eScale = 2;
						h.xlbl = '$E/\mathrm{cm}^{-1}$';
					case '-nmev'
						% nm scale in DFT plots
						ts = 0.658;
						eScale = 3;
						h.xlbl = '$\lambda/\mathrm{nm}$';
					case '-fftshift'
						DFTshift = 1;
					case '-unicol'
						% same colors on all lines
						unicolor = 1;
					case '-resetcolororder'
						resetColorOrder = 1;
					case '-norm'
						h.normalise = 1;
					case '-dist'
						h.distribute = 1;
					case '-chain'
						h.chain = varargin{m+1};
						[varargin{[m,m+1]}] = deal([]);
					case '-state'
						h.state = varargin{m+1};
						[varargin{[m,m+1]}] = deal([]);
					case '-xlim'
						h.xlim = varargin{m+1};
						[varargin{[m,m+1]}] = deal([]);
					case '-ylim'
						h.ylim = varargin{m+1};
						[varargin{[m,m+1]}] = deal([]);
					case '-rowheight'
						h.rowheight = varargin{m+1};
						[varargin{[m,m+1]}] = deal([]);
					case '-rowwidth'
						h.rowwidth = varargin{m+1};
						[varargin{[m,m+1]}] = deal([]);
					case '-patchthickness'
						h.patchthickness = varargin{m+1};
						[varargin{[m,m+1]}] = deal([]);
					case '-invisible'
						h.visible = 0;
					otherwise
						% pass through as direct plot options!
						plotOpt = [plotOpt, varargin(m)];
				end
			end

			% Create figure handles and panels
			if isempty(h.f)
				if isempty(h.ax)
					h.f = figure();
				else
					axes(h.ax(1));					% select axes
					h.f = gcf;						% get their figure
				end
			end
			if isempty(h.ax)
				h.ax = gca;
			end
			h.f.Renderer = 'painters';				% provides much better 2D results generally!
			hold all;
			box on;
			if ~h.visible
				h.f.Visible = 'off';
			end
			
			if resetColorOrder
				h.ax.ColorOrderIndex = 1;
			end
			
			if unicolor
				idx = h.ax.ColorOrderIndex;
			end
			
			switch lower(type)						% all the plot types
				case 'sz'
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = obj.tresults.spin.sz(1:obj.lastIdx);
					h.leglbl = {obj.LegLabel};
					h.ylbl = '$\left<\sigma_z\right>$';
				case 'spin'
					if isempty(obj.spin)
						disp('please calculate spin first');
						return;
					else
						h.xdata = obj.t(1:obj.lastIdx)*ts;
						h.ydata = [obj.spin,sqrt(obj.spin(:,1).^2+obj.spin(:,2).^2)];
						h.leglbl = {'$\left<\sigma_x\right>$','$\left<\sigma_y\right>$','$\left<\sigma_z\right>$','$\sqrt{\left<\sigma_x\right>^2+\left<\sigma_y\right>^2}$'}';
					end
					h.ylbl = '$\left<\sigma\right>$';
				case 'calctime'
					idx = find(obj.para.tdvp.calcTime,1,'last');
					h.xdata  = obj.para.tdvp.t(2:idx+1)*ts;
					h.ydata  = obj.para.tdvp.calcTime(1:idx);
					h.leglbl = {obj.LegLabel};
					h.ylbl   = 'CPU time/h';
				case 'calctime-d'
					% 1st derivative, in hours
					h.xdata  = obj.para.tdvp.t(2:obj.lastIdx-1)*ts;
					h.ydata  = diff(obj.para.tdvp.calcTime(1:obj.lastIdx-1));
					h.leglbl = {obj.LegLabel};
					h.ylbl = 'CPU time/sweep/h';
				case 'calctime-d-sec'
					% 1st derivative, in seconds
					idx = find(obj.para.tdvp.calcTime,1,'last');
					h.xdata = obj.para.tdvp.t(2:idx)*ts;
					h.ydata = 3600*diff(obj.para.tdvp.calcTime(1:idx));
					h.leglbl = {obj.LegLabel};
					h.ylbl = 'CPU time/sweep/s';
				case 'rhoii'
					box on;
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = abs(obj.getData('rhoii'));
					h.leglbl = obj.sysLabel;
					if ~isempty(h.state)
						h.ax.ColorOrderIndex = h.state;		% ensure same colors
						h.ydata = h.ydata(:,h.state);
						h.leglbl = h.leglbl(h.state);
					end
					h.ylbl = '$\rho_{ii}(t)$';
				case 'rhoii-diff'
					%% Plot the Time derivative of the system's populations
					box on;
					h.xdata = obj.t(1:obj.lastIdx-1)*ts;
					h.ydata = diff(abs(obj.getData('rhoii')))/(obj.dt*ts);
					h.leglbl = cellfun(@(x) sprintf('$d\\rho_{%s}/dt$',x(2:end-1)),obj.sysLabel,'UniformOutput',false);
					if ~isempty(h.state)
						h.ax.ColorOrderIndex = h.state;		% ensure same colors
						h.ydata = h.ydata(:,h.state);
						h.leglbl = h.leglbl(h.state);
					end
					h.ylbl = '$\rho_{ii}(t)$';
				case 'rhoij-real'
					if ~isfield(obj.tresults,'rho'), return, end;
					r = obj.getData('rhoij');							% {rhoij, info}
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = real(r{1});
					h.leglbl = arrayfun(@(i) sprintf('$%d\\to%d$',r{2}(1,i),r{2}(2,i)),(1:size(r{2},2))','UniformOutput',false);
					h.ylbl = 'Re($\rho_{ij})$';
				case 'rhoij-imag'
					if ~isfield(obj.tresults,'rho'), return, end;
					r = obj.getData('rhoij');							% {rhoij, info}
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = imag(r{1});
					h.leglbl = arrayfun(@(i) sprintf('$%d\\to%d$',r{2}(1,i),r{2}(2,i)),(1:size(r{2},2))','UniformOutput',false);
					h.ylbl = 'Im($\rho_{ij})$';
				case 'rhoij-abs'
					if ~isfield(obj.tresults,'rho'), return, end;
					r = obj.getData('rhoij');							% {rhoij, info}
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = abs(r{1});
					h.leglbl = arrayfun(@(i) sprintf('$%d\\to%d$',r{2}(1,i),r{2}(2,i)),(1:size(r{2},2))','UniformOutput',false);
					h.ylbl = '$|\rho_{ij}|$';
				case 'rho-current'
					if ~isfield(obj.tresults,'rho'), return, end;
					r = obj.tresults.rho(1:obj.lastIdx,:,:);
					d = size(r);
					rhoDiag = zeros(d([1,2]));
					LegLabD = cell(1,d(2));
					rhoCur = zeros(d(1),d(2)*(d(2)-1)/2-1);				% calculates number of off-diagonals
					LegLabCur = cell(1,size(rhoCur,2));
					iC = 1;
					for jj = 1:d(2)
						for ii = 1:jj
							if ii == jj
								rhoDiag(:,ii) = real(r(:,ii,jj));		% get population
								LegLabD{ii} = sprintf('$%d$',ii);
							else
								rhoCur(:,iC)  = real(r(:,ii,jj));		% get Current
								LegLabCur{iC} = sprintf('$%d\\to%d$',ii,jj);
								iC = iC+1;
							end
						end
					end
					% resize ax and add top axis for current!
					h.ax.Position(4) = h.ax.Position(4)/2;
					box on;
 					pl1 = plot(obj.t(1:obj.lastIdx)*ts,rhoDiag, plotOpt{:},'Displayname',obj.LegLabel);
					h.ax2 = axes('Position',h.ax.Position);
					h.ax2.Position(2) = h.ax2.Position(2)+h.ax2.Position(4);
					pl2 = plot(obj.t(1:obj.lastIdx)*ts,rhoCur, plotOpt{:},'Displayname',obj.LegLabel);
% 					ax2.XTickLabel = [];
					h.ax2.XAxisLocation = 'top';
					h.ax2.YAxisLocation = 'right';
					ylabel('Im($\rho_{ij})$')
					pl = [pl1;pl2];
					set(pl,{'Displayname'},[LegLabD,LegLabCur]');
					h.ylbl = '$\rho_{ii}$';
					axes(h.ax);
				case 'rho-dpmes'
					% Plots rhoii, rhoij-real and rhoij-imag stacked in 3 axes
					% start explicitely new figure?
					LegLabij = {'LE$^- \leftrightarrow $TT','CT$^+ \leftrightarrow $LE$^+$','CT$^- \leftrightarrow $TT','CT$^- \leftrightarrow $LE$^-$'};
					
					% rhoii:
					h.ax(1).Position(4) = h.ax(1).Position(4)/3;
					hsub = obj.plot('rhoii',h.ax(1),varargin{:});
					legend('TT','LE$^+$','LE$^-$','CT$^+$','CT$^-$');
					h.ylbl = '$\rho_{ii}$';
					col = {hsub.pl.Color};
					grid on;
					axis tight;
					% rhoij-real:
					h.ax(2) = axes('Position',h.ax(1).Position); hold all;
					h.ax(2).Position(2) = h.ax(2).Position(2)+h.ax(2).Position(4);
					hsub = obj.plot('rhoij-real',h.ax(2),varargin{:});
					h.ax(2).XTickLabels = '';
					h.ax(2).YTick = h.ax(2).YTick(2:end);
					xlabel('');
					hsub.pl([1,3,5,7,8,10]).delete;
					set(hsub.pl([2,4,6,9]),{'Color'},col(1:4)');			% reset color order
					grid on;
					axis tight;
					h.ax(2).YLim = h.ax(2).YLim + [-1,1]*0.1*abs(diff(h.ax(2).YLim));
					
					% rhoij-imag:
					h.ax(3) = axes('Position',h.ax(2).Position); hold all;
					h.ax(3).Position(2) = h.ax(3).Position(2)+h.ax(3).Position(4);
					hsub = obj.plot('rhoij-imag',h.ax(3),varargin{:});
					h.ax(3).XTickLabels = '';
					h.ax(3).YTick = h.ax(3).YTick(2:end);
					xlabel('');
					hsub.pl([1,3,5,7,8,10]).delete;
					set(hsub.pl([2,4,6,9]),{'Color'},col(1:4)');			% reset color order
					leg = legend(LegLabij{:});
					leg.Position = [0.4,0.6,0.14,0.15];
					grid on;
					axis tight;
					h.ax(3).YLim = h.ax(3).YLim + [-1,1]*0.1*abs(diff(h.ax(3).YLim));
					axes(h.ax(1));						% select first axis again
				case 'rhoii-ft'
					% applies DFT to the population probability of rho
					DFTplot = 1;
					h.data = squeeze(real(obj.getData('rhoii')));		% t x L x nChain
% 					h.data = TDVPData.meanRes(h.data);
					h.data = TDVPData.movAvgRes(h.data,350/obj.dt);
					h.tdata = obj.t(1:obj.lastIdx);
					h.Fs = 1/diff(obj.t(1:2));
					h.window = 'kaiser';
					h.zeroPadFact = 10;
					h.reassign = 1;
					h = TDVPData.FFT(h);
					f = h.xdata;
					h.ylbl = '$FT(\rho_{ii})$';
					h.tlbl = 'Population Fourier';
				case 'rhoii-osc-res'
					% applies DFT to the population probability of rho
					DFTplot = 1;
					m = obj.lastIdx;					% last datapoint to include
					window = hann(m,'periodic');
% 					maxN = m;
					maxN = pow2(nextpow2(10*m));			% if > lastIdx -> zero padding
					f = (0:maxN-1)/obj.t(2)/maxN;
					if isempty(obj.rhoOscRes)
						obj.rhoOscRes = obj.getData('rhoii-osc-res');		% t x nStates
					end
					rhoii = obj.rhoOscRes;
					ft = fft(real(rhoii(1:m,:)).*(window*ones(1,size(rhoii,2))),maxN,1);
% 					ft = fft(real(rhoii(1:m,:)),maxN,1);
					if DFTshift
						ft = fftshift(ft,1);
						f = f-f(end)/2;
					end
					pl = plot(f, abs(ft));
					h.ylbl = '$FT(\rho_{ii})$';
				case 'rhoii-osc-res-med'
					% applies DFT to the population probability of rho
					DFTplot = 1;
					m = obj.lastIdx;					% last datapoint to include
					window = hann(m,'periodic');
% 					maxN = m;
					maxN = pow2(nextpow2(10*m));			% if > lastIdx -> zero padding
					f = (0:maxN-1)/obj.t(2)/maxN;
					if isempty(obj.rhoOscRes)
						obj.rhoOscRes = obj.getData('rhoii-osc-res-med');		% t x nStates
					end
					rhoii = obj.rhoOscRes;
					ft = fft(real(rhoii(1:m,:)).*(window*ones(1,size(rhoii,2))),maxN,1);
% 					ft = fft(real(rhoii(1:m,:)),maxN,1);
					if DFTshift
						ft = fftshift(ft,1);
						f = f-f(end)/2;
					end
					h.xdata = f;
					h.ydata = abs(ft);
					h.ylbl = '$FT(\rho_{ii})$';
				case 'rhoii-osc-res-mean'
					% applies DFT to the population probability of rho
					DFTplot = 1;
					m = obj.lastIdx;					% last datapoint to include
					window = hann(m,'periodic');
% 					maxN = m;
					maxN = pow2(nextpow2(10*m));			% if > lastIdx -> zero padding
					f = (0:maxN-1)/obj.t(2)/maxN;
					if isempty(obj.rhoOscRes)
						obj.rhoOscRes = obj.getData('rhoii-osc-res-mean');		% t x nStates
					end
					rhoii = obj.rhoOscRes;
					ft = fft(real(rhoii(1:m,:)).*(window*ones(1,size(rhoii,2))),maxN,1);
% 					ft = fft(real(rhoii(1:m,:)),maxN,1);
					if DFTshift
						ft = fftshift(ft,1);
						f = f-f(end)/2;
					end
					h.xdata = f;
					h.ydata = abs(ft);
					h.ylbl = '$FT(\rho_{ii})$';
				case 'rho-se-abs'
					r = obj.getData('rho-se');							% {rhoij, info}
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = abs(r{1});
					h.leglbl = r{2};
					h.ylbl = '$|\rho_{SE}|$';
				case 'vne'
					% von Neumann Entropy of the system
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = obj.getData('vne');				% t x 1
					h.leglbl = arrayfun(@(i) sprintf('$%d$',i),(1:size(h.ydata,2))','UniformOutput',false);
					h.ylbl = '$S_{vNE}$';
				case 'chain-para'
					%% Plot the Chain parameters epsilon & t
					% plot into existing axes if there are enough
					NC = obj.nChains;
					[m,n] = TDVPData.bestGrid(NC);
					if (numel(h.ax) < NC)
						htemp = TDVPData.plotGrid(m,n,h.f,'rowheight',h.rowheight,'rowwidth',h.rowwidth);
						h.ax = htemp.ax;
					end
					h.pl = gobjects(m,n,2);
					for ii = 1:NC
						[mm,nn] = ind2sub([m,n],ii);
						hold(h.ax(ii),'on');
						h.pl(mm,nn,1) = plot(h.ax(ii),obj.para.chain{ii}.epsilon,'DisplayName','$\epsilon$'); hold all;
						h.pl(mm,nn,2) = plot(h.ax(ii),obj.para.chain{ii}.t,'LineStyle','--','Color',h.pl(mm,nn,1).Color,'DisplayName','$t$');
						TDVPData.decrColorOrderIndex(h.ax(ii));
						t = title(h.ax(ii),sprintf('Chain %s',obj.chainLabel{ii}));
						t.Units = 'norm';
						if t.Position(2) > 1
							t.Position(2) = t.Position(2)-0.1;
						end
						if mm == m
							xlabel(h.ax(ii),'Site $k$');
						end
						if nn == 1
							ylabel(h.ax(ii),'');
						end
					end
					return
				case 'chain-n-rc'
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = zeros(obj.lastIdx,size(obj.occC,3));
					if isfield(obj.para,'useTreeMPS') && obj.para.useTreeMPS
						for ii = 1:size(obj.occC,3)
							pos = find(obj.occC(2,:,ii),1,'first');					% find pos of rc in chain ii
							h.ydata(:,ii) = real(obj.occC(1:obj.lastIdx,pos,ii));
						end
					else
						h.ydata = squeeze(real(obj.occC(1:obj.lastIdx,2,:)));		% t x L x nChain
					end
					h.leglbl = obj.chainLabel;
					h.ylbl = '$\langle n \rangle$';
				case 'chain-n-rc-diff'
					h.xdata = obj.t(1:obj.lastIdx-1)*ts;
					h.ydata = zeros(obj.lastIdx,size(obj.occC,3));
					if isfield(obj.para,'useTreeMPS') && obj.para.useTreeMPS
						for ii = 1:size(obj.occC,3)
							pos = find(obj.occC(2,:,ii),1,'first');					% find pos of rc in chain ii
							h.ydata(:,ii) = real(obj.occC(1:obj.lastIdx,pos,ii));
						end
					else
						h.ydata = squeeze(real(obj.occC(1:obj.lastIdx,2,:)));		% t x L x nChain
					end
					h.ydata = diff(h.ydata)/(obj.dt*ts);
					h.leglbl = obj.chainLabel;
					h.ylbl = '$\langle n \rangle$';
				case 'chain-n-d-rc'
					%% Plot the diabatic RC Chain occupation
					%
					if ~isempty(h.chain)
						%% this is specific call to subplot!
						h.xdata = obj.t(1:obj.lastIdx)*ts;
						rcCol = find(~all(squeeze(obj.occCd(1:obj.lastIdx,:,1,h.chain))==0),1);
						h.ydata = squeeze(real(obj.occCd(1:obj.lastIdx,rcCol,:,h.chain)));		% t x L x state x nChain
						if ~isempty(h.state)
							h.ax.ColorOrderIndex = h.chain;
							h.ydata = h.ydata(:,h.state);
						end
						h.leglbl = obj.sysLabel;
						h.ylbl = '$\langle n \rangle$';
						h.tlbl = 'Diabatic Occupation';
					else
						%% Plot for all chains: this is the first call to generate grid plot!
						hlp_t_L_state_chain_grid(obj.occCd,'chain-n-d-rc');
						return;				% exit here!
					end
				case 'chain-n-a-rc'
					% h.chain: 1x1; h.state: 1 x n possible
					if ~isempty(h.chain)
						%% this is specific call to subplot!
						h.xdata = obj.t(1:obj.lastIdx)*ts;
						rcCol = find(~all(squeeze(obj.occCa(1:obj.lastIdx,:,1,h.chain))==0),1);
						h.ydata = squeeze(real(obj.occCa(1:obj.lastIdx,rcCol,:,h.chain)));		% t x L x state x nChain
						h.leglbl = arrayfun(@(i) sprintf('$%d$',i),(1:size(h.ydata,2))','UniformOutput',false);
						if ~isempty(h.state)
							h.ax.ColorOrderIndex = h.chain;
							h.ydata = h.ydata(:,h.state);
							h.leglbl = arrayfun(@(i) sprintf('$%d$',i),h.state,'UniformOutput',false);
						end
						h.ylbl = '$\langle n \rangle$';
						h.tlbl = 'Adiabatic Occupation';
					else
						%% Plot for all chains: this is the first call to generate grid plot!
						hlp_t_L_state_chain_grid(obj.occCa,'chain-n-a-rc');
						return;				% exit here!
					end
				case 'chain-n-site2-ft'
					% applies DFT to the occupation of site 2
					DFTplot = 1;
					m = obj.lastIdx;					% last datapoint to include
					window = hann(m,'periodic');
					maxN = pow2(nextpow2(20*m));			% if > lastIdx -> zero padding
					f = (0:maxN-1)/obj.t(2)/maxN;
					h.data = squeeze(real(obj.occC(1:obj.lastIdx,2,:)));		% t x L x nChain
					h.data = TDVPData.movAvgRes(h.data,350/obj.dt);
					ft = fft(real(h.data).*(window*ones(1,size(h.data,2))),maxN,1);
					if DFTshift
						ft = fftshift(ft,1);
						f = f-f(end)/2;
					end
					h.xdata = f;
					h.ydata = abs(ft).^2;
					h.ylbl = '$FFT(\left< n_k \right>)$';
					h.tlbl = 'Chain Occupation Site 2 Fourier';
				case 'chain-n-rc-ft'
					% applies DFT to the occupation of site 2
					DFTplot = 1;
					if isfield(obj.para,'useTreeMPS') && obj.para.useTreeMPS
						for ii = 1:size(obj.occC,3)
							pos = find(obj.occC(2,:,ii),1,'first');					% find pos of rc in chain ii
							h.data(:,ii) = real(obj.occC(1:obj.lastIdx,pos,ii));
						end
					else
						h.data = squeeze(real(obj.occC(1:obj.lastIdx,2,:)));		% t x L x nChain
					end
% 					h.data = TDVPData.meanRes(h.data);
					h.data = TDVPData.movAvgRes(h.data,350/obj.dt);
					h.tdata = obj.t(1:obj.lastIdx);
					h.Fs = 1/diff(obj.t(1:2));
					h.window = 'kaiser';
					h.zeroPadFact = 10;
					h.reassign = 1;
					h = TDVPData.FFT(h);
					f = h.xdata;
					h.ylbl = '$FFT(\left< n_k \right>)$';
					h.tlbl = 'Chain Occupation Site 2 Fourier';
				case 'chain-n-rc-tft'
					% applies DFT to the occupation of site 2
					DFTplot = 1;
					h.data = real(squeeze(real(obj.occC(1:obj.lastIdx,2,:))));		% t x L x nChain
					h.data = TDVPData.movAvgRes(h.data,350/obj.dt);
					h.tdata = obj.t;
					h.Fs = 1/diff(obj.t(1:2));
					h.window = 'kaiser';
					h.zeroPadFact = 100;
					h.reassign = 1;
					h = TDVPData.tFFT(h);
					f = h.xdata;
					h.ylbl = '$FFT(\left< n_k \right>)$';
					h.tlbl = 'Chain Occupation Site 2 Fourier';
				case 'chain-x-rc'
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = zeros(obj.lastIdx,size(obj.xC,3));
					if isfield(obj.para,'useTreeMPS') && obj.para.useTreeMPS
						for ii = 1:size(obj.xC,3)
							pos = find(obj.xC(2,:,ii),1,'first');					% find pos of rc in chain ii
							h.ydata(:,ii) = real(obj.xC(1:obj.lastIdx,pos,ii));
						end
					else
						h.ydata = squeeze(real(obj.xC(1:obj.lastIdx,2,:)));		% t x L x nChain
					end
					h.leglbl = obj.chainLabel;
					h.ylbl = '$\langle x \rangle$';
				case 'chain-x-d-rc'
					if ~isempty(h.chain)
						%% this is specific call to subplot!
						h.xdata = obj.t(1:obj.lastIdx)*ts;
						% find first non-zero column == RC
						rcCol = find(~all(squeeze(obj.xCd(1:obj.lastIdx,:,1,h.chain))==0),1);
						h.ydata = squeeze(real(obj.xCd(1:obj.lastIdx,rcCol,:,h.chain)));		% t x L x state x nChain
						if ~isempty(h.state)
							h.ax.ColorOrderIndex = h.chain;
							h.ydata = h.ydata(:,h.state);
						end
						h.leglbl = obj.sysLabel;
						h.ylbl = '$\langle x \rangle$';
						h.tlbl = 'Diabatic Displacement';
						plot([min(h.xdata),max(h.xdata)],[0,0],'k:','DisplayName','$\langle x \rangle=0$');	% zero baseline
					else
						%% this is the first call to generate grid plot!
						hlp_t_L_state_chain_grid(obj.xCd,'chain-x-d-rc');
						return;				% exit here!
					end
				case 'chain-x-a-rc'
					if ~isempty(h.chain)
						%% this is specific call to subplot!
						h.xdata = obj.t(1:obj.lastIdx)*ts;
						rcCol = find(~all(squeeze(obj.xCa(1:obj.lastIdx,:,1,h.chain))==0),1);
						h.ydata = squeeze(real(obj.xCa(1:obj.lastIdx,rcCol,:,h.chain)));		% t x L x state x nChain
						h.leglbl = arrayfun(@(i) sprintf('$%d$',i),(1:size(h.ydata,2))','UniformOutput',false);
						if ~isempty(h.state)
							h.ax.ColorOrderIndex = h.chain;
							h.ydata = h.ydata(:,h.state);
							h.leglbl = arrayfun(@(i) sprintf('$%d$',i),h.state,'UniformOutput',false);
						end
						h.ylbl = '$\langle x \rangle$';
						h.tlbl = 'Adiabatic Displacement';
						plot([min(h.xdata),max(h.xdata)],[0,0],'k:','DisplayName','$\langle x \rangle=0$');	% zero baseline
					else
						%% this is the first call to generate grid plot!
						hlp_t_L_state_chain_grid(obj.xCa,'chain-x-a-rc');
						return;				% exit here!
					end
				case 'chain-x-site2-ft'
					% applies DFT to the occupation of site 2
					DFTplot = 1;
					m = obj.lastIdx;					% last datapoint to include
					window = hann(m,'periodic');
					maxN = pow2(nextpow2(20*m));			% if > lastIdx -> zero padding
					f = (0:maxN-1)/obj.t(2)/maxN;
					h.data = squeeze(real(obj.xC(1:obj.lastIdx,2,:)));		% t x L x nChain
					h.data = TDVPData.movAvgRes(h.data,760/obj.dt);
					ft = fft(real(h.data).*(window*ones(1,size(h.data,2))),maxN,1);
					if DFTshift
						ft = fftshift(ft,1);
						f = f-f(end)/2;
					end
					h.xdata = f;
					h.ydata = abs(ft);
					h.ylbl = '$FFT(\left< x_k \right>)$';
					h.tlbl = 'Chain Displacement Site 2 Fourier';
				case 'chain-x2-site2-ft'
					% applies DFT to the occupation of site 2
					DFTplot = 1;
					m = obj.lastIdx;					% last datapoint to include
					window = hann(m,'periodic');
					maxN = pow2(nextpow2(20*m));			% if > lastIdx -> zero padding
					f = (0:maxN-1)/obj.t(2)/maxN;
					h.data = squeeze(real(obj.xC2(1:obj.lastIdx,2,:)));		% t x L x nChain
					h.data = TDVPData.movAvgRes(h.data,760/obj.dt);
					ft = fft(real(h.data).*(window*ones(1,size(h.data,2))),maxN,1);
					if DFTshift
						ft = fftshift(ft,1);
						f = f-f(end)/2;
					end
					h.xdata = f;
					h.ydata = abs(ft);
					h.ylbl = '$FFT(\left< n_k \right>)$';
					h.tlbl = 'Chain Displacement squared Site 2 Fourier';
				case 'hshi'
					pl = plot(obj.t(1:obj.lastIdx)*ts,[obj.tresults.hshi(1:obj.lastIdx,:),sum(obj.tresults.hshi(1:obj.lastIdx,:),2)]);
					h.ylbl = '$\left< H_i \right>$';
				case 'stateproj'
					pl = plot(obj.t(1:obj.lastIdx)*ts, [abs(obj.stateProj(1:obj.lastIdx)),real(obj.stateProj(1:obj.lastIdx)),imag(obj.stateProj(1:obj.lastIdx))]);
					h.ylbl = 'Autocorrelation';
				case 'chain-x-t-avg'
					hold all;
					h.ydata = squeeze(sum(real(obj.xC(1:obj.lastIdx,:,:,:)),3));		% t x L x  nChain
					h.ydata = permute(h.ydata,[1 3 2]);				% t x state x L x chain
					pl = plot(obj.t(1:obj.lastIdx)*ts,h.ydata(:,:,2));
					h.ylbl = '$\left< x_k \right>$';					
				case 'linabs'
					% linear absorption as DFT of stateProj autocorrelation function
					DFTplot = 1;
					m = obj.lastIdx;					% last datapoint to include
% 					m = 5000;
					maxN = m;
% 					maxN = pow2(nextpow2(1*m));			% if > lastIdx -> zero padding
					f = (0:maxN-1)/obj.t(2)/maxN;
					linAbs = fft(conj(obj.stateProj(1:m)),maxN);	% do conj to get positive real part!
					if DFTshift
						linAbs = fftshift(linAbs);
						f = f-f(end)/2;
					end
% 					pl = plot(f, [real(linAbs),imag(linAbs),abs(linAbs)]);
					pl = plot(f, real(linAbs));
					h.ylbl = 'Linear Absorption';
				case 'pot-dyn'
					% Plots the Potential landscape as experienced by the wavepakets over time
					% highly experimental and only works with hand-defined hamiltonian
					if strcmp(obj.para.model,'SpinBoson')
						%%
						[sigmaX,~,sigmaZ] = spinop(obj.para.spinbase);			% gives XYZ operators with respect to specific main base
						H = @(x,n) - obj.para.hx./2.*sigmaX - obj.para.hz./2.*sigmaZ + obj.para.chain{1}.t(1).*sigmaZ./2 *x*sqrt(2) + eye(2)*obj.para.chain{1}.epsilon(1) *x^2;
					end
					%%
					h.xdata = obj.t(1:obj.lastIdx)*ts;
% 					rcCol = find(~all(squeeze(obj.occCa(1:obj.lastIdx,:,1,h.chain))==0),1);
					displ = squeeze(real(obj.xCa(1:obj.lastIdx,2,:,1)));		% t x L x state x nChain
					occ = squeeze(real(obj.occCa(1:obj.lastIdx,2,:,1)));		% t x L x state x nChain
					temp = arrayfun(@(i) eig(H(displ(i,2),occ(i,1)))',(1:obj.lastIdx)','UniformOutput',false);
					h.ydata = cell2mat(temp);
				case 'heff'
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					if ~isempty(h.state)
						h.xdata = obj.t(1:obj.lastIdx)*ts;
						h.ydata = real(cell2mat(arrayfun(@(i) sort(eig(squeeze(obj.Heff(i,h.state,:,h.state,:))))',[1:obj.lastIdx]','UniformOutput',false)));
						
						h.ylbl = sprintf('$E_{eff}(%d)$',h.state);
					else
						% this is the first call to generate grid plot!
						nn = ceil(sqrt(size(obj.Heff,2)));				% count the number of adiabatic states to derive grid
						htemp = TDVPData.plotGrid(nn,nn,h.f);
						for ii = 1:size(obj.Heff,2)
							obj.plot('heff',varargin{:},htemp.ax(ii),'-state',ii);
						end
						return;				% exit here!
					end
				case 'heff-pop'
					% heff, but with additional thickness according to population of given potential surface
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					if ~isempty(h.state)
						h.xdata = obj.t(1:obj.lastIdx)*ts;
						
						out = obj.getData('heff-pop','state',h.state);
						D = out{1};								% t x dk_eig
						pop = out{2};							% t x dk_eig
						
						if isempty(h.ylim)
							h.ylim = [min(D(:)),max(D(:))];
						end
						for ii = 1:size(pop,2)
							h.pl = TDVPData.plotVariance(h.xdata,D(:,ii),pop(:,ii),h.ylim, h.ax,'thickness',h.patchthickness);
						end
						
						%  V{3}'*squeeze(obj.Heff(3,1,:,1,:))*V{3}
						h.ylbl = sprintf('$E_{eff}(%d)$',h.state);
						return;
					else
						% this is the first call to generate grid plot!
						nPlots = size(obj.Heff,2);
						[mm,nn] = TDVPData.bestGrid(nPlots);
						htemp = TDVPData.plotGrid(mm,nn,h.f);
						for ii = 1:size(obj.Heff,2)
							obj.plot('heff-pop',varargin{:},htemp.ax(ii),'-state',ii);
						end
						return;				% exit here!
					end
				case 'heff-pop-diab'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					if ~isempty(h.state)
						h.xdata = obj.t(1:obj.lastIdx)*ts;
						
						out = obj.getData('heff-pop-diab','state',h.state);
						D = out{1};									% t x dk_eig
						pop = out{2};								% t x dk x dk_eig
						
						if isempty(h.ylim)
							h.ylim = [min(D(:)),max(D(:))];
						end
						for ii = 1:size(pop,2)
							h.pl = TDVPData.plotVariance(h.xdata,D(:,ii),sum(pop(:,:,ii),2), h.ylim, h.ax, 'subshades',pop(:,:,ii),'thickness',h.patchthickness);
% 							h.pl = TDVPData.plotVarianceLines(h.xdata,D(:,ii),sum(pop(:,:,ii),2), h.ylim, h.ax, 'subshades',pop(:,:,ii),'thickness',h.patchthickness);
						end
						
						h.t = title(sprintf('$E_{eff}(%d)$',h.state));
						h.t.Units  = 'normalized';
						h.t.Position = [0.5,0.9];
						return;
					else
						% this is the first call to generate grid plot!
						nPlots = size(obj.Heff,2);
						[mm,nn] = TDVPData.bestGrid(nPlots);
						htemp = TDVPData.plotGrid(mm,nn,h.f,'rowheight',h.rowheight,'rowwidth',h.rowwidth);
						% for good screen, use: ,'rowheight',12,'rowwidth',45
						% for two-column spanning pubFig use 'rowheight',4,'rowwidth',16.5
						for ii = 1:size(obj.Heff,2)
							obj.plot('heff-pop-diab',varargin{:},htemp.ax(ii),'-state',ii,'-ylim',h.ylim,'-patchthickness',h.patchthickness);
						end
						if ~isempty(h.ylim)
							set(htemp.f.Children,{'YLim'},{h.ylim});
						end
						if ~isempty(h.xlim)
							set(htemp.f.Children,{'XLim'},{h.xlim});
						end
						h = htemp;
						return;				% exit here!
					end
				case 'heff-pop-diab-swap'
					% Partly working!
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% Additionally finds discontinuities arising from switching adiabatic states and swaps accordingly to get smoother dynamics (experimental!)
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					% want to rearrange obj.Heff t x D x dk x D x dk, and obj.sysState: t x dk x D
					close(h.f);
					if ~ismatrix(h.state)
						error('for this, h.state has to be a matrix defining: [tIdx, n1, n2], which D states have to be swapped');
					end
					obj2 = obj.getData('heff-swap','state',h.state);
					
					h = obj2.plot('heff-pop-diab',varargin{:},'-xlim',h.xlim,'-ylim',h.ylim,'-rowwidth',h.rowwidth,'-rowheight',h.rowheight,'-patchthickness',h.patchthickness);
					h.obj = obj2;																		% return object for further processing
					return;
					
					%%
					% Can detect discontinuities in obj.sysState
					% Convert into probabilities to avoid complex numbers
% 					tempProb = obj.sysState(1:obj.lastIdx,:,:).*conj(obj.sysState(1:obj.lastIdx,:,:));		% t x dk x D
% 					
% 					thresh = 10^-1.4;
% 					derivs = diff(tempProb,1,1);						% derivative along time
% 					jumps = derivs.* ((derivs < (-thresh)) | (derivs>thresh));
% 					
% 					[m,n,o] = find(jumps);
% 					tSorted = sortrows([m,n,o],1);
% 					[I,J] = ind2sub([5,5],tSorted(:,2));
% 					tSorted = [tSorted(:,1),I,J,tSorted(:,3)];
% 					%%
% 					[V,D] = arrayfun(@(i) eig(squeeze(obj.Heff(i,h.state,:,h.state,:)),'vector'),[1:obj.lastIdx]','UniformOutput',false);
% 					D = real(cell2mat(D'))';
% 					[D,I] = sort(D,2);																% sort eigenvalues ascending, get I to sort V
% 					V = arrayfun(@(i) V{i}(:,I(i,:)),[1:obj.lastIdx]','UniformOutput',false);		% reorder eigenvectors accordingly
% 					% V{i}: dk x D_eig
% 					Vtemp = cell2mat(V);												% creates (dk*t) x D_eig
% 					Vtemp = reshape(Vtemp,[size(V{1},1),length(V),size(V{1},1)]);		% dk x t x D_eig
% 					Vtemp = permute(Vtemp, [2,3,1]);									% t x D_eig x dk
% 
% 					% obj.sysState: t x dk x D
% 					tempState = permute(obj.sysState(1:obj.lastIdx,:,h.state),[1,3,2]);				% t x 1 x dk
% 					pop = sum(bsxfun(@times,Vtemp ,tempState),3);
% 					pop = pop.*conj(pop);															% t x D_eig x 1
% 
% 					pop = bsxfun(@times,(Vtemp.*conj(Vtemp)),pop);									% t x D_eig x dk
% 					pop = permute(pop,[1,3,2]);														% t x dk x D_eig
% 
% 					if isempty(h.ylim)
% 						h.ylim = [min(D(:)),max(D(:))];
% 					end
% 
% 					% find discontinuities in D and swap!
% 					derivs = diff(D,1,1);
% 					[mIdx,nIdx]    = find(abs(derivs)>10^-1);
% 					
% 					obj2.plot('heff-pop-diab',varargin{:},'-ylim',h.ylim);
% 					return;
				case 'heff-full'
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					% get the fully matricised Heff
					out = obj.getData('heff-full');
					D = out{1};																		% t x D*dk_eig
					V = out{2};																		% t x D*dk x D*dk_eig
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = D;

					h.ylbl = sprintf('$E (eV)$');
% 					return
				case 'heff-full-diab'
					% color surfaces according to the mean diabatic character
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					col = get(0,'defaultaxescolororder');
					% get the fully matricised Heff
					out = obj.getData('heff-full-diab');
					D  = out{1};																		% t x D*dk_eig
					V  = out{2};																		% t x D*dk x D*dk_eig
					VC = out{3};																		% D*dk_eig x dk x t
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = D;
					
					h.ylbl = sprintf('$E/\\mathrm{eV}$');
					
					character = mean(VC,3);
					d = size(character);
					h.color = arrayfun(@(i) character(i,:)*col(1:d(2),:), (1:d(1))','UniformOutput',false);
				case 'heff-full-diab-nondiag'
					% color surfaces according to the mean diabatic character
					% equals 'heff-full-diab-partdiag-v2' since the character of the surfaces does not change
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					col = get(0,'defaultaxescolororder');
					% get the fully matricised Heff
					out = obj.getData('heff-full-diab-nondiag');
					D  = out{1};																		% t x D*dk_eig
					V  = out{2};																		% t x D*dk x D*dk_eig
					VC = out{3};																		% D*dk_eig x dk x t
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = D;
					
					h.ylbl = sprintf('$E/\\mathrm{eV}$');
					
					character = mean(VC,3);
					d = size(character);
					h.color = arrayfun(@(i) character(i,:)*col(1:d(2),:), (1:d(1))','UniformOutput',false);
				case 'heff-full-diab-partdiag'
					% color surfaces according to the mean diabatic character
					% equals 'heff-full-diab-partdiag-v2' since the character of the surfaces does not change
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					col = get(0,'defaultaxescolororder');
					% get the fully matricised Heff
					out = obj.getData('heff-full-diab-partdiag','state',h.state);
					D  = out{1};																		% t x D*dk_eig
					V  = out{2};																		% t x D*dk x D*dk_eig
					VC = out{3};																		% D*dk_eig x dk x t
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = D;
					
					h.ylbl = sprintf('$E/\\mathrm{eV}$');
					
					character = mean(VC,3);
					d = size(character);
					h.color = arrayfun(@(i) character(i,:)*col(1:d(2),:), (1:d(1))','UniformOutput',false);
				case 'heff-full-diab-v2'
					% color surfaces according to the diabatic character per time as dot series
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					col = get(0,'defaultaxescolororder');
					% get the fully matricised Heff
					out = obj.getData('heff-full-diab');
					D  = out{1};																		% t x D*dk_eig
					V  = out{2};																		% t x D*dk x D*dk_eig
					VC = out{3};																		% D*dk_eig x dk x t
					
					d = size(VC);
					h.xdata = bsxfun(@times, obj.t(1:obj.lastIdx)'*ts, ones(1,d(1)));					% t x D*dk_eig
					h.xdata = h.xdata(:);																% t*D*dk_eig
					h.ydata = D(:);
					VC = permute(VC,[3,1,2]);															% t x D*dk_eig x dk
					VC = reshape(VC,[],d(2));															% t*D*dk_eig x dk
					h.cdata = VC*col(1:d(2),:);															% t*D*dk_eig x col
					
					h.pl = scatter(h.xdata,h.ydata,1,h.cdata,'.');
					h.ylbl = sprintf('$E/\\mathrm{eV}$');
					pl = h.pl;
				case 'heff-full-diab-nondiag-v2'
					% color surfaces according to the diabatic character per time as dot series
					% equivalent to 'heff-full-diab-nondiag' which has lines
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					col = get(0,'defaultaxescolororder');
					% get the fully matricised Heff
					out = obj.getData('heff-full-diab-nondiag');
					D  = out{1};																		% t x D*dk_eig
					V  = out{2};																		% t x D*dk x D*dk_eig
					VC = out{3};																		% D*dk_eig x dk x t
					
					d = size(VC);
					h.xdata = bsxfun(@times, obj.t(1:obj.lastIdx)'*ts, ones(1,d(1)));					% t x D*dk_eig
					h.xdata = h.xdata(:);																% t*D*dk_eig
					h.ydata = D(:);
					VC = permute(VC,[3,1,2]);															% t x D*dk_eig x dk
					VC = reshape(VC,[],d(2));															% t*D*dk_eig x dk
					h.cdata = VC*col(1:d(2),:);															% t*D*dk_eig x col
					
					h.pl = scatter(h.xdata,h.ydata,1,h.cdata,'.');
					h.ylbl = sprintf('$E/\\mathrm{eV}$');
					pl = h.pl;
				case 'heff-full-diab-partdiag-v2'
					% color surfaces according to the diabatic character per time as dot series
					% equivalent to 'heff-full-diab-nondiag' which has lines
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					col = get(0,'defaultaxescolororder');
					% get the fully matricised Heff
					out = obj.getData('heff-full-diab-partdiag','state',h.state);
					D  = out{1};																		% t x D*dk_eig
					V  = out{2};																		% t x D*dk x D*dk_eig
					VC = out{3};																		% D*dk_eig x dk x t
					
					d = size(VC);
					h.xdata = bsxfun(@times, obj.t(1:obj.lastIdx)'*ts, ones(1,d(1)));					% t x D*dk_eig
					h.xdata = h.xdata(:);																% t*D*dk_eig
					h.ydata = D(:);
					VC = permute(VC,[3,1,2]);															% t x D*dk_eig x dk
					VC = reshape(VC,[],d(2));															% t*D*dk_eig x dk
					h.cdata = VC*col(1:d(2),:);															% t*D*dk_eig x col
					
					h.pl = scatter(h.xdata,h.ydata,1,h.cdata,'.');
					h.ylbl = sprintf('$E/\\mathrm{eV}$');
					pl = h.pl;
				case 'heff-full-pop'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% wrong surface color & wrong shading color
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					
					out = obj.getData('heff-full-pop');
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
% 					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						h.pl = TDVPData.plotVariance(h.xdata,D(:,ii),pop(:,ii), h.ylim, h.ax, 'thickness',h.patchthickness);
% 						h.pl = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'thickness',h.patchthickness);
					end

% 					h.t = title(sprintf('$E_{eff, full}$'));
% 					h.t.Units  = 'normalized';
% 					h.t.Position = [0.5,0.9];
				case 'heff-full-pop-nondiag'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					
					out = obj.getData('heff-full-pop-nondiag');
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
% 					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						h.pl = TDVPData.plotVariance(h.xdata,D(:,ii),pop(:,ii), h.ylim, h.ax, 'thickness',h.patchthickness);
% 						h.pl = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'thickness',h.patchthickness);
					end

% 					h.t = title(sprintf('$E_{eff, full}$'));
% 					h.t.Units  = 'normalized';
% 					h.t.Position = [0.5,0.9];
				case 'heff-full-pop-partdiag'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					
					out = obj.getData('heff-full-pop-partdiag','state',h.state);
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
% 					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						h.pl = TDVPData.plotVariance(h.xdata,D(:,ii),pop(:,ii), h.ylim, h.ax, 'thickness',h.patchthickness);
% 						h.pl = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'thickness',h.patchthickness);
					end

% 					h.t = title(sprintf('$E_{eff, full}$'));
% 					h.t.Units  = 'normalized';
% 					h.t.Position = [0.5,0.9];
				case 'heff-full-pop-v2'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% colour surfaces according to mixing of diabatic staets
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop');
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
% 					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
%						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax,'thickness',h.patchthickness);
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),pop(:,ii), h.ylim, h.ax,'thickness',h.patchthickness);
						delete(plTemp(1));
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					% plot diabatic coloured scatter ontop
					htemp = obj.plot('heff-full-diab-v2',plotArgs{:},h.ax);
					h.pl(end+1) = htemp.pl;
					h.pl(end).SizeData = 0.05;
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot	
				case 'heff-full-pop-nondiag-v2'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% colour surfaces according to mixing of diabatic staets
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop-nondiag');
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
% 					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
%						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax,'thickness',h.patchthickness);
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),pop(:,ii), h.ylim, h.ax,'thickness',h.patchthickness);
						delete(plTemp(1));
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					% plot diabatic coloured scatter ontop
					htemp = obj.plot('heff-full-diab-nondiag',plotArgs{:},h.ax);
					h.pl = [h.pl,htemp.pl'];
% 					h.pl(end).SizeData = 0.05;
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot	
				case 'heff-full-pop-partdiag-v2'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% colour surfaces according to mixing of diabatic staets
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop-partdiag','state',h.state);
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
% 					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
%						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax,'thickness',h.patchthickness);
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),pop(:,ii), h.ylim, h.ax,'thickness',h.patchthickness);
						delete(plTemp(1));
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					% plot diabatic coloured scatter ontop
					htemp = obj.plot('heff-full-diab-partdiag',plotArgs{:},'-state',h.state,h.ax);
					h.pl = [h.pl,htemp.pl'];
% 					h.pl(end).SizeData = 0.05;
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot	
				case 'heff-full-pop-diab'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% wrong surface color, correct shade color
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					
					out = obj.getData('heff-full-pop');
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'subshades',popDiab(:,:,ii),'thickness',h.patchthickness);
% 							h.pl = TDVPData.plotVarianceLines(h.xdata,D(:,ii),sum(pop(:,:,ii),2), h.ylim, h.ax, 'subshades',pop(:,:,ii),'thickness',h.patchthickness);
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					
% 					h.t = title(sprintf('$E_{eff, full}$',h.state));
% 					h.t.Units  = 'normalized';
% 					h.t.Position = [0.5,0.9];
					pl = h.pl;
% 					return;
				case 'heff-full-pop-diab-nondiag'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% colour surfaces according to mixing of diabatic staets
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop-nondiag');
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'subshades',popDiab(:,:,ii),'thickness',h.patchthickness);
						delete(plTemp(1));
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					% plot diabatic coloured scatter ontop
					htemp = obj.plot('heff-full-diab-nondiag',plotArgs{:},h.ax);
					h.pl = [h.pl,htemp.pl'];
% 					h.pl(end).SizeData = 0.05;
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot
				case 'heff-full-pop-diab-partdiag'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% colour surfaces according to mixing of diabatic staets
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop-partdiag','state',h.state);
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'subshades',popDiab(:,:,ii),'thickness',h.patchthickness);
						delete(plTemp(1));
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					% plot diabatic coloured scatter ontop
					htemp = obj.plot('heff-full-diab-partdiag',plotArgs{:},'-state',h.state,h.ax);
					h.pl = [h.pl,htemp.pl'];
% 					h.pl(end).SizeData = 0.05;
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot
				case 'heff-full-pop-diab-v2'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% colour surfaces according to mixing of diabatic staets
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop');
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'subshades',popDiab(:,:,ii),'thickness',h.patchthickness);
						delete(plTemp(1));
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					% plot diabatic coloured scatter ontop
					htemp = obj.plot('heff-full-diab-v2',plotArgs{:},h.ax);
					h.pl(end+1) = htemp.pl;
					h.pl(end).SizeData = 0.05;
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot
				case 'heff-full-pop-diab-nondiag-v2'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% colour surfaces according to mixing of diabatic staets
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop-nondiag');
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'subshades',popDiab(:,:,ii),'thickness',h.patchthickness);
						delete(plTemp(1));
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					% plot diabatic coloured scatter ontop
					htemp = obj.plot('heff-full-diab-nondiag-v2',plotArgs{:},h.ax);
					h.pl(end+1) = htemp.pl;
					h.pl(end).SizeData = 0.05;
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot
				case 'heff-full-pop-diab-partdiag-v2'
					% heff, but with additional thickness according to population of given potential surface
					% this plots the thickness according to population of diabatic states on each surface
					% colour surfaces according to mixing of diabatic staets
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop-partdiag','state',h.state);
					D       = out{1};								% t x D*dk_eig
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'subshades',popDiab(:,:,ii),'thickness',h.patchthickness);
						delete(plTemp(1));
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					% plot diabatic coloured scatter ontop
					htemp = obj.plot('heff-full-diab-partdiag-v2',plotArgs{:},'-state',h.state,h.ax);
					h.pl(end+1) = htemp.pl;
					h.pl(end).SizeData = 0.05;
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot
				case 'heff-current'
					% Now only for DPMES
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = obj.getData('heff-current');				% t x D' x dk' x D x dk
					h.ydata = squeeze(sum(sum(h.ydata,2),4));			% t x dk' x dk
					h.ylbl  = '$\langle j \rangle$';
					% generate selection of upper right and legend labels
					stateLab = {'TT','LE$^+$','LE$^-$','CT$^+$','CT$^-$'};
					legLab = cell(5,5); select = zeros(5);
					for ii = 1:5
						for jj = 1:5
							legLab{ii,jj} = sprintf('%s $\\to$ %s',stateLab{ii},stateLab{jj});
							if ii > jj 
								select(ii,jj) = 1;						% select only lower left triangular since current is antisymmetric
							end
						end
					end
					h.ydata = reshape(h.ydata, size(h.ydata,1),[]);		% t x D' * dk' x D * dk
					h.ydata = h.ydata(:,logical(select(:)));
					h.ydata = h.ydata/ts;								% convert from energy to quanta
					pl = plot(h.xdata,h.ydata);
					h.leglbl = legLab(logical(select(:)));
% 					set(pl,{'DisplayName'},legLab(logical(select(:))));
					h.pl = pl;
% 					return;
				case 'heff-current-nonzero'
					% Now only for DPMES
					% deletes zero currents
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = obj.getData('heff-current');				% t x D' x dk' x D x dk
					h.ydata = squeeze(sum(sum(h.ydata,2),4));			% t x dk' x dk
					h.ylbl  = '$\langle j \rangle$';
					% generate selection of upper right and legend labels
					stateLab = {'TT','LE$^+$','LE$^-$','CT$^+$','CT$^-$'};
					legLab = cell(5,5); select = zeros(5);
					for ii = 1:5
						for jj = 1:5
							legLab{ii,jj} = sprintf('%s $\\to$ %s',stateLab{ii},stateLab{jj});
							if ii > jj && sum(abs(h.ydata(:,ii,jj))) ~= 0
								select(ii,jj) = 1;
							end
						end
					end
					h.ydata = reshape(h.ydata, size(h.ydata,1),[]);		% t x D' * dk' x D * dk
					h.ydata = h.ydata(:,logical(select(:)));
					pl = plot(h.xdata,h.ydata);
					h.leglbl = legLab(logical(select(:)));
					h.pl = pl;
				case 'heff-current-cumsum'
					% Now only for DPMES
					% output units: should be quanta
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = obj.getData('heff-current');				% t x D' x dk' x D x dk
					h.ydata = squeeze(sum(sum(h.ydata,2),4));			% t x dk' x dk
					h.ylbl  = '$\int\langle j \rangle$/quanta (needs checking)';
					
					% generate selection of upper right and legend labels
					stateLab = {'TT','LE$^+$','LE$^-$','CT$^+$','CT$^-$'};
					legLab = cell(5,5); select = zeros(5);
					for ii = 1:5
						for jj = 1:5
							legLab{ii,jj} = sprintf('%s $\\to$ %s',stateLab{ii},stateLab{jj});
							if ii > jj 
								select(ii,jj) = 1;
							end
						end
					end
					h.ydata = reshape(h.ydata, size(h.ydata,1),[]);		% t x dk' * dk
					h.ydata = cumsum(h.ydata(:,logical(select(:))),1);	% integral approximation. units of time and hbar cancel each other, so no modification needed
					pl = plot(h.xdata,h.ydata);
					h.leglbl = legLab(logical(select(:)));
% 					set(pl,{'DisplayName'},);
					h.pl = pl;
% 					return;
				case 'heff-current-flowfrom-detail'
					% Now only for DPMES
					% sums over all acceptors
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = obj.getData('heff-current');				% t x D' x dk' x D x dk
					h.ydata = squeeze(sum(sum(h.ydata,5),4));			% t x D' x dk'
					h.ylbl  = '$\langle j \rangle$';
					
					% generate selection of upper right and legend labels
					stateLab = {'TT','LE$^+$','LE$^-$','CT$^+$','CT$^-$'};
					legLab = cell(5,5); select = zeros(5);
					for ii = 1:5
						for jj = 1:5
							legLab{ii,jj} = sprintf('$|%d\\rangle$,%s',ii,stateLab{jj});
						end
					end
					h.ydata = reshape(h.ydata, size(h.ydata,1),[]);		% t x D' * dk'
					pl = plot(h.xdata,h.ydata);
					h.leglbl = legLab(:);
					h.pl = pl;
				case 'heff-current-flowfrom-system'
					% Now only for DPMES
					% Shows which states loose population through coherent currents
					% sums over all acceptors and all Bond states
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = obj.getData('heff-current');				% t x D' x dk' x D x dk
					h.ydata = squeeze(sum(sum(sum(h.ydata,5),4),2));	% t x dk'
					h.ylbl  = '$\langle j \rangle$';
					
					% generate selection of upper right and legend labels
					stateLab = {'TT','LE$^+$','LE$^-$','CT$^+$','CT$^-$'};
					pl = plot(h.xdata,h.ydata);
					h.leglbl = stateLab(:);
					h.pl = pl;
				case 'pes-from-tes'
					% get the adiabatic states from the TES
					% use these to get the corresponding PES
					% Plot as dots colored according to diabatic mixture
					
					col = get(0,'defaultaxescolororder');
					% set to tes to avoid errors
					obj = obj.setHeffTo('tes');
					
					% get the diabatic populations as eigenvectors of TES
					out = obj.getData('pes-from-tes');
					D  = out{1};																		% t x D*dk_eig
					VC = out{2};																		% D*dk_eig x dk x t
					
					d = size(VC);
					h.xdata = bsxfun(@times, obj.t(1:obj.lastIdx)'*ts, ones(1,d(1)));					% t x D*dk_eig
					h.xdata = h.xdata(:);																% t*D*dk_eig
					h.ydata = D(:);
					VC = permute(VC,[3,1,2]);															% t x D*dk_eig x dk
					VC = reshape(VC,[],d(2));															% t*D*dk_eig x dk
					h.cdata = VC*col(1:d(2),:);															% t*D*dk_eig x col
					
					h.pl = scatter(h.xdata,h.ydata,40,h.cdata,'.');
					h.ylbl = sprintf('$E/\\mathrm{eV}$');
					pl = h.pl;
				case 'pes-from-tes-pop'
					% similar to 'heff-full-pop-diab-v2', but with PES in adiabatic basis of TES
					% set to tes to avoid errors
					obj = obj.setHeffTo('tes');
					
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop');
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					out = obj.getData('pes-from-tes');
					D       = out{1};								% t x D*dk_eig
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						h.pl = TDVPData.plotVariance(h.xdata,D(:,ii),pop(:,ii), h.ylim, h.ax, 'thickness',h.patchthickness);
% 						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'subshades',popDiab(:,:,ii),'thickness',h.patchthickness);
% 						delete(plTemp(1));
% 						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
% 						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot
				case 'pes-from-tes-pop-diab-v2'
					% similar to 'heff-full-pop-diab-v2', but with PES in adiabatic basis of TES
					% set to tes to avoid errors
					obj = obj.setHeffTo('tes');
					
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop');
					pop     = out{2};								% t x D*dk_eig;			Pop on each surface
					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					out = obj.getData('pes-from-tes');
					D       = out{1};								% t x D*dk_eig
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(pop,2)
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'subshades',popDiab(:,:,ii),'thickness',h.patchthickness);
						delete(plTemp(1));
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					% plot diabatic coloured scatter ontop
					htemp = obj.plot('pes-from-tes',plotArgs{:},h.ax);
					h.pl(end+1) = htemp.pl;
					h.pl(end).SizeData = 0.05;
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot
				case 'tes-from-pes'
					% get the adiabatic states from the PES
					% use these to get the corresponding TES
					% Plot as dots colored according to diabatic mixture
					
					col = get(0,'defaultaxescolororder');
					
					% set to pes to avoid errors
					obj = obj.setHeffTo('pes');
					% get the adiabatic populations as eigenvectors of PES
					out = obj.getData('tes-from-pes');
					D  = out{1};																		% t x D*dk_eig
					VC = out{2};																		% D*dk_eig x dk x t
					
					d = size(VC);
					h.xdata = bsxfun(@times, obj.t(1:obj.lastIdx)'*ts, ones(1,d(1)));					% t x D*dk_eig
					h.xdata = h.xdata(:);																% t*D*dk_eig
					h.ydata = D(:);
					VC = permute(VC,[3,1,2]);															% t x D*dk_eig x dk
					VC = reshape(VC,[],d(2));															% t*D*dk_eig x dk
					h.cdata = VC*col(1:d(2),:);															% t*D*dk_eig x col
					
					h.pl = scatter(h.xdata,h.ydata,40,h.cdata,'.');
					h.ylbl = sprintf('$E/\\mathrm{eV}$');
					pl = h.pl;
				case 'tes-from-pes-pop-diab-v2'
					% similar to 'heff-full-pop-diab-v2', but with TES in adiabatic basis of PES
					
					% set to pes to avoid errors
					obj = obj.setHeffTo('pes');
					
					if isempty(obj.Heff)
						error('Not available, need to extract Observable heff');
						close(h.f);
						return;
					end
					
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ylbl = '$E/\mathrm{eV}$';
					out     = obj.getData('heff-full-pop');
					popDiab = out{3};								% t x dk x D*dk_eig		Diab pop on each surface
					
					out = obj.getData('tes-from-pes');
					D       = out{1};								% t x D*dk_eig
					
					if isempty(h.ylim)
						h.ylim = [min(D(:)),max(D(:))];
					end
					for ii = 1:size(popDiab,3)
						plTemp = TDVPData.plotVariance(h.xdata,D(:,ii),squeeze(sum(popDiab(:,:,ii),2)), h.ylim, h.ax, 'subshades',popDiab(:,:,ii),'thickness',h.patchthickness);
						delete(plTemp(1));
						sel = arrayfun(@(x) ~isa(x,'matlab.graphics.GraphicsPlaceholder') && isvalid(x) ,plTemp);
						h.pl = [h.pl,plTemp(sel)];
					end
					plotArgs = {};
					if eScale
						plotArgs = [plotArgs,{'-fsev'}];
					end
					
					% plot diabatic coloured scatter ontop
					htemp = obj.plot('tes-from-pes',plotArgs{:},h.ax);
					h.pl(end+1) = htemp.pl;
					h.pl(end).SizeData = 0.05;
					pl = h.pl;
					h.ydata = [];			% delete to finish without replot
				case 'heff-full-main-bonds'
					out = obj.getData('heff-full-main-bonds');
					h.xdata = obj.t(1:obj.lastIdx)*ts;
					h.ydata = out{1};								% t x D*dk
					h.ydata = h.ydata.*conj(h.ydata);				% convert to probabilities
					h.ylbl = '$P(Sys-Bond state)$';
					h.leglbl = out{2};
					h.tlbl = 'Largest Populated Surface - System-Bond State breakdown';
				case 'ttm-norm'
					% plot norm of transfer tensors
					h.xdata = obj.t(2:obj.lastIdx)*ts;
					h.ydata = obj.tresults.TTM.Tnorm(1:obj.lastIdx-1)./(obj.dt^2);
					h.ylbl = '$|T|/\Delta t^2$';
					h.leglbl = {obj.LegLabel};
				case 'ttm-norm-hull'
					% plot norm of transfer tensors
					h.xdata = obj.t(2:obj.lastIdx)*ts;
					h.ydata = obj.tresults.TTM.Tnorm(1:obj.lastIdx-1)./(obj.dt^2);
					h.ydata = TDVPData.maxFilter(h.ydata,[1,121]);
					h.ylbl = '$|T|/\Delta t^2$';
					h.leglbl = {obj.LegLabel};
				case 'spin-ttm'
					% extrapolate spin observables using transfer tensors
					spin = obj.getData('spin-ttm','factor',10);
					h.xdata = spin.t;
					h.ydata = [spin.sx,spin.sy,spin.sz];
				otherwise
					error('TDVPData:plot','PlotType not avaliable');
			end
			
			if ~isempty(h.ydata) && isempty(h.pl)
				if h.normalise
					h.ydata = bsxfun(@rdivide,h.ydata,max(h.ydata,[],1));
				end
				h.pl = plot(h.ax,h.xdata, h.ydata, plotOpt{:});							% TODO: make fail-safe for axes arrays!
				pl = h.pl;
			end
			if ~isempty(pl) && isempty(h.pl)
				h.pl = pl;
			end
			
			if h.distribute
				if  h.normalise
					% distribute across one axis
					for kk = 1:length(h.pl)
						set(h.pl(kk),'ydata', h.pl(kk).YData+length(h.pl)-kk);
					end
				else
					% create axis for each plot
					N = length(h.pl);
					heightTot = h.ax.Position(4);
					heightSingle = heightTot/N;
					h.ax.Position(4) = heightSingle;
					tempDel = h.pl(2:end);				% temphandle for deletion
					box on; axis tight;
					% from bottom to top:
% 					for i = 2:N							% from bottom to top
% 						h.ax(i) = axes; box on; axis tight;
% 						h.ax(i).Position([1,3,4,2]) = [h.ax(i-1).Position([1,3,4]),sum(h.ax(i-1).Position([2,4]))];
% 						h.ax(i).XTickLabel = '';
% 						h.pl(i) = copyobj(h.pl(i),h.ax(i));
% 				end
% 					axes(h.ax(1));
					% from top to bottom
					h.ax.Position(2) = h.ax.Position(2)+h.ax.Position(4)*(N-1);			% shift first plot to the top
					for i = 2:N
						h.ax(i) = axes; box on; axis tight;
						h.ax(i).Position([1,3,4,2]) = [h.ax(i-1).Position([1,3,4]),-diff(h.ax(i-1).Position([2,4]))];
						h.ax(i).XTickLabel = '';
						h.pl(i) = copyobj(h.pl(i),h.ax(i));
					end
					h.ax(N).XTickLabelMode = 'auto';
					h.ax(1).XTickLabel = '';
					tempDel.delete();
				end
			end
			
			if numel(h.pl) == numel(h.leglbl)
				set(h.pl(:),{'Displayname'},h.leglbl(:));
			end
			
			if DFTplot == 1
				if ts ~= 1
					f = f/0.658;
				end
				if eScale == 1
					f = f*4.135; % in eV
					h.xlbl = '$E/\mathrm{eV}$';
				elseif eScale == 2
					f = f*4.135*8065.73; % in cm
				elseif eScale == 3
					f = 1239.84193./(f*4.135);		% in nm: hc/(E in eV);
				end
				for k = 1:length(h.pl)
					h.pl(k).XData = f(1:length(h.pl(k).XData));
				end
			end
			xlabel(h.xlbl);
			ylabel(h.ylbl);
			if ~isempty(h.xlim)
				xlim(h.xlim);
			end
			if ~isempty(h.ylim)
				ylim(h.ylim);
			end
			
			
			if ~isempty(h.color)
				set(pl,{'Color'},h.color);
			end
			if unicolor
				arrayfun(@(x) set(x,'Color',pl(1).Color), pl);
				h.ax.ColorOrderIndex = idx + 1;
			end
			
			function hlp_t_L_state_chain_grid(data,plottype)
				%% hlp_t_L_state_chain_grid(data,plottype)
				% creates grid plots for PlotTypes:
				%	chain-n-d-rc
				%	chain-n-a-rc
				%	chain-x-d-rc
				%	chain-x-a-rc
				%
				%	data: t x L x state x nChain
				nPlots = size(data,4);
				nStates = size(data,3);
				[m,n] = TDVPData.bestGrid(nPlots);
				htemp = TDVPData.plotGrid(m,n,h.f,'rowheight',h.rowheight,'rowwidth',h.rowwidth);
				h.ax = htemp.ax;
				if ~isempty(h.xlim)
					varargin = [varargin,{'-xlim',h.xlim}];
				end
				if ~isempty(h.ylim)
					varargin = [varargin,{'-ylim',h.ylim}];
				end
				if ~isempty(h.state)
					varargin = [varargin,{'-state',h.state}];
					nStates = length(h.state);
				end
				h.pl = gobjects(m,n,nStates);
				for ii = 1:numel(h.ax)
					if ii > nPlots
						delete(h.ax(ii));
						h.ax(ii) = [];
						continue;
					end
					[mm,nn] = ind2sub([m,n],ii);
					htemp = obj.plot(plottype,varargin{:},h.ax(ii),'-chain',ii);
					h.pl(mm,nn,:) = htemp.pl;
					t = title(h.ax(ii),sprintf('Chain %s',obj.chainLabel{ii}));
					t.Units = 'norm';
					if t.Position(2) > 1
						t.Position(2) = t.Position(2)-0.1;
					end
					if ii == nPlots && mod(ii,mm) ~= 0
						h.ax(ii).XTickLabelMode = 'auto';
					end	
				end
			end
		end
			
		function h   = plotSld1D(obj, type, varargin)
			%% Plot with Slider to flip through different 2D plots
			% default: Plot versus time
			h = struct();				% struct containing all infos & handles
			h.ts = 1;					% scale time axis
			h.logY = 0;					% linear z-axis by default
			h.evTocm = 0;				% convert all eV to cm^-1
			h.xdata = (obj.t(1:obj.lastIdx)*h.ts).';
			h.ydata = [];				% T x L x NC x ...  if not, needs to be reshaped!
			h.noSldDims = 1;			% default: only plot 1 dim
			h.xSize = length(h.xdata)*ones(1,obj.nChains);
			
			h.xlbl = '$t$';
			h.ylbl = '';
			h.tlbl = '';
			
			h.f  = [];
			h.ax = [];
			
			for m = 1:nargin-2
				if isobject(varargin{m})
					if isa(varargin{m},'matlab.ui.Figure')
						h.f  = varargin{m};
					elseif isa(varargin{m},'matlab.graphics.axis.Axes')
						h.ax = varargin{m};
					end
					continue;
				end
					
				switch lower(varargin{m})
					case '-fsev'
						% fs timescale for H in eV
						h.ts = 0.658;
						h.xlbl = '$t/\mathrm{fs}$';
						h.xdata = h.xdata*h.ts;
					case '-log'
						% Logarithmic scale in z
						h.logY = 1;
					case '-cm'
						% use cm^-1 as Units
						h.evTocm = 1;
				end
			end
			
			% Create figure handles and panels
			if isempty(h.f)
				if isempty(h.ax)
					h.f = figure();
				else
					axes(h.ax);					% select axes
					h.f = gcf;					% get their figure
				end
			end
			if isempty(h.ax)
				h.ax = gca;
			end
			% fix axes size and add space to the right
			h.ax.Units = 'pixels';
			h.f.Position(3) = h.f.Position(3)+100;
			h.controlPanel = uipanel('Tag','Panel','Units','pixels','Position',[h.f.Position(3)-100,1,100,h.f.Position(4)]);
				
			% set normalised units to enable scaling with window size
			h.ax.Units = 'norm';
			h.controlPanel.Units = 'norm';
			
% 			if ~isempty(strfind(lower(type),'star'))
% 				% need to find lengths of chains in star picture -> length of x-axis
% 				if strcmp(obj.para.chain{1}.mapping, 'LanczosTriDiag')		% take a representative
% % 					h.xSize = cellfun(@(x) length(x.xi), obj.para.chain);
% 				else
% % 					h.xSize = size(obj.omega,1)*ones(1,h.nChains);
% 				end
% 			end
			
			switch lower(type)
				case 'chain-n-t'
					% Plot n(t), slider in k & NC
					h.ydata = real(obj.occC(1:obj.lastIdx,:,:));		% t x L x nChain
					h.ylbl = '$\left< n_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Occupation vs t';
				case 'chain-n-d-t'
					% Plot n(t), slider in k & NC
					h.ydata = real(obj.occCd(1:obj.lastIdx,:,:,:));		% t x L x state x nChain
					h.ydata = permute(h.ydata,[1 3 2 4]);				% t x state x L x chain
					h.noSldDims = 2;
					h.ylbl = '$\left< n_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Occupation vs t - diabatic';
				case 'chain-n-a-t'
					% Plot n(t), slider in k & NC
					h.ydata = real(obj.occCa(1:obj.lastIdx,:,:,:));		% t x L x state x nChain
					h.ydata = permute(h.ydata,[1 3 2 4]);				% t x state x L x chain
					h.noSldDims = 2;
					h.ylbl = '$\left< n_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Occupation vs t - adiabatic';
				case 'chain-n-c-t'
					% Plot n(t), slider in k & NC
					h.ydata = real(obj.occCc(1:obj.lastIdx,:,:,:));		% t x L x state x nChain
					h.ydata = permute(h.ydata,[1 3 2 4]);				% t x state x L x chain
					h.noSldDims = 2;
					h.ylbl = '$\left< n_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Occupation vs t - coherence';
				case 'chain-n-d-t-sum'
					%% Plot n(t), slider in k & NC
					% sum over L
					h.ydata = real(obj.occCd(1:obj.lastIdx,:,:,:));		% t x L x state x nChain
					h.ydata = sum(h.ydata,2);
					h.ydata = permute(h.ydata,[1 3 4 2]);				% t x state x chain x 1
					h.noSldDims = 2;
					h.ylbl = '$\sum_k \left< n_k \right>$';
					h.sldlbl = {'Chain'};
					h.tlbl = 'Chain Occupation vs t - diabatic';
				case 'chain-n-k'
					% Plot n(k), slider in t & NC
					h.ydata = real(obj.occC(1:obj.lastIdx,:,:));		% t x L x nChain
					h.ydata = permute(h.ydata,[2,1,3]);					% L x t x nChain
					h.ylbl = '$\left< n_k \right>$';
					h.xdata = (1:obj.para.L).';
					h.xSize = obj.para.L*ones(1,obj.nChains);
					h.xlbl = 'Site $k$';
					h.sldlbl = {'Time t =','Chain'};
					h.tlbl = 'Chain Occupation vs k';
				case 'chain-n-d-k'
					%% Plot n(k), slider in t & NC
					h.ydata = real(obj.occCd(1:obj.lastIdx,:,:,:));		% t x L x state x nChain
					h.ydata = permute(h.ydata,[2 3 1 4]);				% L x state x t x chain
					h.ylbl = '$\left< n_k \right>$';
					h.xdata = (1:obj.para.L).';
					h.xSize = obj.para.L*ones(1,obj.nChains);
					h.xlbl = 'Site $k$';
					h.sldlbl = {'State','Time t =','Chain'};
					h.tlbl = 'Chain Occupation vs k - diabatic';
				case 'chain-x-t'
					h.ydata = real(obj.xC(1:obj.lastIdx,:,:,:));		% t x L  x nChain
					h.ylbl = '$\left< x_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Displacement';
				case 'chain-x-d-t'
					h.ydata = real(obj.xCd(1:obj.lastIdx,:,:,:));		% t x L x state x nChain
					h.ydata = permute(h.ydata,[1 3 2 4]);				% t x state x L x chain
					h.noSldDims = 2;
					h.ylbl = '$\left< x_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Displacement - diabatic';
				case 'chain-x-a-t'
					h.ydata = real(obj.xCa(1:obj.lastIdx,:,:,:));		% t x L x state x nChain
					h.ydata = permute(h.ydata,[1 3 2 4]);				% t x state x L x chain
					h.noSldDims = 2;
					h.ylbl = '$\left< x_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Displacement - adiabatic';
				case 'chain-x-c-t'
					h.ydata = abs(obj.xCc(1:obj.lastIdx,:,:,:));		% t x L x nChain
					h.ydata = permute(h.ydata,[1 3 2 4]);				% t x chain x L
					h.noSldDims = 2;
					h.ylbl = '$\left< x_k \right>$';
					h.sldlbl = {'Site k ='};
					h.tlbl = 'Chain Displacement - coherence';
				case 'chain-x2-t'
					h.ydata = real(obj.xC2(1:obj.lastIdx,:,:,:));		% t x L  x nChain
					h.ylbl = '$\left< x^2_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Displacement squared';
				case 'chain-x2-d-t'
					h.ydata = real(obj.xC2d(1:obj.lastIdx,:,:,:));		% t x L x state x nChain
					h.ydata = permute(h.ydata,[1 3 2 4]);				% t x state x L x chain
					h.noSldDims = 2;
					h.ylbl = '$\left< x^2_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Displacement squared - diabatic';
				case 'chain-x2-a-t'
					h.ydata = real(obj.xC2a(1:obj.lastIdx,:,:,:));		% t x L x state x nChain
					h.ydata = permute(h.ydata,[1 3 2 4]);				% t x state x L x chain
					h.noSldDims = 2;
					h.ylbl = '$\left< x^2_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Displacement squared - adiabatic';
				case 'chain-dx-t'
					h.ydata = real(obj.xC2(1:obj.lastIdx,:,:,:)-obj.xC(1:obj.lastIdx,:,:,:).^2);		% t x L  x nChain
					h.ylbl = '$\left< x^2_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Displacement Uncertainty';
				case 'chain-dx-d-t'
					h.ydata = real(obj.xC2d(1:obj.lastIdx,:,:,:)-obj.xCd(1:obj.lastIdx,:,:,:).^2);		% t x L x state x nChain
					h.ydata = permute(h.ydata,[1 3 2 4]);				% t x state x L x chain
					h.noSldDims = 2;
					h.ylbl = '$\left< x^2_k \right>$';
					h.sldlbl = {'Site k =','Chain'};
					h.tlbl = 'Chain Displacement Uncertainty - diabatic';
				case 'chain-x-t-avg'
					hold all;
					h.ydata = squeeze(sum(real(obj.xC(1:obj.lastIdx,:,:,:)),3));		% t x L x  nChain
					h.ydata = permute(h.ydata,[1 3 2]);				% t x state x L x chain
					h.noSldDims = 2;
					h.ylbl = '$\left< x_k \right>$';
					h.sldlbl = {'Site k ='};
					h.tlbl = 'Chain Displacement';
				case 'star-n-t'
					h.ydata = real(obj.occS(1:obj.lastIdx,:,:));		% t x w x nChain
					h.yPlots = 1;
					h.ylbl = '$\left< n(\omega) \right>$';
					if ~h.evTocm
% 						h.xdata = obj.omega;							% w x nChain
% 						h.xlbl = 'Mode $\omega$ in $eV$';
					else
% 						h.xdata = obj.omega./1.23984e-4;				% w x nChain
% 						h.xlbl = 'Mode $\omega$ in $cm^{-1}$';
					end
					h.sldlbl = {'w =','Chain'};
					h.tlbl = 'Star Occupation';
				case 'star-x-t'
					hold all;
					h.ydata = real(obj.xS(1:obj.lastIdx,:,:,:));		% t x w x state x chain
					h.ydata = permute(h.ydata,[1 3 2 4]);				% t x state x w x chain
					h.noSldDims = 2;									% first 2 dims without Slider
					h.ylbl = '$\left< f_k \right>$';
					if ~h.evTocm
% 						h.xdata = obj.omega;							% w x nChain
% 						h.xlbl = 'Mode $\omega$ in $eV$';
					else
% 						h.xdata = obj.omega./1.23984e-4;				% w x nChain
% 						h.xlbl = 'Mode $\omega$ in $cm^{-1}$';
					end
					h.sldlbl = {'w =','Chain'};
					h.tlbl = 'Star Displacement';
				case 'state'
					hold all;
					h.ydata = obj.sysState(1:obj.lastIdx,:,:).*conj(obj.sysState(1:obj.lastIdx,:,:));			% t x dk x D; calc the probability
					h.ydata = permute(h.ydata,[1,3,2]);
					h.ydata(:,end+1,:) = sum(h.ydata,2);
					h.ylbl = '$|\Psi_k|^2$';
					h.noSldDims = 2;
					h.tlbl = 'Adiabatic State evolution';
					h.sldlbl = {'State #'};
				case 'state-adiab'
					hold all;
					h.ydata = obj.sysState(1:obj.lastIdx,:,:).*conj(obj.sysState(1:obj.lastIdx,:,:));			% t x dk x D; calc the probability
					h.ydata(:,end+1,:) = sum(h.ydata,2);
					h.ylbl = '$|\Psi_k|^2$';
					h.noSldDims = 2;
					h.tlbl = 'Adiabatic State evolution';
					h.sldlbl = {'Bond #'};
				otherwise
					return;
			end
			
			h.f.Name = h.tlbl;
			
			if h.logY
				h.ydata = sign(h.ydata) .* log10(abs(h.ydata));
				h.ylbl = ['$\log_{10}',h.ylbl(2:end)];
			end
			
			h.nSld = ndims(h.ydata)-h.noSldDims;					% number of sliders needed
			h.SldIdx = num2cell(ones(1,h.nSld),1);					% indices for each slider dimension, indicating currently displayed slice
			h.ysize = size(h.ydata);
			
			if h.noSldDims == 1
				h.pl{1} = plot(h.xdata(1:h.xSize(1),1),h.ydata(:,h.SldIdx{:}));
			elseif h.noSldDims == 2
				hold all;
				for ii = 1:size(h.ydata,2)
					h.pl{ii} = plot(h.xdata(1:h.xSize(1),1),h.ydata(:,ii,h.SldIdx{:}));
				end
			end
			axis tight
			grid on
			xlabel(h.xlbl);
			ylabel(h.ylbl);
			
			% Create Slider(s) inside h.controlPanel
			h.sld_w = 70; h.sld_h = 20;
			h.cP_padIn = 5;		% controlPanel padding inside
			for ii = 1:h.nSld
				h.sldText{ii} = uicontrol('Parent',h.controlPanel,'Style','text','String',sprintf('%s %d',h.sldlbl{ii},1),...
										  'HorizontalAlignment','left','FontSize',10,'FontName',h.ax.FontName);
				h.sld{ii} = javax.swing.JScrollBar(0,1,1,1,h.ysize(h.noSldDims+ii)+1);		%JScrollBar(int orientation, int value, int extent, int min, int max)
				[~,h.sldContainer{ii}] = javacomponent(h.sld{ii}, [5, 5,h.sld_w,h.sld_h], h.controlPanel);	% position defined in posDisplay()
				h.sld{ii}.setUnitIncrement(1); h.sld{ii}.setBlockIncrement(3);
				h.hsld{ii} = handle(h.sld{ii},'CallbackProperties');
				set(h.hsld{ii},'AdjustmentValueChangedCallback',@(source,callbackdata) callback_1D_Sld(source, callbackdata,ii));%(source,callbackdata)  set(h.ax,'UserData',round(source.Value)));
			end
			posDisplay();
			set(h.f, 'ResizeFcn',@posDisplay);
			
			function callback_1D_Sld(source, callbackdata, n)
				% Callback for 2D Slider plots
				% needs to be sub-function due to shared handle h.
				%  n: number of dimension in z being changed
				h.SldIdx{n} = round(source.Value);
				h.sldText{n}.String = sprintf('%s %d', h.sldlbl{n},round(source.Value));
				for kk = 1:length(h.pl)
					if length(h.pl) == 1
						set(h.pl{1} ,'ydata', h.ydata(:,h.SldIdx{:}));
					else
						set(h.pl{kk},'ydata', h.ydata(:,kk,h.SldIdx{:}));
					end
				end
			end
			
			function posDisplay(varargin)
				for kk = 1:h.nSld
					h.sldText{kk}.Position      = [h.cP_padIn, h.f.Position(4)-(h.cP_padIn*kk+2*(kk-0.5)*h.sld_h) ,h.sld_w,h.sld_h];
					h.sldContainer{kk}.Position = [h.cP_padIn, h.sldText{kk}.Position(2)-h.sld_h,h.sld_w,h.sld_h];
				end
			end
			
		end
		
		function [h,pl]   = plotSld1DFT(obj, type, varargin)
			%% Plot with Slider to change DFT range and window / padding
			% No slider in dataset here. Only in DFT parameters
			% default: Plot versus freq
			if length(obj) > 1
				% deal with object array
				pl = gobjects(length(obj),1);
				for ii = 1:length(obj)
					[hTemp,temp] = obj(ii).plotSld1DFT(type,varargin{:});
					temp = reshape(temp,1,[]);
					if size(pl,2) ~= size(temp,2)
						pl(end,size(temp,2)) = gobjects(1,1);
					end
					pl(ii,:) = temp;
					h(ii) = hTemp;
				end
				return;
			end	
			h = struct();				% struct containing all infos & handles
			h.freqScale = 1;			% scale freq axis
			h.fInvert = 0;				% need 1/f for nm scale?
			h.logY = 0;					% linear y-axis by default
			h.evTocm = 0;				% convert all eV to cm^-1
			h.normalise = 0;
			h.distribute = 0;			% distribute plot lines along y; only use if h.normalise = 1, otherwise gets messy
			h.smoothRes = 0;			% apply to residual after smoothing?
			h.meanRes = 0;				% apply to residual after mean
			h.xdata = obj.t(1:obj.lastIdx).';	% the time axis
			h.ydata = [];				% T x L x NC x ...  if not, needs to be reshaped!
			h.noSldDims = 1;			% default: only plot 1 dim
			h.xSize = length(h.xdata)*ones(1,obj.nChains);
			
			
			h.xlbl = '$f$';
			h.ylbl = '';
			h.tlbl = '';
			h.llbl = {};
			
			h.f  = [];
			h.ax = [];
			
			% FFT settings
			h.data = [];								% FFT done in 1st dimension
			h.dataRange = obj.lastIdx;					% last (real) datapoint to include
			h.zeroPadFactor = 4;						% multiples of dataset to pad
			h.power2 = 1;								% FFT in power of 2?
			h.useWindowFcn = 1;							% whether to use the (1) Hann or (2) Kaiser window function
			h.FFTshift = 0;
			
			for m = 1:nargin-2
				if isobject(varargin{m})
					if isa(varargin{m},'matlab.ui.Figure')
						h.f  = varargin{m};
					elseif isa(varargin{m},'matlab.graphics.axis.Axes')
						h.ax = varargin{m};
					end
					continue;
				end
				if isnumeric(varargin{m})
					continue;
				end
				switch lower(varargin{m})
					case '-ev'
						% ev scale for H in eV
						h.freqScale = h.freqScale/0.658*4.135;
						h.xlbl = '$E/\mathrm{eV}$';
					case '-cmev'
						% use cm^-1 as Units
						h.freqScale = h.freqScale/0.658*4.135*8065.73;
						h.xlbl = '$E/\mathrm{cm}^{-1}$';
					case '-nmev'
						% nm scale, H in ev
						h.freqScale = h.freqScale/0.658*4.135/1239.84193;
						h.fInvert = 1;
						h.xlbl = '$\lambda/\mathrm{nm}$';
					case '-log'
						% Logarithmic scale in z
						h.logY = 1;
					case '-fftshift'
						h.FFTshift = 1;
					case '-norm'
						h.normalise = 1;
					case '-dist'
						h.distribute = 1;
					case '-smoothres'
						if nargin-2 > m && isnumeric(varargin{m+1})
							h.smoothRes = varargin{m+1};
						else
							h.smoothRes = 1;
						end
					case '-meanres'
						h.meanRes = 1;
					case '-kaiser'
						h.useWindowFcn = 2;
				end
			end
			
			% Create figure handles and panels
			if isempty(h.f)
				if isempty(h.ax)
					h.f = figure();
				else
					axes(h.ax);					% select axes
					h.f = gcf;					% get their figure
				end
			end
			if isempty(h.ax)
				h.ax = gca;
			end
			% fix axes size and add space to the right
			h.ax.Units = 'pixels';
			h.cP_width = 150;						% control Panel width
			h.f.Position(3) = h.f.Position(3)+h.cP_width;
			h.controlPanel = uipanel('Tag','Panel','Units','pixels','Position',[h.f.Position(3)-h.cP_width,1,h.cP_width,h.f.Position(4)]);
			
			% set normalised units to enable scaling with window size
			h.ax.Units = 'norm';
			h.controlPanel.Units = 'norm';
			
			% define dataset for FFT
			switch lower(type)
				case 'linabs'
					% linear absorption as DFT of stateProj autocorrelation function
					h.data = conj(obj.stateProj(1:obj.lastIdx));
					h.useWindowFcn = 0;								% needs to be switched off!!
				case 'rhoii'
					% applies DFT to the population probability of rho
					rhoii = gettRhoiiSystem(obj);		% t x nStates
					h.data = real(rhoii(1:obj.lastIdx,:));
					h.noSldDims = 2;					% plot 2nd dim simultaneously
				case 'rhoii-osc-res'
					rhoii = obj.getData('rhoii-osc-res');
					h.data = rhoii;
					h.noSldDims = 2;
% 					h.useWindowFcn = 0;
				case 'rhoij-imag'
					r = obj.getData('rhoij');							% {rhoij, info}
					h.llbl = arrayfun(@(i) sprintf('$%d\\to%d$',r{2}(1,i),r{2}(2,i)),(1:size(r{2},2))','UniformOutput',false);
					h.data = imag(r{1});
					h.noSldDims = 2;
% 					h.useWindowFcn = 0;
				case 'chain-n'
					h.data = squeeze(real(obj.occC(1:obj.lastIdx,2,:)));		% t x L x nChain
					h.ylbl = '$FFT(\left< n_k \right>)$';
					h.tlbl = 'Chain Occupation Site 2 Fourier';
					h.llbl = arrayfun(@(i) sprintf('Chain $%d$',i),(1:size(h.data,2))','UniformOutput',false);
					h.noSldDims = 2;
% 					h.useWindowFcn = 0;
				case 'chain-x2'
					h.data = squeeze(real(obj.xC2(1:obj.lastIdx,2,:)));		% t x L x nChain
					h.ylbl = '$\left< n_k \right>$';
					h.tlbl = 'Chain displacement squared Site 2 Fourier';
					h.llbl = arrayfun(@(i) sprintf('Chain $%d$',i),(1:size(h.data,2))','UniformOutput',false);
					h.noSldDims = 2;
% 					h.useWindowFcn = 0;
                case 'chain-x-c'
                    h.data = abs(obj.xCc(1:obj.lastIdx,2,:,:));		% t x L x nChain
					h.data = permute(h.data,[1 3 2 4]);				% t x chain x L
					h.noSldDims = 2;
					h.ylbl = '$\left< x_k \right>$';
					h.sldlbl = {'Site k ='};
                    h.llbl = arrayfun(@(i) sprintf('Chain $%d$',i),(1:size(h.data,2))','UniformOutput',false);
					h.tlbl = 'Chain Displacement - coherence';
				otherwise
					warning('TDVPData:plotSld1DFT','Plot type not existent!');
					return;
			end
			
			if h.smoothRes == 1
				h.data = TDVPData.movAvgRes(h.data,760/obj.dt);					% 170: 300cm or 350: 145cm or 760: 66cm accuracy
			elseif h.smoothRes ~= 0
				h.data = TDVPData.movAvgRes(h.data,h.smoothRes/obj.dt);					% 170: 300cm or 350: 145cm or 760: 66cm accuracy
			end
			
			if h.meanRes == 1
				h.data = TDVPData.meanRes(h.data);
			end
			
			calcFFT();				% fills h.ydata with FFT
			
			h.f.Name = h.tlbl;
			
			if h.logY
				h.ydata = sign(h.ydata) .* log10(abs(h.ydata));
				h.ylbl = ['$\log_{10}',h.ylbl(2:end)];
			end
			
			if ~isvector(h.ydata)
				h.nSld = 2+ndims(h.ydata)-h.noSldDims;				% number of sliders needed; +1 for DFT
			else
				h.nSld = 2;
			end
			h.ysize = size(h.data);
			h.sldlbl = {'tMax','nPad'};
			h.sldInitVal= [h.ysize(1),0];
			h.sldLimits = [1,h.ysize(1); 0, 100];					% [min1, max1; min2, max2; ...]
			h.SldIdx = num2cell(ones(1,h.nSld),1);					% indices for each slider dimension, indicating currently displayed slice
			
			h.pl = plot(h.xdata, plotSpec());		% plotspec returns the actual data from h.ydata type-specific
			if ~isempty(h.llbl)
				set(h.pl,{'DisplayName'},h.llbl);					% h.llbl needs to be same shape as h.pl
			end
			
			if h.distribute
				if  ~h.normalise
					% create axis for each plot
					N = length(h.pl);
					heightTot = h.ax.Position(4);
					heightSingle = heightTot/N;
					h.ax.Position(4) = heightSingle;
					tempDel = h.pl(2:end);				% temphandle for deletion
% 					box on; axis tight;
					% from bottom to top:
% 					for i = 2:N							% from bottom to top
% 						h.ax(i) = axes; box on; axis tight;
% 						h.ax(i).Position([1,3,4,2]) = [h.ax(i-1).Position([1,3,4]),sum(h.ax(i-1).Position([2,4]))];
% 						h.ax(i).XTickLabel = '';
% 						h.pl(i) = copyobj(h.pl(i),h.ax(i));
% 					end
% 					axes(h.ax(1));
					% from top to bottom
					h.ax.Position(2) = h.ax.Position(2)+h.ax.Position(4)*(N-1);			% shift first plot to the top
					for i = 2:N
						h.ax(i) = axes; box on; axis tight;
						h.ax(i).Position([1,3,4,2]) = [h.ax(i-1).Position([1,3,4]),-diff(h.ax(i-1).Position([2,4]))];
						h.ax(i).XTickLabel = '';
						h.pl(i) = copyobj(h.pl(i),h.ax(i));
					end
					h.ax(N).XTickLabelMode = 'auto';
					h.ax(1).XTickLabel = '';
					tempDel.delete();
				end
			end
			
			axis tight
			grid on
			xlabel(h.xlbl);
			ylabel(h.ylbl);
			
			% Create Slider(s) inside h.controlPanel
			h.cP_padIn = 5;		% controlPanel padding inside
			h.sld_w = (h.cP_width-h.cP_padIn); h.sld_h = 20;
			for ii = 1:h.nSld
				h.sldText{ii} = uicontrol('Parent',h.controlPanel,'Style','text','String',sprintf('%s %d',h.sldlbl{ii},1),...
										  'HorizontalAlignment','left','FontSize',10,'FontName',h.ax(1).FontName);
				h.sld{ii} = javax.swing.JScrollBar(0,h.sldInitVal(ii),1,h.sldLimits(ii,1),h.sldLimits(ii,2)+1);									%JScrollBar(int orientation, int value, int extent, int min, int max)
				h.SldIdx{ii} = h.sldInitVal(ii)-1;
				[~,h.sldContainer{ii}] = javacomponent(h.sld{ii}, [5, 5,h.sld_w,h.sld_h], h.controlPanel);	% position defined in posDisplay()
				h.sld{ii}.setUnitIncrement(round(diff(h.sldLimits(ii,:))/100)); h.sld{ii}.setBlockIncrement(round(diff(h.sldLimits(ii,:))/10));
				h.hsld{ii} = handle(h.sld{ii},'CallbackProperties');
				set(h.hsld{ii},'AdjustmentValueChangedCallback',@(source,callbackdata) callback_1D_Sld(source, callbackdata,ii));%(source,callbackdata)  set(h.ax,'UserData',round(source.Value)));
			end
			
			posDisplay();
			callback_1D_Sld();
			set(h.f, 'ResizeFcn',@posDisplay);
			
			pl = h.pl;
			
			function callback_1D_Sld(source, callbackdata, n)
				% Callback for 1D Slider plots
				% needs to be sub-function due to shared handle h.
				%  n: number of dimension in z being changed
				if nargin > 0
					h.SldIdx{n} = round(source.Value);
					h.sldText{n}.String = sprintf('%s %d', h.sldlbl{n},round(source.Value));
					h.dataRange = h.SldIdx{1};
					h.zeroPadFactor = h.SldIdx{2};
					calcFFT();
				end
				out = plotSpec();
				for kk = 1:length(h.pl)
					if length(h.pl) == 1
						set(h.pl(1) ,'xdata',h.xdata,'ydata', out);
					else
						if h.normalise && h.distribute
							set(h.pl(kk),'xdata',h.xdata,'ydata', out(:,kk)+length(h.pl)-kk);
						else
							set(h.pl(kk),'xdata',h.xdata,'ydata', out(:,kk));
						end
					end
				end
			end
			
			function calcFFT()
				% performs the FFT. Will be called after each slider adjustment
				maxN = h.dataRange*(1+h.zeroPadFactor);
				if h.power2
					maxN = pow2(nextpow2(maxN));
				end
				h.xdata = (0:maxN-1)/obj.t(2)/maxN;			% get frequency
				if h.useWindowFcn
					if h.useWindowFcn == 1
						window = hann(h.dataRange,'periodic');
					elseif h.useWindowFcn == 2
						window = kaiser(h.dataRange,2.5);
					end
					h.ydata = fft(h.data(1:h.dataRange,:).*(window*ones(1,size(h.data,2))),maxN,1);				% do FFT in 1st dimension
				else
					h.ydata = fft(h.data,maxN,1);			% do FFT in 1st dimension
				end
				if h.FFTshift
					h.ydata = fftshift(h.ydata,1);			% shift to have -f/2 : f/2
					h.xdata = h.xdata - h.xdata(end)/2;
				end
				h.xdata = h.xdata * h.freqScale;
				if h.fInvert
					h.xdata = 1./h.xdata;					% necessary for wavelength plots
				end
				h.xSize = length(h.xdata);
			end
			
			function out = plotSpec()
				% define, which components of the DFT will be plotted, real imag or abs
				switch lower(type)
					case 'linabs'
						% linear absorption as DFT of stateProj autocorrelation function
						out = [real(h.ydata),imag(h.ydata)];
					case 'rhoii'
						% applies DFT to the population probability of rho
						out = abs(h.ydata);
					case 'rhoii-osc-res'
						% applies DFT to the residual of the populations of rho
						out = abs(h.ydata);
					case 'rhoij-imag'
						out = abs(h.ydata);
					case 'chain-n'
						out = abs(h.ydata);
					case 'chain-x2'
						out = abs(h.ydata);
                    case 'chain-x-c'
                        out = abs(h.ydata);
					otherwise
						warning('TDVPData:plotSld1DFT','Plot type not existent!');
						return;
				end
				if h.normalise
					out = bsxfun(@rdivide,out,max(out,[],1));
				end
			end
			
			function posDisplay(varargin)
				if isvalid(h.sldText{1})
					for kk = 1:h.nSld
						h.sldText{kk}.Position      = [h.cP_padIn, h.f.Position(4)-(h.cP_padIn*kk+2*(kk-0.5)*h.sld_h) ,h.sld_w,h.sld_h];
						h.sldContainer{kk}.Position = [h.cP_padIn, h.sldText{kk}.Position(2)-h.sld_h,h.sld_w,h.sld_h];
					end
				end
			end
			
		end
		
		function h   = plotSld2D(obj, type, varargin)
			%% Plot with Slider to flip through different 2D plots
			h = struct();				% struct containing all infos & handles
			h.ts = 1;					% scale time axis
			h.logZ = 0;					% linear z-axis by default
			h.evTocm = 0;				% convert all eV to cm^-1
			h.xdata = (1:obj.para.L).';
			h.ydata = obj.t(1:obj.lastIdx)*h.ts;
			h.zdata = [];					% T x L x NC x ...  if not, needs to be reshaped!
			h.xSize = obj.para.L*ones(1,obj.nChains);	% size for each slider value
			h.imagesc = 0;
			
			if obj.nChains > 1 && obj.para.useTreeMPS
				h.xSize = cellfun(@(x) x.L,obj.para.chain)+1;
				h.xdata = (1:max(h.xSize))'*ones(1,obj.nChains);
			end
			
			h.xlbl = 'Site $k$';
			h.ylbl = '$t$';
			h.zlbl = '';
			h.tlbl = '';
			
			for m = 1:nargin-2
				switch lower(varargin{m})
					case '-fsev'
						% fs timescale for H in eV
						h.ts = 0.658;
						h.ylbl = '$t/\mathrm{fs}$';
						h.ydata = h.ydata*h.ts;
					case '-log'
						% Logarithmic scale in z
						h.logZ = 1;
					case '-cm'
						% use cm^-1 as Units
						h.evTocm = 1;
					case '-imagesc'
						h.imagesc = 1;
				end
			end
			
			% Create figure handles and panels
			h.f = figure();
			h.ax = gca;
			% fix axes size and add space to the right
			h.ax.Units = 'pixels';
			h.f.Position(3) = h.f.Position(3)+100;
			h.controlPanel = uipanel('Tag','Panel','Units','pixels','Position',[h.f.Position(3)-100,1,100,h.f.Position(4)]);
			
			% set normalised units to enable scaling with window size
			h.ax.Units = 'norm';
			h.controlPanel.Units = 'norm';
			
			if ~isempty(strfind(lower(type),'star'))
				% need to find lengths of chains in star picture -> length of x-axis
% 				if strcmp(obj.para.chain{1}.mapping, 'LanczosTriDiag')		% take a representative
% 					h.xSize = cellfun(@(x) length(x.xi), obj.para.chain);
% 				else
					h.xSize = size(obj.omega,1)*ones(1,obj.nChains);
					h.ydata = obj.tS;
					h.ySize = length(h.ydata);
% 				end
			end
			
			switch lower(type)											% all the plot choices
				case 'chain-n'
					h.zdata = real(obj.occC(1:obj.lastIdx,:,:));		% t x L x chain
					h.zlbl = '$\left< n_k \right>$';
					h.sldlbl = {'Chain'};
					h.tlbl = 'Chain Occupation';
				case 'chain-n-d'
					h.zdata = real(obj.occCd(1:obj.lastIdx,:,:,:));		% t x L x state x chain
					h.zlbl = '$\left< n_k \right>$';
					h.sldlbl = {'State','Chain'};
					h.tlbl = 'Chain Occupation - diabatic';
				case 'chain-n-a'
					h.zdata = real(obj.occCa(1:obj.lastIdx,:,:,:));		% t x L x state x chain
					h.zlbl = '$\left< n_k \right>$';
					h.sldlbl = {'State','Chain'};
					h.tlbl = 'Chain Occupation - adiabatic';
				case 'chain-x'
					h.zdata = real(obj.xC(1:obj.lastIdx,:,:));			% t x L x chain
					h.zlbl = '$\left< x_k \right>$';
					h.sldlbl = {'Chain'};
					h.tlbl = 'Chain Displacement';
				case 'chain-x-d'
					h.zdata = real(obj.xCd(1:obj.lastIdx,:,:,:));		% t x L x state x chain
					h.zlbl = '$\left< x_k \right>$';
					h.sldlbl = {'State','Chain'};
					h.tlbl = 'Chain Displacement - diabatic';
				case 'chain-x-a'
					h.zdata = real(obj.xCa(1:obj.lastIdx,:,:,:));		% t x L x state x chain
					h.zlbl = '$\left< x_k \right>$';
					h.sldlbl = {'State','Chain'};
					h.tlbl = 'Chain Displacement - adiabatic';
				case 'chain-x-avg'
					h.zdata = squeeze(sum(real(obj.xCa(1:obj.lastIdx,:,:,:)),3));		% t x L x state x chain
					h.zlbl = '$\left< x_k \right>$';
					h.sldlbl = {'Chain'};
					h.tlbl = 'Chain Displacement';
				case 'chain-x-ft'
					h.maxN = pow2(nextpow2(obj.lastIdx));				% FT window length
					h.ydata = (0:h.maxN/2)/h.ydata(end);				% fs * k/N = k/T where k=0... N/2
					h.ylbl = '$f$ in $1/\mathrm{fs}$';
					if h.evTocm && h.ts ~=1
						h.ydata = h.ydata*4.135*8065.73;
						h.ylbl = '$f/\mathrm{cm}^{-1}$';
					end
					h.zdata = fft(real(obj.xC(1:obj.lastIdx,:,:,:)),h.maxN,1);		% f x L x state x chain
					h.zdata = 2*abs(h.zdata(1:h.maxN/2+1,:,:,:));
					h.zlbl = '$FT(\left< x_k \right>)$';
					h.sldlbl = {'State','Chain'};
					h.tlbl = 'Chain Displacement';
				case 'chain-x2'
					h.zdata = real(obj.xC2(1:obj.lastIdx,:,:));			% t x L x chain
					h.zlbl = '$\left< x^2_k \right>$';
					h.sldlbl = {'Chain'};
					h.tlbl = 'Chain Displacement squared';
				case 'chain-x2-d'
					h.zdata = real(obj.xC2d(1:obj.lastIdx,:,:,:));		% t x L x chain
					h.zlbl = '$\left< x^2_k \right>$';
					h.sldlbl = {'State','Chain'};
					h.tlbl = 'Chain Displacement squared - diabatic';
				case 'chain-x2-a'
					h.zdata = real(obj.xC2a(1:obj.lastIdx,:,:,:));		% t x L x chain
					h.zlbl = '$\left< x^2_k \right>$';
					h.sldlbl = {'State','Chain'};
					h.tlbl = 'Chain Displacement squared - adiabatic';
				case 'chain-x-unc'										% uncertainty in x
					x = real(obj.xC(1:obj.lastIdx,:,:));
					x2 = real(obj.xC2(1:obj.lastIdx,:,:));
					h.zdata = sqrt(x2-x.^2);							% t x L x chain
					h.zdata(:,1,:) = h.zdata(:,1,:) + 0.5;				% balance spin site
					h.zlbl = '$\Delta x_k$';
					h.sldlbl = {'Chain'};
					h.tlbl = 'Chain Displacement - uncertainty';
				case 'star-n'
					h.zdata = real(obj.occS(1:h.ySize,:,:));			% t x w x nChains
					h.zlbl = '$\left< n(\omega) \right>$';
					if ~h.evTocm
						h.xdata = obj.omega;							% w x nChain
						h.xlbl = 'Mode $\omega$ in $eV$';
					else
						h.xdata = obj.omega./1.23984e-4;				% w x nChain
						h.xlbl = 'Mode $\omega$ in $cm^{-1}$';
					end
					h.sldlbl = {'Chain'};
					h.tlbl = 'Star Occupation';
				case 'star-x'
					h.zdata = real(obj.xS(1:obj.lastIdx,:,:,:));		% t x w x state x chain
					h.zlbl = '$\left< f_k \right>$';
					if ~h.evTocm
						h.xdata = obj.omega;							% w x nChain
						h.xlbl = 'Mode $\omega$ in $eV$';
					else
						h.xdata = obj.omega./1.23984e-4;				% w x nChain
						h.xlbl = 'Mode $\omega$ in $cm^{-1}$';
					end
					h.sldlbl = {'State','Chain'};
					h.tlbl = 'Star Displacement';
				case 'state-adiab'
					if ~isempty(obj.sysState)
						h.zdata = obj.sysState(1:obj.lastIdx,:,:).*conj(obj.sysState(1:obj.lastIdx,:,:));			% t x dk x D; calc the probability
					else
						error('Not available, need to extract Observable ss');
						close(h.f);
						return;
					end
					h.zlbl = '$|\Psi_k|^2$';
					h.xlbl = 'Diabatic States';
					h.xdata = (1:size(obj.sysState,2)).';
					h.xSize = size(obj.sysState,2)*ones(1,obj.nChains);
					h.tlbl = 'Adiabatic State evolution';
					h.sldlbl = {'Bond State'};
				case 'heff-offdiag'
					h.zdata = obj.Heff.*conj(obj.Heff);;				% t x D' x dk' x D x dk
					d = size(h.zdata);
					h.zdata = reshape(h.zdata, d(1),d(2)*d(3),[]);		% t x D' * dk' x D * dk
					h.zdata = permute(h.zdata, [2,3,1]);				% D' * dk' x D * dk x t
					for ii = 1:size(h.zdata,3)
						h.zdata(:,:,ii) = h.zdata(:,:,ii) - diag(diag(h.zdata(:,:,ii)));	% remove the offdiagonals
					end
					h.zlbl = '$|H_{eff}|$';
					h.xdata = (1:size(h.zdata,2))';						% D'*dk' x D*dk
					h.xSize = length(h.xdata);
					h.ydata = (1:size(h.zdata,1))';						% D'*dk' x D*dk
					h.ySize = length(h.ydata);
					h.xlbl = 'Initial State D*dk';
					h.ylbl = 'Final State D*dk';
					h.sldlbl = {'t'};
					h.tlbl = 'Effective Hamiltonian paths';
				case 'heff-current'
					h.zdata = obj.getData('heff-current');				% t x D' x dk' x D x dk
					d = size(h.zdata);
					h.zdata = reshape(h.zdata, d(1),d(2)*d(3),[]);		% t x D' * dk' x D * dk
					h.zlbl = '$\left< j \right>$';
					h.xdata = (1:size(h.zdata,2))';						% D'*dk' x D*dk
					h.xSize = length(h.xdata);
					h.xlbl = 'Final State D*dk';
					h.sldlbl = {'In. State'};
					h.tlbl = 'Currents';
				case 'heff-current-v2'
					h.zdata = obj.getData('heff-current');				% t x D' x dk' x D x dk
					d = size(h.zdata);
					h.zdata = reshape(h.zdata, d(1),d(2)*d(3),[]);		% t x D' * dk' x D * dk
					h.zdata = permute(h.zdata,[3,2,1]);
					h.zlbl = '$\left< j \right>$';
					h.xdata = (1:size(h.zdata,2))';						% D'*dk' x D*dk
					h.ydata = h.xdata;
					h.xSize = length(h.xdata);
					h.xlbl = 'Initial State D*dk';
					h.ylbl = 'Final State D*dk';
					h.sldlbl = {'t_i'};
					h.tlbl = 'Currents';
				otherwise
					return;
			end
			
			
			h.f.Name = h.tlbl;
			
			if h.logZ
				h.zdata = sign(h.zdata) .* log10(abs(h.zdata));
				h.zlbl = ['$\log_{10}',h.zlbl(2:end)];
			end
			
			h.nSld = ndims(h.zdata)-2;					% number of sliders needed
			h.SldIdx = num2cell(ones(1,h.nSld),1);		% indices for each slider dimension, indicating currently displayed slice
			h.ysize = size(h.ydata);
			h.zsize = size(h.zdata);
			
			if ~h.imagesc
				h.pl = surf(h.xdata(1:h.xSize(1),1),h.ydata,h.zdata(:,1:h.xSize(1),h.SldIdx{:}));
				shading interp
				rotate3d on
			else
				h.pl = imagesc(h.xdata(1:h.xSize(1),1),h.ydata,h.zdata(:,1:h.xSize(1),h.SldIdx{:}));
				h.ax.YAxis.Direction = 'normal';
			end
			axis tight
			h.ax.View = [0 90];
			h.ax.TickDir = 'out';
			xlabel(h.xlbl);
			ylabel(h.ylbl);
			zlabel(h.zlbl);
			h.cb = colorbar; h.cb.Title.Interpreter = 'latex';
			h.cb.Title.String = h.zlbl;

			% Create Slider(s) inside h.controlPanel
			h.sld_w = 70; h.sld_h = 20;
			h.cP_padIn = 5;		% controlPanel padding inside
			for ii = 1:h.nSld
				h.sldText{ii} = uicontrol('Parent',h.controlPanel,'Style','text','String',sprintf('%s %d',h.sldlbl{ii},1),...
										  'HorizontalAlignment','left','FontSize',10,'FontName',h.ax.FontName);
				h.sld{ii} = javax.swing.JScrollBar(0,1,1,1,h.zsize(2+ii)+1);		%JScrollBar(int orientation, int value, int extent, int min, int max)
				[~,h.sldContainer{ii}] = javacomponent(h.sld{ii}, [5, 5,h.sld_w,h.sld_h], h.controlPanel);	% position defined in posDisplay()
				h.sld{ii}.setUnitIncrement(1); h.sld{ii}.setBlockIncrement(3);
				h.hsld{ii} = handle(h.sld{ii},'CallbackProperties');
				set(h.hsld{ii},'AdjustmentValueChangedCallback',@(source,callbackdata) callback_2D_Sld(source, callbackdata,ii));%(source,callbackdata)  set(h.ax,'UserData',round(source.Value)));
			end
			posDisplay();
			set(h.f, 'ResizeFcn',@posDisplay);
			
			function callback_2D_Sld(source, callbackdata, n)
				% Callback for 2D Slider plots
				% needs to be sub-function due to shared handle h.
				%  n: number of dimension in z being changed
				h.SldIdx{n} = round(source.Value);
				h.sldText{n}.String = sprintf('%s %d', h.sldlbl{n},round(source.Value));
% 				fprintf('%d',h.SldIdx{:});
				if size(h.xdata,2) > 1 && size(h.xdata,2) == size(h.zdata,n+2)
					% high probability that current slider is going through chains
					mc = round(source.Value);
					set(h.pl,'xdata', h.xdata(1:h.xSize(mc),mc));
					set(h.pl,'zdata', h.zdata(:,1:h.xSize(mc),h.SldIdx{:}));
					h.ax.CLim = h.ax.ZLim;
				else
					xs = length(get(h.pl,'xdata'));
					if isa(h.pl, 'matlab.graphics.primitive.Image')
						set(h.pl,'cdata', h.zdata(:,1:xs,h.SldIdx{:}));
% 						h.ax.CLim = [min(h.pl.CData(:)), max(h.pl.CData(:))];
% 						set(h.pl,'cdatamapping','scaled');					% refresh Clim
					else
						set(h.pl,'zdata', h.zdata(:,1:xs,h.SldIdx{:}));
						h.ax.CLim = h.ax.ZLim;
					end
				end
	% 			fprintf('Source.Value: %g, Sld.Value: %g',source.Value,h.sld.getValue);
	% 			callbackdata
			end
			
			function posDisplay(varargin)
				for kk = 1:h.nSld
					h.sldText{kk}.Position      = [h.cP_padIn, h.f.Position(4)-(h.cP_padIn*kk+2*(kk-0.5)*h.sld_h) ,h.sld_w,h.sld_h];
					h.sldContainer{kk}.Position = [h.cP_padIn, h.sldText{kk}.Position(2)-h.sld_h,h.sld_w,h.sld_h];
				end
			end
			
		end
		
		function h   = plotSldBloch(obj, type, varargin)
			h = struct();
			h.xlbl = '$\left<sigma_x\right>$';
			h.ylbl = '$\left<sigma_y\right>$';
			h.zlbl = '$\left<sigma_z\right>$';
			h.tlbl = '';
			h.linewidth = 5;
			if isempty(obj.rho)
				error('need reduced density matrix for this function!');
			end
			assert(size(obj.rho,2) == 2, 'Can only do this for 2-lvl systems right now');
			
			% Observables to compute
			[sigmaX,sigmaY,sigmaZ] = spinop('Z');
			h.sigmaX = reshape(sigmaX.',[],1);
			h.sigmaY = reshape(sigmaY.',[],1);
			h.sigmaZ = reshape(sigmaZ.',[],1);
			
			% reshape rho into convenient shape
			h.d = size(obj.rho);		% t x i x j x k (states)
			if length(h.d) == 4
				h.Nstates = h.d(4);
				h.rho = reshape(obj.rho,h.d(1),[],h.d(4));
				h.sv = zeros(h.d(1),h.Nstates-1);
				h.tlbl = 'Adiabatic Spin evolution';
			else
				h.Nstates = 1;
				h.rho = reshape(obj.rho,h.d(1),[]);
				h.sv = [];			% not used since no sv existing
				h.tlbl = 'Spin Evolution';
			end
			
			for ii = 2:h.Nstates
				% extract the SV^2 to plot bloch vectors with correct length
				h.sv(:,ii-1) = real(diag(squeeze(h.rho(:,:,1))*squeeze(h.rho(:,:,ii))'));	% achieves the same
			end
			
			h.xdata = zeros(h.d(1),h.Nstates);
			h.xsize = size(h.xdata);
			h.ydata = zeros(h.d(1),h.Nstates);
			h.ysize = size(h.ydata);
			h.zdata = zeros(h.d(1),h.Nstates);
			h.zsize = size(h.zdata);
			for ii = 1:h.Nstates
				% compute <sigma> of all rho
				h.xdata(:,ii) = real(squeeze(h.rho(:,:,ii))*h.sigmaX);
				h.ydata(:,ii) = real(squeeze(h.rho(:,:,ii))*h.sigmaY);
				h.zdata(:,ii) = real(squeeze(h.rho(:,:,ii))*h.sigmaZ);
			end
			
			% Create figure handles and panels
			h.f = figure();
			h.f.Name = h.tlbl;
			h.ax = gca; hold all;
			% fix axes size and add space to the right
			h.ax.Units = 'pixels';
			h.f.Position(3) = h.f.Position(3)+100;
			h.controlPanel = uipanel('Tag','Panel','Units','pixels','Position',[h.f.Position(3)-100,1,100,h.f.Position(4)]);
% 			set(h.f,'DefaultLineLineWidth',h.linewidth);
			h.ax.Color = 'none';

			
			% set normalised units to enable scaling with window size
			h.ax.Units = 'norm';
			h.controlPanel.Units = 'norm';
			
			% Slider options
			h.nSld = 1;					% Use slider only in t
			h.SldIdx = num2cell(ones(1,h.nSld),1);		% indices for each slider dimension, indicating currently displayed slice
			h.sldlbl = {'Time'};
			
			% Create Plot for t = 0
			h.sphereon = true;
			if h.sphereon
				sphere
				daspect([1 1 1])
				alpha(0.2)
% 				h.ax.Children.Linestyle = ':';			% thinner mesh
				set(get(gca,'children'),'linestyle',':','linewidth',0.5)
			end
% 			scatter3(tresults.spin.sx,tresults.spin.sy,tresults.spin.sz,20,'filled');
			h.pl = cell(h.Nstates,1);
			for ii = 1:h.Nstates
				if ii == 1
% 					h.pl{ii} = plot3([0,h.xdata(h.SldIdx{:},ii)],[0,h.ydata(h.SldIdx{:},ii)],[0,h.zdata(h.SldIdx{:},ii)]);
					h.pl{ii} = quiver3(0,0,0,h.xdata(h.SldIdx{:},ii),h.ydata(h.SldIdx{:},ii),h.zdata(h.SldIdx{:},ii),1,...
								'LineWidth',h.linewidth);
				else
					h.pl{ii} = quiver3(0,0,0,0,0,0,1,'LineWidth',h.linewidth*0.5);
% 					h.pl{ii} = plot3([0,0],[0,0],[0,0],'LineWidth',h.linewidth*0.5);		% dummy initialise! assume always pure initial state
				end
			end
			h.ax.XLim = [-1,1];
			h.ax.YLim = [-1,1];
			h.ax.ZLim = [-1,1];
			h.ax.View = [-29,16];
			rotate3d on
			
			%
			% Create Slider(s) inside h.controlPanel
			h.sld_w = 70; h.sld_h = 20;
			h.cP_padIn = 5;		% controlPanel padding inside
			for ii = 1:h.nSld		% only have 1 slider!!
				h.sldText{ii} = uicontrol('Parent',h.controlPanel,'Style','text','String',sprintf('%s %d',h.sldlbl{ii},1),...
										  'HorizontalAlignment','left','FontSize',10,'FontName',h.ax.FontName);
				h.sld{ii} = javax.swing.JScrollBar(0,1,1,1,h.xsize(1)+1);		%JScrollBar(int orientation, int value, int extent, int min, int max)
				[~,h.sldContainer{ii}] = javacomponent(h.sld{ii}, [5, 5,h.sld_w,h.sld_h], h.controlPanel);	% position defined in posDisplay()
				h.sld{ii}.setUnitIncrement(1); h.sld{ii}.setBlockIncrement(3);
				h.hsld{ii} = handle(h.sld{ii},'CallbackProperties');
				set(h.hsld{ii},'AdjustmentValueChangedCallback',@(source,callbackdata) callback_2D_Sld(source, callbackdata,ii));%(source,callbackdata)  set(h.ax,'UserData',round(source.Value)));
			end
			posDisplay();
			set(h.f, 'ResizeFcn',@posDisplay);
			
			function callback_2D_Sld(source, callbackdata, n)
				% Callback for 2D Slider plots
				% needs to be sub-function due to shared handle h.
				%  n: number of dimension in z being changed
				h.SldIdx{n} = round(source.Value);
				h.sldText{n}.String = sprintf('%s %d', h.sldlbl{n},round(source.Value));
% 				fprintf('%d',h.SldIdx{:});
				tt = round(source.Value);
				for kk = 1:h.Nstates
					if kk == 1
						% update the full vector
						h.pl{1}.UData = h.xdata(tt,1);
						h.pl{1}.VData = h.ydata(tt,1);
						h.pl{1}.WData = h.zdata(tt,1);
					else
						% update the part vectors
						% get their origin from the previous vector
						if kk == 2
% 							x0 = 0; y0 = 0; z0 = 0;
							% do nothing
						else
							% set origin to end of previous quiver
							h.pl{kk}.XData = h.pl{kk-1}.XData+h.pl{kk-1}.UData;
							h.pl{kk}.YData = h.pl{kk-1}.YData+h.pl{kk-1}.VData;
							h.pl{kk}.ZData = h.pl{kk-1}.ZData+h.pl{kk-1}.WData;
						end
						% calculate their destination point
						h.pl{kk}.UData = h.xdata(tt,kk)*h.sv(tt,kk-1);
						h.pl{kk}.VData = h.ydata(tt,kk)*h.sv(tt,kk-1);
						h.pl{kk}.WData = h.zdata(tt,kk)*h.sv(tt,kk-1);
% 						set(h.pl{kk},'xdata', [x0,x1]);
% 						set(h.pl{kk},'ydata', [y0,y1]);
% 						set(h.pl{kk},'zdata', [z0,z1]);
					end
				end
	% 			fprintf('Source.Value: %g, Sld.Value: %g',source.Value,h.sld.getValue);
	% 			callbackdata
			end
			
			function posDisplay(varargin)
				for kk = 1:h.nSld
					h.sldText{kk}.Position      = [h.cP_padIn, h.f.Position(4)-(h.cP_padIn*kk+2*(kk-0.5)*h.sld_h) ,h.sld_w,h.sld_h];
					h.sldContainer{kk}.Position = [h.cP_padIn, h.sldText{kk}.Position(2)-h.sld_h,h.sld_w,h.sld_h];
				end
			end
			
		end
		
		function out = chainParaInfo(obj)
			% returns the parameters of each chain.
			out = [];
			for ii = 1:obj.para.nChains
				out(ii).xi = obj.para.chain{ii}.xi;
				out(ii).Gamma = obj.para.chain{ii}.Gamma;
				out(ii).epsilon = obj.para.chain{ii}.epsilon;
				out(ii).t = obj.para.chain{ii}.t;
			end
		end
		
		function exportCSV(obj,type,loc,varargin)
			% function used for exporting dataset of 'type' in obj-array to location loc: cell array
			%
			if length(obj) > 1
				for ii = 1:length(obj)
					obj(ii).exportCSV(type,loc{ii});
				end
				return;
			end
			
			% from now: loc is a string, obj is single object
			
			switch lower(type)
				case 'sz'
					x = reshape(obj.t(1:obj.lastIdx),[],1);
					y = reshape(obj.tresults.spin.sz(1:obj.lastIdx),[],1);
					out = [x,y];
				case 'chain-n'
					x = reshape(1:obj.para.L,[],1);
					y = reshape(obj.t(1:obj.lastIdx),[],1);
					z = obj.occC(1:obj.lastIdx,:,1);
					out = [0,x';y,z];
				case 'star-n'
					x = reshape(obj.omega,[],1);
					y = reshape(obj.tS,[],1);
					z = obj.occS(:,:,1);								% t x w x nChains
					% only pick first chain for now.
					out = [0,x';y,z];
				case 'star-x'
					x = reshape(obj.omega,[],1);
					y = reshape(obj.tS,[],1);
					z = real(obj.xS(:,:,1));								% t x w x nChains
					% only pick first chain for now.
					out = [0,x';y,z];
				otherwise
					error('Dataset type not supported yet!');
			end
			dir = fileparts(loc);
			if ~exist(dir,'dir')
				mkdir(dir);
			end
			csvwrite(loc,out);
			
        end
        
        function out = getSysEnvObs(obj,type,varargin)
			%% function out = getSysEnvObs(obj,type,varargin)
			%	Contract the environment states with an operator specified by
			%	type:
			%	 'n' : occupation
			%    'x' : displacement
			%
			%	Only for bosons on the 2nd site for now.
			%	
			%	out: NC x 1 cell
			%		with Op_(r',r) right observable containing the contracted operator
			NC = obj.para.nChains;
			out = cell(obj.lastIdx, NC);		% t x NC
			for mc = 1:NC
				% compute for all times simultaneously
				dk = cellfun(@(x) size(x{mc},1), obj.Vmat(1:obj.lastIdx,2));	% dk(t)
				bp = arrayfun(@(x) bosonop(x,0,'n'),dk,'UniformOutput',false);
				switch lower(type)
					% select the appropriate operator
					case 'n'
						op = cellfun(@(x) x*x', bp,'UniformOutput',false);
					case 'x'
						op = cellfun(@(x) (x+x')./2, bp,'UniformOutput',false);
				end
				out(:,mc) = arrayfun(@(t) updateCright([],obj.mps{t,2}{mc}, obj.Vmat{t,2}{mc},op{t},obj.mps{t,2}{mc}, obj.Vmat{t,2}{mc}),...
					(1:obj.lastIdx)','UniformOutput',false);
				
			end
		end
		
		function obj = getSpin(obj)
			h = struct();
			% Observables to compute
			[sigmaX,sigmaY,sigmaZ] = spinop('Z');
			h.sigmaX = reshape(sigmaX.',[],1);
			h.sigmaY = reshape(sigmaY.',[],1);
			h.sigmaZ = reshape(sigmaZ.',[],1);
			
			% reshape rho into convenient shape
			h.d = size(obj.rho);		% t x i x j x k (states)
			if length(h.d) == 4
				h.Nstates = h.d(4);
				h.rho = reshape(obj.rho,h.d(1),[],h.d(4));
				h.sv = zeros(h.d(1),h.Nstates-1);
				h.tlbl = 'Adiabatic Spin evolution';
			else
				h.Nstates = 1;
				h.rho = reshape(obj.rho,h.d(1),[]);
				h.sv = [];			% not used since no sv existing
				h.tlbl = 'Spin Evolution';
			end
			
			obj.spin = zeros(obj.lastIdx,3,h.Nstates,'single');
			for ii = 1:h.Nstates
				% compute <sigma> of all rho
				obj.spin(:,1,ii) = real(squeeze(h.rho(1:obj.lastIdx,:,ii))*h.sigmaX);
				obj.spin(:,2,ii) = real(squeeze(h.rho(1:obj.lastIdx,:,ii))*h.sigmaY);
				obj.spin(:,3,ii) = real(squeeze(h.rho(1:obj.lastIdx,:,ii))*h.sigmaZ);
			end
		end
		
		function out = getRDMDiff(obj, obj2)
			%% out = getRDMDiff(obj, obj2)
			% calculates the norm difference between two t-RDM
			% helpful for giving convergence estimates
			%
			assert(all(size(obj.rho) == size(obj2.rho)),'Please use comparable simulations with same length and timesteps');
			d = size(obj.rho);
% 			out = reshape(obj.rho,d(1),[]) - reshape(obj2.rho,d(1),[]);		% t x dk*dk
% 			out = 1 - reshape(obj2.rho,d(1),[])./reshape(obj.rho,d(1),[]);		% t x dk*dk
			out = diff(reshape(obj.rho,d(1),[])) - diff(reshape(obj2.rho,d(1),[]));		% t x dk*dk
			out = sqrt(sum(out.*conj(out)./d(2)^2,2));
		end
		
		function out = getDkUsage(obj,varargin)
			%% getDkUsage()
			%	get the maximum number of relevantly occupied Dk for each boson
			%	only applicable to treeMPS
			assert(~isempty(obj.treeMPS),'No treeMPS stored, unable to give answer!');
			out = TDVPData.dkUsage(obj.treeMPS);						% forward to the static function
		end
	end
	
	% Static methods are functions that do not require 'obj' but are
    % relevant to methods that do.
	% Also some helper functions
    % They are called using TDVPData.myFun()
    methods(Static = true)
		function foldLib = getTDVPLib()
			filePattern = '20*\results-Till*.mat';
			foldLib = [rdir(['..\TDVP\',filePattern]);rdir(['..\TDVP\Data\',filePattern]);rdir(['..\cacheComputations\',filePattern]);...
				       rdir(['..\TDVP-Git\',filePattern]);rdir(['..\TDVP-Git\Data\',filePattern])];
			% remove 'incomplete' files
			foldLib = foldLib(arrayfun(@(x) isempty(strfind(x.name,'incomplete')),foldLib));
		end
		
		function matches = getMatches(lib, dirPat, filPat)
			%% function matches = getMatches(lib, dirPat, filPat)
			%	find files with regexp
			matches = lib(arrayfun(@(x) ~isempty(regexp(x.name,[dirPat,'\\',filPat],'once')),lib));
		end
		
		function B = movAvg(A,m)
			%% function A = movAvg(A,m)
			%	computes the moving average with window size m
			%	smooth along first dimension
			d = size(A);
			B = reshape(A,d(1),[]);
			for i = 1:size(B,2)
				B(:,i) = smooth(B(:,i),m,'moving');
			end
% 			figure(); hold all
% 			plot(reshape(A,d(1),[]),'b')
% 			plot(B,'r')
			B = reshape(B,d);
		end
		
		function B = movAvgRes(A,m)
			%% function A = movAvgRes(A,m)
			%	computes the residual of a moving average with window size m
			B = A - TDVPData.movAvg(A,m);
		end
		
		function B = meanRes(A)
			%% function B = meanRes(A,m)
			%	subtracts the mean along the first dimension
			B = bsxfun(@minus,A,mean(A,1));
		end
		
		function y = medianFilter(x,windowWidth)
            % function imOut = medianFilter(imIn,windowWidth)
            %
			% Applies Median filter to 2D Matrices only
			% windowWidth is [m n] indicating widths in both directions
            % The original code was posted on matlab central file exchange.
            % I'm unable to find the original to give proper credits.
            % This version was slightly modified by Florian Schroeder, April 2017.
            
            if mod(windowWidth, 2 ) == 0
                disp('Window has even size');
                return
			end
			
			d = size(x);
            y = x;
            
            windowBreadth = (windowWidth - 1)/2;
            
            for m = 1:d(1)
                for n = 1:d(2)
                    r1 = max(m-windowBreadth(1),1):min(m+windowBreadth(1),d(1));
                    r2 = max(n-windowBreadth(2),1):min(n+windowBreadth(2),d(2));
                    t1 = x(r1,r2);
                    t2 = median(t1(:));
                    y(m,n) = t2;
                end
            end
            
		end
		
		function y = maxFilter(x,windowWidth)
            % function imOut = maxFilter(imIn,windowWidth)
            %
			% Applies Maximum filter to 2D Matrices
			% windowWidth is [m n] indicating widths in both directions
            % The original code was posted on matlab central file exchange.
            % I'm unable to find the original to give proper credits.
            % This version was slightly modified by Florian Schroeder, May 2017.
            
            if mod(windowWidth, 2 ) == 0
                disp('Window has even size');
                return
			end
			
			d = size(x);
            y = x;
            
            windowBreadth = (windowWidth - 1)/2;
            
            for m = 1:d(1)
                for n = 1:d(2)
                    r1 = max(m-windowBreadth(1),1):min(m+windowBreadth(1),d(1));
                    r2 = max(n-windowBreadth(2),1):min(n+windowBreadth(2),d(2));
                    t1 = x(r1,r2);
                    t2 = max(t1(:));
                    y(m,n) = t2;
                end
            end
            
		end
		
		function resetColorOrder(ax)
			
			if isempty(ax)
				ax = gca;
			end
			
			col = get(0,'defaultAxesColorOrder');
			col = mat2cell(col,ones(7,1),3);
			
			pl_handle = [ax.Children];
			idx = fliplr(1:length(pl_handle));
			idx = mod(idx-1,7)+1;
			set(pl_handle,{'Color'},col(idx,:));
			
		end
		
		function h = plotGrid(rows,cols,varargin)
			%% function h = plotGrid(rows,cols)
			%	creates a grid of axes
			%	onto a 8.5 cm x 6.4 cm figure suitable for 2-column format.
			%
			%	h = plotGrid(rows,cols, 'rowheight',h, 'rowwidth', w)
			%		allows specifying different values in cm. Default: h = 1.0, w = 7.2;
			%
			%	h = plotGrid(rows,cols, figHandle,...)
			%		creates the grid in defined figure 
			%
			p = inputParser;
			addRequired(p,'rows',@isnumeric);
			addRequired(p,'cols',@isnumeric);
			addOptional(p,'f',[],@(x) isa(x,'matlab.ui.Figure'));	%||isa(x,'matlab.graphics.axis.Axes')); would be nice to include for substacking
			addParameter(p,'rowheight',1.0,@isnumeric);
			addParameter(p,'rowwidth' ,7.2,@isnumeric);
			parse(p,rows,cols,varargin{:});
			
			if ~isempty(p.Results.f)
				h.f = p.Results.f;
				figure(h.f);
				clf;
			else
				h.f = figure();
			end
			h.rows = rows;
			h.cols = cols;			% TODO: extend for more columns
			
			% important measures:
				% 2016a: ax.Position = [1.1041   0.70432   6.5822   5.2184]   for only left and bottom axis label
			h.axW = p.Results.rowwidth/h.cols; h.axH = p.Results.rowheight; 
			h.axdX = 0; h.axdY = 0; h.axdYgroup = 0.3;	% maximum sizes and default paddings between axes
			h.axX = 1.1; h.axY = 0.9;		% spacing bottom left for ticks and labels
			
			% Pre-process axes
			h.f.Units = 'centimeters'; h.f.Position([3,4]) = [1.3+p.Results.rowwidth,1+h.axH*h.rows];
			% f.Renderer = 'painters';
			set(h.f,'DefaultAxesFontSize', 8,...
					'DefaultLegendInterpreter','latex',...
					'DefaultLineLineWidth',1);
% 			% Lowest 
% 			i=h.rows*h.cols; h.ax(i) = axes; box on;
% 			xlabel('xlabel');ylabel('ylabel');
% 			h.ax(i).Units = 'centimeters'; h.ax(i).Position([2,3,4]) = [0.9,h.axW,h.axH];
			
			for kk = 1:h.rows*h.cols			% linear index
				[ii,jj] = ind2sub([h.rows,h.cols],kk);
				h.ax(kk) = axes; box on;
				h.ax(kk).Units = 'centimeters'; 
				h.ax(kk).Position([1,2,3,4]) = [h.axX+(jj-1)*h.axW,h.axY+(h.rows-ii)*h.axH,h.axW,h.axH];
				if ii ~= h.rows
					h.ax(kk).XTickLabel = '';
				elseif jj == 1
					xlabel('xlabel');ylabel('ylabel');
				end
				if jj ~= 1
					h.ax(kk).YTickLabel = '';
				end
			end
			[h.ax.Units] = deal('norm');
			h.f.Units = 'norm';
		end
		
		function h = plotRowGroups(rows,varargin)
			%% function h = plotRowGroups(rows[,f,varargin])
			%	creates rows of axes
			%	onto a 8.5 cm x ? cm figure
			%	rows = [5,7]
			%		gives 2 groups of rows
			%	rows = [1,2,3]
			%		give 3 groups of rows
			% Name-Value Pair arguments:
			%	'rowheight'		1x1 real
			
			p = inputParser;
			addRequired(p,'rows',@isnumeric);
			addOptional(p,'f',[],@(x) isa(x,'matlab.ui.Figure'));	%||isa(x,'matlab.graphics.axis.Axes')); would be nice to include for substacking
			addParameter(p,'rowheight',1.0,@isnumeric);
			addParameter(p,'rowwidth' ,7.2,@isnumeric);
			parse(p,rows,varargin{:});
			
			if ~isempty(p.Results.f)
				h.f = p.Results.f;
				figure(h.f);
				clf;
			else
				h.f = figure();
			end
			
			h.rows = rows;
			h.split = cumsum(rows); h.split = h.split(1:end-1);	% which axes have to be shifted upwards for the split
			% important measures:
				% 2016a: ax.Position = [1.1041   0.70432   6.5822   5.2184]   for only left and bottom axis label
			h.axW = p.Results.rowwidth; h.axH = p.Results.rowheight; h.axdX = 0; h.axdY = 0; h.axdYgroup = 0.3;	% maximum sizes and default paddings between axes
			
			% Pre-process axes
			h.f.Units = 'centimeters'; h.f.Position([3,4]) = [8.5,1+h.axH*sum(h.rows)+h.axdYgroup*(length(h.rows)-1)];
			% f.Renderer = 'painters';
			set(h.f,'DefaultAxesFontSize', 8,...
					'DefaultLegendInterpreter','latex',...
					'DefaultLineLineWidth',1);
			% Lowest 
			i=sum(h.rows); h.ax(i) = axes; box on;
			xlabel('xlabel');ylabel('ylabel');
			h.ax(i).Units = 'centimeters'; h.ax(i).Position([2,3,4]) = [0.9,h.axW,h.axH];

			for i = (sum(h.rows)-1):-1:1
				h.ax(i) = axes; box on;
				h.ax(i).Units = 'centimeters'; h.ax(i).Position([1,3,4,2]) = [h.ax(i+1).Position([1,3,4]),sum(h.ax(i+1).Position([2,4]))];
				h.ax(i).XTickLabel = '';
				if any(i == h.split)
					h.ax(i).Position(2) = h.ax(i).Position(2)+h.axdYgroup;		% shift upwards for group split
				end
			end

			[h.ax.Units] = deal('norm');
			movegui(h.f,'center');
		end
		
		function h = FFT(h)
			%% function h = FFT(h)
			%	performs the FFT using the periodogram function
			%	h.data = x;		-> h.ydata = FFT(x)
			%	h.tdata = t;	-> h.xdata = freq
			if ~isfield(h,'reassign')
				h.reassign = 0;
			end
			if ~isfield(h,'window')
				h.window = '';
			end
			if ~isfield(h,'zeroPadFact')
				h.zeroPadFact = 0;
			end
			assert(length(h.tdata) == size(h.data,1),'FFT input has wrong dimensions');
			
			h.m = length(h.tdata);
			
			if strcmp(h.window,'hann')
				h.wind = hann(h.m,'periodic');
			elseif strcmp(h.window,'kaiser')
				h.wind = kaiser(h.m,5);
			else
				h.wind = [];
			end
			
			if h.reassign
				[pre,ft,pxx,fc] = periodogram(h.data,h.wind,pow2(nextpow2((1+h.zeroPadFact)*h.m)),h.Fs,'power','reassigned');
% 				h.xdata = fc;		% center-of-mass frequencies, plots only for stem()
				h.xdata = ft;
				h.ydata = pre;
			else
				[pxx,f] = periodogram(h.data,h.wind,pow2(nextpow2((1+h.zeroPadFact)*h.m)),h.Fs,'power');
				h.xdata = f;
				h.ydata = pxx;
			end
			
		end
		
		function h = tFFT(h)
			%% function h = tFFT(h)
			%	performs the time-resolved FFT using the spectrogram function
			%	h.data = x;		-> h.ydata = FFT(x)
			%	h.tdata = t;	-> h.xdata = freq
			if ~isfield(h,'reassign')
				h.reassign = 0;
			end
			if ~isfield(h,'window')
				h.window = '';
			end
			if ~isfield(h,'zeroPadFact')
				h.zeroPadFact = 0;
			end
			assert(length(h.tdata) == size(h.data,1),'FFT input has wrong dimensions');
			
			h.m = length(h.tdata);
			
			if strcmp(h.window,'hann')
				h.wind = hann(h.m,'periodic');
			elseif strcmp(h.window,'kaiser')
				h.wind = kaiser(h.m,5);
			else
				h.wind = [];
			end
			
			if h.reassign
								  spectrogram(h.data(:,4),kaiser(1024,5),1000,(h.m*7),10,'reassigned','power','MinThreshold',-40)
				[pre,ft,pxx,fc] = spectrogram(h.data,h.wind,pow2(nextpow2((1+h.zeroPadFact)*h.m)),h.Fs,'power','reassigned');
% 				h.xdata = fc;		% center-of-mass frequencies, plots only for stem()
				h.xdata = ft;
				h.ydata = pre;
			else
				[pxx,f] = spectrogram(h.data,h.wind,pow2(nextpow2((1+h.zeroPadFact)*h.m)),h.Fs,'power');
				h.xdata = f;
				h.ydata = pxx;
			end
			
		end
		
		function out = printSvNETreeMPS(treeMPS)
			% this function can only be applied to treeMPS to obtain the SvNE along each bond
			
			% 1. retrieve SvNE to parent or system
			out = sprintf('SvNE([%s]):',num2str(treeMPS.treeIdx,'%-4u'));
			if treeMPS.level == 1
				% get from svd on MPS matrix
				A = reshape(treeMPS.mps{1},[],treeMPS.dk);
				[~,s,~] = svd2(A);
			else
				s = diag(treeMPS.Amat_sv{1});
			end
			out = sprintf('%s %.2f\n',out,vonNeumannEntropy(s));
			
			% 2. concatenate results from children
			if treeMPS.degree > 0
				for ii = 1:treeMPS.degree
					out = [out, TDVPData.printSvNETreeMPS(treeMPS.child(ii))];
				end
			end
		end
		
		function pl = plotVariance(x,y,std,ylim,ax,varargin)
			% plot a line defined by x,y with shaded variance given by std
			% plot into axes given in ax
			% need to give expected ylim = [min,max] in order to determine correct line thickness
			% if giving subshades then: sum(subshades,2) = std
			
			colOrder = get(0,'DefaultaxesColorOrder');	% only used if std ~= n by 1
			
			% process varargin with parser
			p = inputParser;
			addParameter(p,'thickness',0.02,@isnumeric);
			addParameter(p,'alpha' ,0.5,@isnumeric);
			addParameter(p,'subshades',[],@isnumeric);		% t x n array for n subshades to plot
			parse(p,varargin{:});
			
			max_thickness = p.Results.thickness;						% now: 2% of yrange should be maximum std thickness
			fill_alpha = p.Results.alpha;
			
			% make all row vectors
			x   = reshape(x,1,[]); 
			y   = reshape(y,1,[]);
			std = reshape(std,1,[]);
			upper = y + diff(ylim)*std/2 * max_thickness;
			lower = y - diff(ylim)*std/2 * max_thickness;
			
			pl(1) = plot(ax, x,y);		% first plot the mean
			if isempty(p.Results.subshades)
				if verLessThan('matlab','9.3')
					pl(2) = fill([x,x(end:-1:1)], [upper,lower(end:-1:1)], pl(1).Color*fill_alpha + (1-fill_alpha)*[1,1,1]);
				else
					pl(2) = fill([x,x(end:-1:1)], [upper,lower(end:-1:1)], pl(1).Color);
					pl(2).FaceAlpha = fill_alpha;
				end
				uistack(pl(2), 'bottom');
				set(pl(2),'EdgeColor','none');
				TDVPData.decrColorOrderIndex(ax);
			else
				std = p.Results.subshades';			% n x t
				for ii = 1:size(std,1)
					upper = lower + diff(ylim)*std(ii,:) * max_thickness;
					if max(std(ii,:)) > 5e-3
						% Only plot if max patch height > 1%
						if verLessThan('matlab','9.3')
							pl(ii+1) = fill([x,x(end:-1:1)], [upper,lower(end:-1:1)], colOrder(ii,:)*fill_alpha + (1-fill_alpha)*[1,1,1]);
						else
							pl(ii+1) = fill([x,x(end:-1:1)], [upper,lower(end:-1:1)], colOrder(ii,:));
							pl(ii+1).FaceAlpha = fill_alpha;
						end
						uistack(pl(ii+1),'bottom');
						set(pl(ii+1),'EdgeColor','none')
						if ax.ColorOrderIndex == 1
							ax.ColorOrderIndex = 7;
						else
							ax.ColorOrderIndex = ax.ColorOrderIndex-1;
						end
					end
					lower = upper;
				end
			end
		end
		
		function decrColorOrderIndex(ax)
			%% decrColorOrderIndex(ax)
			%	decreases the ColorOrderIndex of ax by one
			if ax.ColorOrderIndex == 1
				ax.ColorOrderIndex = 7;
			else
				ax.ColorOrderIndex = ax.ColorOrderIndex-1;
			end
		end
		
		function pl = plotVarianceLines(x,y,std,ylim,ax,varargin)
			% plot a line defined by x,y with shaded variance given by std
			% plot into axes given in ax
			% need to give expected ylim = [min,max] in order to determine correct line thickness
			% if giving subshades then: sum(subshades,2) = std
			
			colOrder = get(0,'DefaultaxesColorOrder');	% only used if std ~= n by 1
			
			% process varargin with parser
			p = inputParser;
			addParameter(p,'thickness',0.02,@isnumeric);
			addParameter(p,'alpha' ,0.5,@isnumeric);
			addParameter(p,'subshades',[],@isnumeric);		% t x n array for n subshades to plot
			parse(p,varargin{:});
			
			max_thickness = p.Results.thickness;						% now: 2% of yrange should be maximum std thickness
			fill_alpha = p.Results.alpha;
			
			% make all row vectors
			x   = reshape(x,1,[]); 
			y   = reshape(y,1,[]);
			std = reshape(std,1,[]);
			upper = y + diff(ylim)*std/2 * max_thickness;
			lower = y - diff(ylim)*std/2 * max_thickness;
			
			pl(1) = plot(ax, x,y);		% first plot the mean
			if isempty(p.Results.subshades)
				pl(2) = fill([x,x(end:-1:1)], [upper,lower(end:-1:1)], pl(1).Color*fill_alpha + (1-fill_alpha)*[1,1,1]);
				uistack(pl(2), 'bottom');
				set(pl(2),'EdgeColor','none')
				TDVPData.decrColorOrderIndex(ax);
			else
				std = p.Results.subshades';			% n x t
				for ii = 1:size(std,1)
					upper = lower + diff(ylim)*std(ii,:) * max_thickness;
					if max(std(ii,:)) > 5e-3
						% Only plot if max patch height > 1%
						pl(ii+1) = plot([x,x(end:-1:1)], [upper,lower(end:-1:1)], 'Color', colOrder(ii,:)*fill_alpha + (1-fill_alpha)*[1,1,1]);
						uistack(pl(ii+1),'bottom');
						TDVPData.decrColorOrderIndex(ax);
					end
					lower = upper;
				end
			end
		end
		
		function [m,n] = bestGrid(nPlots)
			%% [m,n] = bestGrid(nPlots)
			% find the best grid m x n to display n plots
			% require n>=m
			nAx = nPlots;
			
			% Simple stupid approach
			if nPlots == 1
				m = 1; n = 1;
				return
			elseif nPlots == 2
				m = 1; n = 2;
				return
			elseif nPlots <= 4
				m = 2; n = 2;
				return
			elseif nPlots <= 6
				m = 2; n = 3;
				return
			elseif nPlots <= 8
				m = 2; n = 4;
				return
			elseif nPlots <= 9
				m = 3; n = 3;
			elseif nPlots <= 12
				m = 3; n = 4; 
			end
		end
		
		function out = D2Str(treeMPS)
			%% function out = D2Str(treeMPS)
			%	generates string output of bond dimensions returned in a flattened cell array for each line to display
			if treeMPS.height == 0
				out = {num2str(treeMPS.D,'%-5u')};		% converts a chain to string; pack into cell
				return;
			end
			% else: hasChild!
			childStr = {};
			for i = 1:treeMPS.degree
				childStr{i} = TDVPData.D2Str(treeMPS.child(i));
			end

			for i = 1:treeMPS.degree
				% str1 for the first line of each child
				if i == 1
					str1 = sprintf('%-5u%s',treeMPS.D(1));
				elseif i == treeMPS.degree
					str1 = sprintf('  \\- ');
				else
					str1 = sprintf('  |- ');
				end
				if iscell(childStr{i})
					childStr{i}{1} = sprintf('%s%s',str1,childStr{i}{1});
					for j = 2:length(childStr{i})
						if i == treeMPS.degree
							childStr{i}{j} = sprintf('     %s',childStr{i}{j});
						else
							childStr{i}{j} = sprintf('  |  %s',childStr{i}{j});
						end
					end
				end
			end
			childStr = [childStr{:}];		% flatten cell array
			if treeMPS.isRoot
				out = sprintf('%s\n',childStr{:});
			else
				out = childStr;
			end
			
		end
		
		function out = dk2Str(treeMPS)
			%% function out = dk2Str(treeMPS)
			%	generates string output of bond dimensions returned in a flattened cell array for each line to display
			if treeMPS.height == 0
				out = {num2str(treeMPS.dk,'%-5u')};		% converts a chain to string; pack into cell
				return;
			end
			% else: hasChild!
			childStr = {};
			for i = 1:treeMPS.degree
				childStr{i} = TDVPData.dk2Str(treeMPS.child(i));
			end

			for i = 1:treeMPS.degree
				% str1 for the first line of each child
				if i == 1
					str1 = sprintf('%-5u%s',0);
				elseif i == treeMPS.degree
					str1 = sprintf('  \\- ');
				else
					str1 = sprintf('  |- ');
				end
				if iscell(childStr{i})
					childStr{i}{1} = sprintf('%s%s',str1,childStr{i}{1});
					for j = 2:length(childStr{i})
						if i == treeMPS.degree
							childStr{i}{j} = sprintf('     %s',childStr{i}{j});
						else
							childStr{i}{j} = sprintf('  |  %s',childStr{i}{j});
						end
					end
				end
			end
			childStr = [childStr{:}];		% flatten cell array
			if treeMPS.isRoot
				out = sprintf('%s\n',childStr{:});
			else
				out = childStr;
			end
			
		end
		
		function out = dkUsage(treeMPS)
			%% function out = dkUsage(treeMPS)
			%	generates string output of bond dimensions returned in a flattened cell array for each line to display
			thresh = 1e-1;								% minimum amplitude to be relevant
			if treeMPS.height == 0
				%% iterate through chain and estimate relevant dk usage
				dk = zeros(1,treeMPS.L);
				for ii = 1:treeMPS.L
					est = sqrt(sum(abs(treeMPS.Vmat{ii}).^2,2));				% measure importance by row-norm in dk! (column norm is always 1)
					est = find(est<thresh,1,'first');
					if isempty(est)
						dk(ii) = treeMPS.dk(ii);
					else
						dk(ii) = est;
					end
				end
				out = {num2str(dk,'%-5u')};		% converts a chain to string; pack into cell
				return;
			end
			% else: hasChild!
			childStr = {};
			for i = 1:treeMPS.degree
				childStr{i} = TDVPData.dkUsage(treeMPS.child(i));
			end

			for i = 1:treeMPS.degree
				% str1 for the first line of each child
				if i == 1
					str1 = sprintf('%-5u%s',0);
				elseif i == treeMPS.degree
					str1 = sprintf('  \\- ');
				else
					str1 = sprintf('  |- ');
				end
				if iscell(childStr{i})
					childStr{i}{1} = sprintf('%s%s',str1,childStr{i}{1});
					for j = 2:length(childStr{i})
						if i == treeMPS.degree
							childStr{i}{j} = sprintf('     %s',childStr{i}{j});
						else
							childStr{i}{j} = sprintf('  |  %s',childStr{i}{j});
						end
					end
				end
			end
			childStr = [childStr{:}];		% flatten cell array
			if treeMPS.isRoot
				out = sprintf('%s\n',childStr{:});
			else
				out = childStr;
			end
			
		end
		
		function out = subsample(data,fromDt, toDt)
			%% function out = subsample(data,fromDt, toDt)
			%	subsample data taken at intervals fromDt down to an interval toDt
			%
			%	for now only 2D-array supported
			if mod(toDt/fromDt,1) == 0 || mod(toDt/fromDt,1) > 0.95
				step = round(toDt/fromDt);
				if step == 1
					% no subsampling needed
					out = data;
				else
					out = data(1:step:end,:);
				end
			else
				warning('Incommensurate subsampling performed by interpolation');
				out = 0;	%% TODO: has to be implemented; use interp?
			end
		end
		
		function res = getTensorVNE(A,n)
			%% function res = getTensorVNE(A,n)
			%	calculate the von Neumann entropy for any possible partition of tensor A
			%	as a proxy for the entanglement and thus most optimal partition into subtensors.
			%   n gives highest order of bunched legs
			
			N = ndims(A);
			dA = size(A);
			res = {};
			k = 1;
			if nargin == 1
				n = 4;			% default setting -> up to quarts
			end

			% First look at single partitions
			% the ones with smallest vNE could be singled out alone
			for ii = 1:N
				[B, dB] = tensShape(A,'unfoldiso', ii,dA);
				[~,S,~] = svd2(B);
				res{k,1} = ii;
				res{k,2} = vonNeumannEntropy(S);
				res{k,3} = diag(S);
				fprintf('dim %d: %g\n',res{k,1:2});
				k = k+1;
			end
			
			if n == 1
				return;
			end
			
			% Second look at pairs
			for ii = 1:N
				for jj = ii+1:N
					[B, dB] = tensShape(A,'unfoldiso', [ii,jj],dA);
					[~,S,~] = svd2(B);
					res{k,1} = [ii,jj];
					res{k,2} = vonNeumannEntropy(S);
					res{k,3} = diag(S);
					pairSum = sum([res{res{k,1},2}]);
					pairDiff = abs(diff([res{res{k,1},2}]));
					fprintf('dim [%d,%d]: %5.3g;\t diff: %5.2g; \t sum: %5.3g; \t rel.Compression: %5.1f\n',res{k,1:2}, pairDiff, pairSum, (1-res{k,2}/pairSum)*100);
					k = k+1;
				end
			end
			
			if n == 2
				return;
			end
			
			% Third look into triads
			for ii = 1:N
				for jj = ii+1:N
					for kk = jj+1:N
						[B, dB] = tensShape(A,'unfoldiso', [ii,jj,kk],dA);
						[~,S,~] = svd2(B);
						res{k,1} = [ii,jj,kk];
						res{k,2} = vonNeumannEntropy(S);
						res{k,3} = diag(S);
						triadSum = sum([res{res{k,1},2}]);
						fprintf('dim [%d,%d,%d]: %5.3g;\t individual sum: %5.3g; \t rel.Compression: %5.1f\n',res{k,1:2}, triadSum, (1-res{k,2}/triadSum)*100);
						k = k+1;
					end
				end
			end

			if n == 3
				return;
			end
			
			% Fourth look into quarts
			for ii = 1:N
				for jj = ii+1:N
					for kk = jj+1:N
						for ll = kk+1:N
							[B, dB] = tensShape(A,'unfoldiso', [ii,jj,kk,ll],dA);
							[~,S,~] = svd2(B);
							res{k,1} = [ii,jj,kk,ll];
							res{k,2} = vonNeumannEntropy(S);
							res{k,3} = diag(S);
							triadSum = sum([res{res{k,1},2}]);
							fprintf('dim [%d,%d,%d,%d]: %5.3g;\t individual sum: %5.3g; \t rel.Compression: %5.1f\n',res{k,1:2}, triadSum, (1-res{k,2}/triadSum)*100);
							k = k+1;
						end
					end
				end
			end
			if n == 4
				return;
			end
		end
		
		function res = getTensorVNEAdjacency(A)
			%% function res = getTensorVNEAdjacency(A)
			%	calculate Adjacency matrix of the von Neumann entropy for any possible partition of tensor A
			
			N = ndims(A);
			out = TDVPData.getTensorVNE(A,2);
			
			res = zeros(N);						% N x N adjacency matrix
			
			for kk = N+1:size(out,1)
				ii = out{kk,1}(1); jj = out{kk,1}(2);
				pairSum = sum([out{out{kk,1},2}]);
				res(ii,jj) = (pairSum - out{kk,2})/2;
% 				res(ii,jj) = (pairSum - out{kk,2})/out{kk,2};	% weighted by residual entropy. 
			end
			return;
		end
		
		function res = getTensorRDMEntanglementAdjacency(A)
			%% function res = getTensorRDMEntanglementAdjacency(A)
			%	computes the entanglement between pairs of legs of A from the RDM of the legs
			%	created 11/08/2017 by FS
			
			N = ndims(A);
			dA = size(A);
			res = zeros(N,N);
			entropies = zeros(N,N,3);
			
			%% Contract into RDM of ii-jj to compute entanglement
			for jj = 2:N
				for ii = 1:(jj-1)
					ind1 = 1:N;
					ind1([ii,jj]) = [];
					rhoAB = contracttensors(conj(A),N,ind1, A,N,ind1);
					rhoA  = contracttensors(conj(A),N,[ind1,jj], A,N,[ind1,jj]);	% ii subsystem
					rhoB  = contracttensors(conj(A),N,[ind1,ii], A,N,[ind1,ii]);	% jj subsystem
					rhoAB = reshape(rhoAB,size(rhoAB,1)*size(rhoAB,2),[]);
					% quantum purity:
% 					entropies(ii,jj,1) = sum(diag(rhoA*rhoA));
% 					entropies(ii,jj,2) = sum(diag(rhoB*rhoB));
% 					entropies(ii,jj,3) = sum(diag(rhoAB*rhoAB));
					% entanglement entropy
					entropies(ii,jj,1) = -sum(diag( rhoA*log( rhoA)));
					entropies(ii,jj,2) = -sum(diag( rhoB*log( rhoB)));
					D = eig(rhoAB);											% entropy is only correct from eigenvalues
					entropies(ii,jj,3) = -D.'*log(D);						% basically same as getTensorVNE
				end
			end
			res = real(entropies);
		end
	end
end