/*
	STPlot v2.2
	API for plotting with gnuplot.

	Created by Stilianos Louca on 2008.11.23
	Last modified on 2012.06.07
	
	Tested using GCC 4.2.1 and gnuplot 4.4.2.
	
	License agreement:
	You are free to use and modify STPlot for personal or educational purposes as you like. 
	You may not redistribute it, modified or not, without giving credits to the original author. 
	The author would furthermore be grateful if he is cited in any scientific works making use of STPlot.
	
	See the 'User Guide' that came with this source for a thorough documentation of STPlot.
	
 */
 
#ifndef STPLOT_DEF
#define STPLOT_DEF
#include "STPlot.h"

namespace ST{


string STPLOT_GNUPLOT_PATH; 



#pragma mark - 
#pragma mark GnuplotPipe
#pragma mark - 


//copy constructor
GnuplotPipe::GnuplotPipe(const GnuplotPipe &original): 
STPipe(original) // copy all pipe-related parts
{
	// adopt the original's settings
	gnuplotPath 	= original.gnuplotPath;
	scriptPath 		= original.scriptPath;
	_foundGnuplot	= original._foundGnuplot;
	if(original.scriptStream!=NULL){
		scriptStream = fopen(scriptPath.c_str(), "a");
	}else{
		scriptStream = NULL;
	}
}



//assignment operator
GnuplotPipe &GnuplotPipe::operator=(const GnuplotPipe &original){
	STPipe::operator=(original); // assign all pipe-related parts

	// close any existing output streams
	if(scriptStream!=NULL){
		fclose(scriptStream);
		remove(scriptPath.c_str());
		scriptStream = NULL;
	}

	// adopt the original's settings
	gnuplotPath 	= original.gnuplotPath;
	_foundGnuplot	= original._foundGnuplot;
	scriptPath 		= original.scriptPath;
	if(original.scriptStream!=NULL){
		scriptStream = fopen(scriptPath.c_str(), "a");
	}
	
	return *this;
}
	


void GnuplotPipe::open(){
	setRequirePipe(false);
	tryToFindGnuplot();
	STPipe::open(ST_DEFAULT_PERSIST_PLOT ? gnuplotPath + " -persist" : gnuplotPath); 
}

void GnuplotPipe::open(string _gnuplotPath){
	setRequirePipe(false);
	tryToFindGnuplot(_gnuplotPath);
	STPipe::open(ST_DEFAULT_PERSIST_PLOT ? gnuplotPath + " -persist" : gnuplotPath); 
}


bool GnuplotPipe::tryToFindGnuplot(const string &firstGuess){
	gnuplotPath = firstGuess;
	if((gnuplotPath=="") || (!STPlot_aux_fileExists(gnuplotPath))){
		gnuplotPath = STPLOT_GNUPLOT_PATH;
	}else{ return (_foundGnuplot=true); }
	
	if((gnuplotPath=="") || (!STPlot_aux_fileExists(gnuplotPath))){
		if(!getProgramPath(string("gnuplot")+ST_GNUPLOT_IDEAL_VERSION, gnuplotPath)){
			gnuplotPath = "";
		}
	}else{ return (_foundGnuplot=true); }
	
	if((gnuplotPath=="") || (!STPlot_aux_fileExists(gnuplotPath))){
		if(!getProgramPath("gnuplot", gnuplotPath)){
			gnuplotPath = "";
		}
	}else{ return (_foundGnuplot=true); }
	
	if((gnuplotPath=="") || (!STPlot_aux_fileExists(gnuplotPath))){
		gnuplotPath = ST_DEFAULT_GNUPLOT_PATH;
	}else{ return (_foundGnuplot=true); }
	
	if((gnuplotPath=="") || (!STPlot_aux_fileExists(gnuplotPath))){
		gnuplotPath = "gnuplot";
	}else{ return (_foundGnuplot=true); }
	
	if((gnuplotPath=="") || (!STPlot_aux_fileExists(gnuplotPath))){
		return (_foundGnuplot=false);
	}else{ return (_foundGnuplot=true); }
}


bool GnuplotPipe::write(const char *format, ...){
	va_list args, logArgumentList;
	bool result;
	va_start(args, format);
	if(isInScriptMode()){
		va_start(logArgumentList, format);
		result=(vfprintf(scriptStream, format, args) >= 0);
		if(logStream != NULL){
			if(logStream->good()){
				vstreamprintf(*logStream, format, logArgumentList);
			}
		}
		va_end(logArgumentList);
	}else{
		result = STPipe::write(format, args);
	}
	va_end(args);
	return result;
}



bool GnuplotPipe::write(const string &message){
	if(isInScriptMode()){
		bool result = (fprintf(scriptStream, message.c_str()) >= 0);
		if((logStream != NULL) && (logStream->good())) (*logStream) << message;
		return result;
	}else{
		return STPipe::write(message);
	}	
}



void GnuplotPipe::transformToScript(){
	if(system(NULL) == 0) return; //no shell access, so script can not be executed.
	scriptPath = STPlot_aux_getUniqueTempFilePath();
	//tmpnam(NULL); // deprecated
	scriptStream = fopen(scriptPath.c_str(), "w");
}


void GnuplotPipe::transformToPipe(){
	if(scriptStream != NULL){
		fclose(scriptStream);
		remove(scriptPath.c_str());
		scriptStream = NULL;
	}
}


bool GnuplotPipe::executeScript(){
	if(!isInScriptMode()) return false;
	if(!foundGnuplot()) return false;
	if(fflush(scriptStream) == EOF) return false;
	system((gnuplotPath+" '" + scriptPath + "'" + (SUPPRESS_GNUPLOT_ERRORS ? " 2>/dev/null" : "")).c_str());
	return true;
}

bool GnuplotPipe::isInScriptMode() const{
	return (scriptStream != NULL);
}


string GnuplotPipe::currentGnuplotVersion() const{
	string version;
	if(!ST::getCommandOutput(gnuplotPath + " --version", true, version)){
		return "unknown";
	}
	int p0 = strlen("gnuplot ");
	int p1 = version.find(" patchlevel");
	int p2 = p1 + strlen(" patchlevel");
	return version.substr(p0, p1-p0);
}





#pragma mark -
#pragma mark Escaping & Greek symbols
#pragma mark 

//Used internally for replacing greek symbols with appropriate gnuplot symbols
//For more details see: http://t16web.lanl.gov/Kawano/gnuplot/label-e.html
const char* Symbols[] = 	{		"$Alpha$",		"$Beta$",		"$Gamma$",		"$Delta$",		"$Epsilon$",		"$Zeta$",		"$Eta$",			"$Theta$",		"$Iota$",		"$Kappa$",		"$Lambda$",		"$Mu$",			"$Nu$",			"$Xi$",			"$Omicron$",		"$Pi$",			"$Rho$",			"$Sigma$",		"$Tau$",			"$Upsilon$",		"$Phi$",			"$Chi$",			"$Psi$",			"$Omega$",
									"$alpha$",		"$beta$",		"$gamma$",		"$delta$",		"$epsilon$",		"$zeta$",		"$eta$",			"$theta$",		"$iota$",		"$kappa$",		"$lambda$",		"$mu$",			"$nu$",			"$xi$",			"$omicron$",		"$pi$",			"$rho$",			"$sigma$",		"$tau$",			"$upsilon$",		"$phi$",			"$chi$",			"$psi$",			"$omega$",
									"$vartheta$",	"$varphi$", "$overline$", "$underline$", "$partial$", "$Sum$", "$Prod$", "$infty$", "$times$"};
const char* GnuplotSymbols[] =	{	"{/Symbol A}",	"{/Symbol B}",	"{/Symbol G}",	"{/Symbol D}",	"{/Symbol E}",		"{/Symbol Z}",	"{/Symbol H}",		"{/Symbol Q}",	"{/Symbol I}",	"{/Symbol K}",	"{/Symbol L}",	"{/Symbol M}",	"{/Symbol N}",	"{/Symbol X}",	"{/Symbol O}",		"{/Symbol P}",	"{/Symbol R}",		"{/Symbol S}",	"{/Symbol T}",		"{/Symbol U}",		"{/Symbol F}",		"{/Symbol C}",		"{/Symbol Y}",		"{/Symbol W}",
									"{/Symbol a}",	"{/Symbol b}",	"{/Symbol g}",	"{/Symbol d}",	"{/Symbol e}",		"{/Symbol z}",	"{/Symbol h}",		"{/Symbol q}",	"{/Symbol i}",	"{/Symbol k}",	"{/Symbol l}",	"{/Symbol m}",	"{/Symbol n}",	"{/Symbol x}",	"{/Symbol o}",		"{/Symbol p}",	"{/Symbol r}",		"{/Symbol s}",	"{/Symbol t}",		"{/Symbol u}",		"{/Symbol f}",		"{/Symbol c}",		"{/Symbol y}",		"{/Symbol w}",
									"{/Symbol \\112}", "{/Symbol \\152}", "{/Symbol \\140}", "{/Symbol \\137}", "{/Symbol \\266}", "{/Symbol \\345}", "{/Symbol \\325}", "{/Symbol \\245}", "{/Symbol \\264}"};
const unsigned int symbolCount = 57;



string replaceSymbols(string s){
	for(unsigned int i=0; i<symbolCount; ++i){
		STPlot_aux_replaceAll(s, Symbols[i], GnuplotSymbols[i]);
	}
	return s;
}

string STPlot_newLineEscapedVersion(const string &s){
	string newVersion(s);
	STPlot_aux_replaceAll(newVersion, "\n", "\\n");
	return newVersion;
}


#pragma mark -
#pragma mark Range
#pragma mark -


Range::Range(){
	unsetRangeMin();
	unsetRangeMax();
}


Range::Range(double rangeMin, double rangeMax){
	setRangeMin(rangeMin);
	setRangeMax(rangeMax);
}


string Range::gnuplotExpression() const{
	ostringstream range;
	
	//determine precision based on the range's position & span
	int e = 2; 
	if(!rangeMinAuto){
		e = max(e, int(1 + log(rangeMin)/log(10)));
	}
	if(!rangeMaxAuto){
		e = max(e, int(1 + log(rangeMax)/log(10)));
	}
	if((!rangeMinAuto) && (!rangeMaxAuto)){
		e += max(0, int( 4 - log(rangeMax - rangeMin)/log(10.0))); //at least 4 orders of magnitude preciser than the minimum needed to separate the range bounds.
		range << setprecision(e);
	}
	
	if(rangeMinAuto){
		range << "[*";
	}else{
		range << "[" << std::scientific<< rangeMin;
	}
	if(rangeMaxAuto){
		range << ":*]";
	}else{
		range << ":" << std::scientific << rangeMax << "]";
	}
	return range.str();	
}


void Range::scaleRangeByFactor(double factor){
	if(rangeMinAuto || rangeMaxAuto) return;
	rangeMax += (factor - 1.0) * (rangeMax - rangeMin) / 2.0;
	rangeMin -= (factor - 1.0) * (rangeMax - rangeMin) / 2.0;
}


void Range::mergeWithRange(Range range){
	if(rangeMinAuto && (!range.rangeMinAuto)) setRangeMin(range.rangeMin);
	if(rangeMaxAuto && (!range.rangeMaxAuto)) setRangeMax(range.rangeMax);
}




Range Range::getNicePlotRange(const vector<double> &values, long start, long end, bool zeroIsReference, bool logarithmic, double epsilon){
	double minValue, maxValue;
	if(end<start) return Range();
	minValue=maxValue=values[0];
	for(long ti=start; ti<=end; ++ti){ 
		maxValue = max(maxValue, values[ti]); 
		minValue = min(minValue, values[ti]); 
	}
	return getNicePlotRange(minValue, maxValue, zeroIsReference, logarithmic, epsilon);
}


Range Range::getNicePlotRange(double minValue, double maxValue, bool zeroIsReference, bool logarithmic, double epsilon){
	if(logarithmic && (maxValue>0)){
		Range range; 
		range.setRangeMax(maxValue*1.2);
		return range;
	}else{
		const double median = (maxValue + minValue)/2;
		epsilon = abs(epsilon);
		if(maxValue-minValue<=epsilon){
			if(epsilon==0){ return (minValue==0 ? Range(-1,1) : Range(0.9*minValue, 1.1*minValue)); }
			else{
				return (abs(maxValue)<epsilon ? Range(-epsilon,+epsilon) : Range(0.9*minValue,1.1*maxValue));
			}
		}else if(zeroIsReference){
			return Range((minValue<0 ? 1.2*minValue : 0), (maxValue>0 ? 1.2*maxValue : 0));			
		}else if(minValue<0 && maxValue>0){
			return Range(1.2*minValue, 1.2*maxValue);
		}else{
			return Range((minValue<0 ? median-(median-minValue)*1.4 : max(0.0,median-(median-minValue)*1.4)), (maxValue>0 ? median+(maxValue-median)*1.4 : min(0.0,median+(maxValue-median)*1.4)));
		}
	}
}




#pragma mark -
#pragma mark MultiRange
#pragma mark -


MultiRange::MultiRange(){
	dim = 2;
}

MultiRange::MultiRange(unsigned int _dim){
	dim = min(3u, max(2u,_dim));
}



MultiRange::MultiRange(Range _xRange, Range _yRange){
	dim = 2;
	xRange = _xRange;
	yRange = _yRange;
}


MultiRange::MultiRange(Range _xRange, Range _yRange, Range _zRange){
	dim = 3;
	xRange = _xRange;
	yRange = _yRange;
	zRange = _zRange;
}



MultiRange::MultiRange(double _xRangeMin, double _xRangeMax, double _yRangeMin, double _yRangeMax){
	dim = 2;
	xRange = Range(_xRangeMin, _xRangeMax);
	yRange = Range(_yRangeMin, _yRangeMax);
}

MultiRange::MultiRange(double _xRangeMin, double _xRangeMax, double _yRangeMin, double _yRangeMax, double _zRangeMin, double _zRangeMax){
	dim = 3;
	xRange = Range(_xRangeMin, _xRangeMax);
	yRange = Range(_yRangeMin, _yRangeMax);
	zRange = Range(_zRangeMin, _zRangeMax);
}
	
string MultiRange::gnuplotExpression() const{
	switch(dim){
	case 1: return xRange.gnuplotExpression();
	case 2: return xRange.gnuplotExpression() + " " + yRange.gnuplotExpression();
	case 3: return xRange.gnuplotExpression() + " " + yRange.gnuplotExpression() + " " + zRange.gnuplotExpression();
	default:return "";
	}		
}


string MultiRange::gnuplotCommand() const{
	switch(dim){
	case 1: return "set xrange " + xRange.gnuplotExpression() + ";";
	case 2: return "set xrange " + xRange.gnuplotExpression() + ";\nset yrange " + yRange.gnuplotExpression() + ";";
	case 3: return "set xrange " + xRange.gnuplotExpression() + ";\nset yrange " + yRange.gnuplotExpression() + ";\nset zrange " + zRange.gnuplotExpression() + ";";
	default:return "";
	}
}


void MultiRange::scaleRangeByFactor(double factor){
	xRange.scaleRangeByFactor(factor);
	yRange.scaleRangeByFactor(factor);
	zRange.scaleRangeByFactor(factor);
}


void MultiRange::mergeWithRange(MultiRange range){
	xRange.mergeWithRange(range.xRange);
	yRange.mergeWithRange(range.yRange);
	zRange.mergeWithRange(range.zRange);
}


#pragma mark -
#pragma mark MultiLabel
#pragma mark -

MultiLabel::MultiLabel(){
	_dim = 2;
	_xLabel = "x";
	_yLabel = "y";
}

MultiLabel::MultiLabel(string xLabel, string yLabel){
	_dim = 2;
	_xLabel = xLabel;
	_yLabel = yLabel;
}

MultiLabel::MultiLabel(string xLabel, string yLabel, string zLabel){
	_dim = 3;
	_xLabel = xLabel;
	_yLabel = yLabel;
	_zLabel = zLabel;
}
	
string MultiLabel::gnuplotCommand(bool interpretSymbols, bool escapeUnderscores) const{
	string formattedXLabel = STPlot_newLineEscapedVersion((interpretSymbols ? replaceSymbols(_xLabel) : _xLabel));
	string formattedYLabel = STPlot_newLineEscapedVersion((interpretSymbols ? replaceSymbols(_yLabel) : _yLabel));
	string formattedZLabel = STPlot_newLineEscapedVersion((interpretSymbols ? replaceSymbols(_zLabel) : _zLabel));
	if(escapeUnderscores){
		STPlot_aux_replaceAll(formattedXLabel, "_", "\\\\_"); STPlot_aux_replaceAll(formattedYLabel, "_", "\\\\_"); STPlot_aux_replaceAll(formattedZLabel, "_", "\\\\_");
	}
	switch(_dim){
	case 2:
		return "set xlabel \"" + formattedXLabel + "\";\nset ylabel \"" + formattedYLabel + "\"" + (_yLabel.size() > ST_DEFAULT_ZLABEL_MAX_SIZE_BEFORE_ROTATION ? " rotate by 90;" : ";");
	default:
		return "set xlabel \"" + formattedXLabel + "\";\nset ylabel \"" + formattedYLabel + "\";\nset zlabel \"" + formattedZLabel + "\"" + (_zLabel.size() > ST_DEFAULT_ZLABEL_MAX_SIZE_BEFORE_ROTATION ? " rotate by 90;" : ";"); 
	}
}




#pragma mark -
#pragma mark Contours
#pragma mark -


Contours::Contours(bool enable){ 
	drawContours 	= enable;
	automatic 		= true;
	fixedCount		= false;
	incremental 	= false;
	placement 		= ST_DEFAULT_CONTOUR_PLACEMENT;
	individualContours	= ST_DEFAULT_INDIVIDUAL_CONTOURS;
}

Contours::Contours(const vector<double> &_levels){ 
	drawContours	= (! _levels.empty()); 
	automatic 		= false;
	incremental 	= false; 
	levels 			= _levels; 
	placement 		= ST_DEFAULT_CONTOUR_PLACEMENT;
	individualContours	= ST_DEFAULT_INDIVIDUAL_CONTOURS;
}

Contours::Contours(double singleLevel){ 
	drawContours	= true; 
	automatic 		= false; 
	incremental 	= false; 
	levels.push_back(singleLevel); 
	placement 		= ST_DEFAULT_CONTOUR_PLACEMENT;
	individualContours	= ST_DEFAULT_INDIVIDUAL_CONTOURS;
}

Contours::Contours(unsigned int _count){
	drawContours	= true;
	automatic		= true;
	incremental		= false;
	fixedCount		= true;
	count			= _count;
	placement		= ST_DEFAULT_CONTOUR_PLACEMENT;
	individualContours	= ST_DEFAULT_INDIVIDUAL_CONTOURS;
}

Contours::Contours(double _start, double _step, unsigned int _count){ 
	drawContours 			= true; 
	automatic 				= false; 
	incremental 			= true;  
	fixedCount				= true; 
	start 					= _start; 
	step 					= _step;  
	count					= _count;
	placement 				= ST_DEFAULT_CONTOUR_PLACEMENT;
	individualContours		= ST_DEFAULT_INDIVIDUAL_CONTOURS;
}

Contours::Contours(double _start, double _step, double _end){ 
	drawContours 			= true; 
	automatic 				= false; 
	incremental 			= true;  
	fixedCount 				= false; 
	start 					= _start; 
	step 					= _step; 
	end 					= _end;
	placement 				= ST_DEFAULT_CONTOUR_PLACEMENT;
	individualContours		= ST_DEFAULT_INDIVIDUAL_CONTOURS;
}



string Contours::contourPlacement2string(ContourPlacement placement){
	switch(placement){
	case ContourPlacementSurface:		return "surface";
	case ContourPlacementBottom:		return "base";
	case ContourPlacementBottomSurface:	return "both";
	default:							return "";
	}
}



string Contours::gnuplotCommand() const{
	ostringstream contours;
	if(! drawContours){
		return "unset contour;";
	}else{
		if(automatic && fixedCount){
			contours << "set contour " << contourPlacement2string(placement) << "; set cntrparam levels " << count << ";";
		}else if(automatic && (!fixedCount)){
			contours << "set contour " << contourPlacement2string(placement) << "; set cntrparam levels auto;";
		}else if((!automatic) && incremental && fixedCount){
			contours << "set contour " + contourPlacement2string(placement) + "; set cntrparam levels " << count << "; set cntrparam levels incremental " << start << ", " << step << ";";
		}else if((!automatic) && incremental && (!fixedCount)){
			contours << "set contour " + contourPlacement2string(placement) + "; set cntrparam levels incremental " << start << ", " << step << ", " << end << ";";
		}else if((!automatic) && (!incremental)){
			contours << "set contour " + contourPlacement2string(placement) + "; set cntrparam levels discrete " << levels[0]; 
			for(unsigned int i=1; i<levels.size(); ++i){
				contours << ", " << levels[i];
			}	
			contours << ";";
		}
	}
	contours << (individualContours ? " set clabel;" : " unset clabel;");
	return contours.str(); 
}



#pragma mark -
#pragma mark Tics
#pragma mark -


Tics::Tics(){
	automatic 	= true;
	startSet 	= false;
	stepSet		= false;
	endSet		= false;
	rotation	= 0;
}


Tics::Tics(unsigned int ticCount,...){
	automatic 	= false;
	rotation	= 0;
	va_list values;
	va_start(values, ticCount);
	add(ticCount, values);
	va_end(values);
	
}

void Tics::add(unsigned int ticCount, va_list values){
	for(unsigned int i=0; i<ticCount; ++i){
		customTicLabels.push_back("");
		customTicValues.push_back(va_arg(values, double));
		customTicWithCustomLabel.push_back(false);
	}
}

void Tics::add(const string &label, double value){
	customTicLabels.push_back(label);
	customTicValues.push_back(value);
	customTicWithCustomLabel.push_back(true);
}

void Tics::add(double value){
	customTicLabels.push_back("");
	customTicValues.push_back(value);
	customTicWithCustomLabel.push_back(false);
}

void Tics::add(unsigned int count, ...){
	va_list values;
	va_start(values, count);
	add(count, values);
	va_end(values);
}


void Tics::clearCustomTics(){
	customTicLabels.clear();
	customTicValues.clear();
	customTicWithCustomLabel.clear();
}

void Tics::disableAutomaticTics(){
	automatic = false;
}

void Tics::enableAutomaticTics(){
	automatic 		= true;
	startSet 		= false;
	stepSet			= false;
	endSet			= false;
}

void Tics::enableAutomaticTics(double _step){
	automatic 		= true;
	startSet 		= false;
	stepSet			= true;
	endSet			= false;
	step			= _step;
}


void Tics::enableAutomaticTics(double _start, double _step){
	automatic 		= true;
	startSet 		= true;
	stepSet			= true;
	endSet			= false;
	start			= _start;
	step			= _step;
}


void Tics::enableAutomaticTics(double _start, double _step, double _end){
	automatic 		= true;
	start 			= true;
	stepSet			= true;
	endSet			= true;
	start			= _start;
	step			= _step;
	end				= _end;
}


void Tics::setToRegularIntegerValues(double rangeMin, double rangeMax){
	disableAutomaticTics();
	clearCustomTics();
	vector<long> grid;
	STPlot_getRegularIntegerGrid(ceil(rangeMin), floor(rangeMax), grid);
	for(int j=0; j<grid.size(); ++j){
		add(double(grid[j]));
	}
}



string Tics::gnuplotCommand(string ticsName, bool interpretSymbols, bool escapeUnderscores) const{
	return gnuplotCommand(ticsName, interpretSymbols, escapeUnderscores, 1.0);
}

	
string Tics::gnuplotCommand(string ticsName, bool interpretSymbols, bool escapeUnderscores, double rescale) const{
	ostringstream expression;
	string ticname;
	if(automatic){
		if(!stepSet){
			expression << "set " << ticsName << " autofreq";
		}else if(!startSet){
			expression << "set " << ticsName << " " << step;
		}else if(!endSet){
			expression << "set " << ticsName << " " << start << ", " << step;
		}else{
			expression << "set " << ticsName << " " << start << ", " << step << ", " << end;
		}
	}
	
	if((!automatic) && customTicValues.empty()){
		expression << "unset " << ticsName;
	}else if((!automatic) && (!customTicValues.empty())){
		expression << "set " << ticsName << " (";
		for(unsigned int i=0; i<customTicValues.size(); ++i){
			if(customTicWithCustomLabel[i]){
				ticname = STPlot_newLineEscapedVersion((interpretSymbols ? replaceSymbols(customTicLabels[i]) : customTicLabels[i]));
				if(escapeUnderscores) STPlot_aux_replaceAll(ticname, "_", "\\\\_");
				expression << (i == 0 ? "" : ", ") << "\"" << ticname << "\" " << customTicValues[i];
			}else{
				expression << (i == 0 ? "" : ", ") << customTicValues[i];
			}
		}
		expression << ")";
	}else if(automatic && (!customTicValues.empty())){
		expression << "\nset " << ticsName << " add (";
		for(unsigned int i=0; i<customTicValues.size(); ++i){
			if(customTicWithCustomLabel[i]){
				ticname = STPlot_newLineEscapedVersion((interpretSymbols ? replaceSymbols(customTicLabels[i]) : customTicLabels[i]));
				if(escapeUnderscores) STPlot_aux_replaceAll(ticname, "_", "\\\\_");
				expression << (i == 0 ? "" : ", ") << "\"" << ticname << "\" " << customTicValues[i];
			}else{
				expression << (i == 0 ? "" : ", ") << customTicValues[i];
			}
		}
		expression << ")";
	}
	
	expression << " scale " << rescale;
	
	// rotation and alignment
	const double rotationM = STPlot_aux_moduloInterval(rotation,-180,180);
	if(rotationM!=0){
		expression << " rotate by " << rotation << " ";
		if(ticsName=="xtics"){
			expression << (rotationM>=0 ? "right" : "left");
		}else if(ticsName=="x2tics"){
			expression << (rotationM>=0 ? "left" : "right");
		}else if(ticsName=="ytics"){
			expression << (rotationM>=-90 && rotationM<=90 ? "right" : "left");
		}else if(ticsName=="y2tics" || ticsName=="cbtics"){
			expression << (rotationM>=-90 && rotationM<=90 ? "left" : "right");
		}
	}
	expression << ";";
	return expression.str();
}



void Tics::setRotation(double angle){
	rotation = angle;
}




#pragma mark -
#pragma mark Label
#pragma mark

string Label::gnuplotCommand(long tag, bool interpretSymbols, int dim, const string &plotFont, bool escapeUnderscores) const{
	ostringstream command;
	string formattedText = STPlot_newLineEscapedVersion((interpretSymbols ? replaceSymbols(text) : text));
	if(escapeUnderscores) STPlot_aux_replaceAll(formattedText, "_", "\\\\_");
	command << "set label " << tag << " \"" << formattedText << "\" at" << (wrtCoordinates ? " first" : " graph") << " " << x;
	if(dim>1) command << ", " << y;
	if(dim>2) command << ", " << z;
	switch(alignment){
	case LabelAlignmentLeft: command << " left"; break;
	case LabelAlignmentRight: command << " right"; break;
	case LabelAlignmentCenter: command << " center"; break;
	}
	command << " front";
	if(fontSize>=0){
		command << " font '" << (plotFont=="" ? "*" : plotFont) << "," << fontSize << "'";
	}
	return command.str();
}


Label::Label(){ 
	x=y=z=0; wrtCoordinates=false; alignment = LabelAlignmentLeft; fontSize = -1; 
}


Label::Label(const string &_text, bool _wrtCoordinates, double _x, double _y, LabelAlignment _alignment){
	text 			= _text;
	wrtCoordinates 	= _wrtCoordinates;
	x 				= _x;
	y 				= _y;
	alignment 		= _alignment;
	fontSize 		= -1;
}



Label::Label(const string &_text, bool _wrtCoordinates, double _x, double _y, double _z, LabelAlignment _alignment){
	text 			= _text;
	wrtCoordinates 	= _wrtCoordinates;
	x 				= _x;
	y 				= _y;
	z 				= _z;
	alignment 		= _alignment;
	fontSize 		= -1;
}





#pragma mark -
#pragma mark Plot Source
#pragma mark -


PlotSource::PlotSource(){
	plotType		= PlotTypeNone;
	colorType		= PlotColorTypeAuto;
}



PlotSource::PlotSource(PlotType _plotType, string _title, string _filePath, unsigned int _index, ...){
	plotType	= _plotType;
	filePath 	= _filePath;
	index 		= _index;
	title		= _title;
	columns.resize(columnCountForPlotType(plotType));
	va_list argumentList;
	va_start(argumentList, _index);
	for(unsigned long i=0; i<columns.size(); ++i){
		columns[i] = va_arg(argumentList, unsigned int);	
	}
	va_end(argumentList);
	
	xStart 			= 0;
	yStart 			= 0;
	xStep			= 1;
	yStep			= 1;
	xCount			= 0;
	yCount			= 0;
	minX 			= 0;
	maxX	 		= 1;
	minY			= 0;
	maxY			= 1;
	lineWidth		= 1;
	lineWidthAuto 	= true;
	color.setRGB(1,0,0);
	colorType 		= PlotColorTypeAuto;
	lineType		= -1;
	lineTypeAuto	= true;
	pointSize		= 1;
	pointSizeAuto	= true;
	pointType		= 1;
	pointTypeAuto	= true;
	vectorHeadSize 	= ST_DEFAULT_VECTOR_HEAD_SIZE;
	vectorHeadStyle	= ST_DEFAULT_VECTOR_HEAD_STYLE;	
	filledCurveStyle= FilledCurveStyleAbove;
	fillReference 	= 0;	
}





PlotSource::~PlotSource(){
	// nothing to be done for now
}


/*
//copy constructor
PlotSource::PlotSource(const PlotSource &original){
	plotType			= original.plotType;
	filePath	 		= original.filePath;
	index 				= original.index;
	title 				= original.title;
	xStart				= original.xStart;
	yStart				= original.yStart;
	xStep				= original.xStep;
	yStep				= original.yStep;
	xCount				= original.xCount; 
	yCount				= original.yCount;
	minX				= original.minX;
	minY				= original.minY;
	maxX				= original.maxX;
	maxY				= original.maxY;
	lineWidth			= original.lineWidth;
	lineWidthAuto		= original.lineWidthAuto;
	color				= original.color;
	colorValue			= original.colorValue;
	colorType			= original.colorType;
	lineType			= original.lineType;
	lineTypeAuto		= original.lineTypeAuto;
	pointSize			= original.pointSize;
	pointSizeAuto		= original.pointSizeAuto;
	pointType			= original.pointType;
	pointTypeAuto		= original.pointTypeAuto;
	vectorHeadSize		= original.vectorHeadSize;
	vectorHeadStyle		= original.vectorHeadStyle;
	xTransformationFunction = original.xTransformationFunction;
	yTransformationFunction = original.yTransformationFunction;
	zTransformationFunction = original.zTransformationFunction;
	wTransformationFunction = original.wTransformationFunction;
	filledCurveStyle	= original.filledCurveStyle;
	fillReference 		= original.fillReference;
	columns				= original.columns;

}


//assignment operator
PlotSource &PlotSource::operator=(const PlotSource &original){
	plotType			= original.plotType;
	filePath	 		= original.filePath;
	index 				= original.index;
	title 				= original.title;
	xStart				= original.xStart;
	yStart				= original.yStart;
	xStep				= original.xStep;
	yStep				= original.yStep;
	xCount				= original.xCount; 
	yCount				= original.yCount;
	minX				= original.minX;
	minY				= original.minY;
	maxX				= original.maxX;
	maxY				= original.maxY;
	lineWidth			= original.lineWidth;
	lineWidthAuto		= original.lineWidthAuto;
	color				= original.color;
	colorValue			= original.colorValue;
	colorType			= original.colorType;
	lineType			= original.lineType;
	lineTypeAuto		= original.lineTypeAuto;
	pointSize			= original.pointSize;
	pointSizeAuto		= original.pointSizeAuto;
	pointType			= original.pointType;
	pointTypeAuto		= original.pointTypeAuto;
	vectorHeadSize		= original.vectorHeadSize;
	vectorHeadStyle		= original.vectorHeadStyle;
	xTransformationFunction = original.xTransformationFunction;
	yTransformationFunction = original.yTransformationFunction;
	zTransformationFunction = original.zTransformationFunction;
	wTransformationFunction = original.wTransformationFunction;
	filledCurveStyle	= original.filledCurveStyle;
	fillReference 		= original.fillReference;
	columns				= original.columns;

	return *this;
}
*/

unsigned int PlotSource::columnCountForPlotType(PlotType type){
	switch(type){
	case PlotTypeNone: 
		return 0;
	case PlotType1DPoints:
		return 1;
	case PlotType1DColoredPoints:
		return 4;
	case PlotType1DVectors:
		return 2;
	case PlotType1DColoredVectors:
		return 5;
	case PlotType2DPoints:
		return 2;
	case PlotType2DColoredPoints:
		return 5;
	case PlotType2DVectors:
		return 4;
	case PlotType2DColoredVectors:
		return 7;
	case PlotType2DVariableSizedSpheres:
		return 3;
	case PlotType2DColoredVariableSizedSpheres:
		return 6;
	case PlotType2DCurve:
		return 2;
	case PlotType2DCurvePoints:
		return 2;
	case PlotType2DGraphLinear:
		return 2;
	case PlotType2DGraphLinearPoints:
		return 2;
	case PlotType2DGraphSplines:
		return 2;
	case PlotType2DGraphSplinesPoints:
		return 2;
	case PlotType2DFilledCurve:
		return 2;
	case PlotType2DFilledCurves:
		return 3;
	case PlotType2DAnalytic:
		return 0;
	case PlotType3DPoints:
		return 3;
	case PlotType3DColoredPoints:
		return 6;
	case PlotType3DVectors:
		return 6;
	case PlotType3DColoredVectors:
		return 9;
	case PlotType3DCurve:
		return 3;
	case PlotType3DCurvePoints:
		return 3;
	case PlotType1DHistogram:
		return 2;
	case PlotType1DEmpiricalDensity:
		return 2;
	case PlotType2DScatteredValueMap:
		return 3;
	case PlotType2DGridValueMap:
		return 0;
	case PlotType2DScannedValueMap:
		return 3;
	case PlotType2DColorImage:
		return 0;
	case PlotType2DValueImage:
		return 0;
	case PlotType3DScatteredSurface:
		return 3;
	case PlotType3DGridSurface:
		return 0;
	case PlotType3DScannedSurface:
		return 3;
	case PlotType3DValuedScannedSurface:
		return 4;
	default:
		return 0;
	}
}



void PlotSource::setLineWidth(double width){
	lineWidthAuto 	= false;
	lineWidth		= max(0.0, width);
}


void PlotSource::setLineWidthAuto(){
	lineWidthAuto 	= true;
}


void PlotSource::setColor(const Color &c){
	colorType	= PlotColorTypeAsColor;
	color		= c;
}


void PlotSource::setColorAsValue(double value){
	colorType	= PlotColorTypeAsValue;
	colorValue	= value;
}


void PlotSource::setColorAuto(){
	colorType	= PlotColorTypeAuto;
}



void PlotSource::setLineType(int type){
	lineTypeAuto	= false;
	lineType		= type;
}



void PlotSource::setLineTypeAuto(){
	lineTypeAuto = true;
}



void PlotSource::setPointSize(double size){
	pointSizeAuto	= false;
	pointSize		= size;
}



void PlotSource::setPointSizeAuto(){
	pointSizeAuto = true;
}


void PlotSource::setPointType(int type){
	pointTypeAuto	= false;
	pointType		= type;
}



void PlotSource::setPointTypeAuto(){
	pointTypeAuto = true;
}


void PlotSource::setColorLineTypeLineWidth(const Color &c, int type, double width){
	setColor(c); setLineType(type); setLineWidth(width);
}



string PlotSource::gnuplotLineWidthExpression() const{	
	ostringstream expression;
	if(!lineWidthAuto){
		expression << " linewidth " << lineWidth;
	}
	return expression.str();
}




string PlotSource::gnuplotColorExpression() const{
	switch(plotType){
	case PlotType1DColoredPoints:
	case PlotType1DColoredVectors:
	case PlotType2DColoredPoints:
	case PlotType2DColoredVectors:
	case PlotType2DColoredVariableSizedSpheres:
	case PlotType3DColoredPoints:
	case PlotType3DColoredVectors:
		return "lc rgb variable";
	
	default:
		ostringstream expression;
		if(colorType==PlotColorTypeAsColor){
			expression << "lc rgb '#" << color.getHex() << "'";
		}else if(colorType==PlotColorTypeAsValue){
			expression << "lc palette cb " << colorValue;		
		}
		return expression.str(); 
	}
}



/* // DEPRECATED, WRONG WAY
string PlotSource::gnuplotArtificialColorColumnExpression() const{
	if(colorType==PlotColorTypeAsValue){
		return ":("+STPlot_makeString(colorValue)+")";
	}else{
		return "";
	}
}
*/


string PlotSource::gnuplotLineTypeExpression() const{
	ostringstream expression;
	if(!lineTypeAuto){
		expression << "linetype " << lineType;
	}
	return expression.str();
}



string PlotSource::gnuplotPointSizeExpression() const{
	ostringstream expression;
	if(!pointSizeAuto){
		expression << "pointsize " << pointSize;
	}
	return expression.str();
}


string PlotSource::gnuplotPointTypeExpression() const{
	ostringstream expression;
	if(!pointTypeAuto){
		expression << "pointtype " << pointType;
	}
	return expression.str();
}


string PlotSource::gnuplotFilledCurveStyleExpression() const{
	switch(filledCurveStyle){
	case FilledCurveStyleAbove: return "above";
	case FilledCurveStyleBelow: return "below";
	default:					return "";
	}
}

	
unsigned int PlotSource::plotDimension() const{
	switch(plotType){
	case PlotTypeNone:
		return 0;
	case PlotType1DPoints:
	case PlotType1DColoredPoints:
	case PlotType1DVectors:
	case PlotType1DColoredVectors:
	case PlotType2DPoints:
	case PlotType2DColoredPoints:
	case PlotType2DVectors:
	case PlotType2DColoredVectors:
	case PlotType2DCurve:
	case PlotType2DCurvePoints:
	case PlotType2DGraphLinear:
	case PlotType2DGraphLinearPoints:
	case PlotType2DGraphSplines:
	case PlotType2DGraphSplinesPoints:
	case PlotType2DVariableSizedSpheres:
	case PlotType2DColoredVariableSizedSpheres:
	case PlotType2DFilledCurve:
	case PlotType2DFilledCurves:
	case PlotType2DAnalytic:
	case PlotType1DHistogram:
	case PlotType1DEmpiricalDensity:
	case PlotType2DScatteredValueMap:
	case PlotType2DGridValueMap:
	case PlotType2DScannedValueMap:
	case PlotType2DColorImage:
	case PlotType2DValueImage:
		return 2;
	case PlotType3DPoints:
	case PlotType3DColoredPoints:
	case PlotType3DVectors:
	case PlotType3DColoredVectors:
	case PlotType3DCurve:
	case PlotType3DCurvePoints:
	case PlotType3DScatteredSurface:
	case PlotType3DGridSurface:
	case PlotType3DScannedSurface:
	case PlotType3DValuedScannedSurface:
		return 3;
	default:
		return 2;
	}
}


bool PlotSource::usesSourceFile() const{
	return (plotType!=PlotType2DAnalytic);
}


bool PlotSource::isPlotWithEPSBlurProblems() const{
	return ((plotType==PlotType2DValueImage) || (plotType==PlotType2DColorImage) || (plotType==PlotType2DGridValueMap));
}


bool PlotSource::isPlotWhichShowsColorbox() const{
	switch(plotType){
	case PlotType2DScatteredValueMap:
	case PlotType2DGridValueMap:
	case PlotType2DScannedValueMap:
	case PlotType2DValueImage:
	case PlotType3DScatteredSurface:
	case PlotType3DGridSurface:
	case PlotType3DScannedSurface:
	case PlotType3DValuedScannedSurface:
		return true;
	default:
		return (colorType==PlotColorTypeAsValue);
	}

}


#pragma mark -
#pragma mark Plot
#pragma mark -


Plot::Plot()
{
	showKeys 			= false;
	showColorBox		= false;
	gridColoringType	= ST_DEFAULT_GRID_COLORING_TYPE;
	gridSmoothingType	= ST_DEFAULT_GRID_SMOOTHING_TYPE;
	showColorMap		= true;
	xDistanceScale 		= 1;
	yDistanceScale 		= 1;
	xGridSize			= ST_DEFAULT_ST_GNUPLOT_GRIDSIZE;
	yGridSize			= ST_DEFAULT_ST_GNUPLOT_GRIDSIZE;
	colorPalette		= ColorPaletteGrey;
	xLogarithmic		= yLogarithmic = zLogarithmic = false;
	backgroundColor		= Color::white();
}


Plot::Plot(const string &_title)
{
	title				= _title;
	showKeys 			= false;
	showColorBox		= false;
	gridColoringType	= ST_DEFAULT_GRID_COLORING_TYPE;
	gridSmoothingType	= ST_DEFAULT_GRID_SMOOTHING_TYPE;
	showColorMap		= true;
	xDistanceScale 		= 1;
	yDistanceScale 		= 1;
	xGridSize			= ST_DEFAULT_ST_GNUPLOT_GRIDSIZE;
	yGridSize			= ST_DEFAULT_ST_GNUPLOT_GRIDSIZE;
	colorPalette		= ColorPaletteGrey;
	xLogarithmic		= yLogarithmic = zLogarithmic = false;
	backgroundColor		= Color::white();
}


Plot::Plot(	const vector<PlotSource> 	&_Sources,
			const MultiLabel 			&_axisLabels, 
			const MultiRange 			&_Range, 
			const string 				&_Title)
{
	sources 	= _Sources;
	axisLabels	= _axisLabels;
	range		= _Range;
	title 		= _Title;
	showKeys 			= true;
	showColorBox		= false;
	gridColoringType	= ST_DEFAULT_GRID_COLORING_TYPE;
	gridSmoothingType	= ST_DEFAULT_GRID_SMOOTHING_TYPE;
	showColorMap		= true;
	xDistanceScale 		= 1;
	yDistanceScale 		= 1;
	xGridSize			= ST_DEFAULT_ST_GNUPLOT_GRIDSIZE;
	yGridSize			= ST_DEFAULT_ST_GNUPLOT_GRIDSIZE;
	colorPalette		= ColorPaletteGrey;
	xLogarithmic		= yLogarithmic = zLogarithmic = false;
	backgroundColor		= Color::white();
}

	
Plot::Plot(	const PlotSource 	&_Source, 
			const MultiLabel 	&_axisLabels, 
			const MultiRange 	&_Range, 
			const string 		&_Title)
{
	sources.push_back(_Source);
	axisLabels	= _axisLabels;
	range 		= _Range;
	title 		= _Title;
	showKeys 			= false;
	showColorBox		= false;
	gridColoringType	= ST_DEFAULT_GRID_COLORING_TYPE;
	gridSmoothingType	= ST_DEFAULT_GRID_SMOOTHING_TYPE;
	showColorMap		= true;
	xDistanceScale 		= 1;
	yDistanceScale 		= 1;
	xGridSize			= ST_DEFAULT_ST_GNUPLOT_GRIDSIZE;
	yGridSize			= ST_DEFAULT_ST_GNUPLOT_GRIDSIZE;
	colorPalette		= ColorPaletteGrey;
	xLogarithmic		= yLogarithmic = zLogarithmic = false;
	backgroundColor		= Color::white();
}



unsigned int Plot::compatibilityGroup(PlotType type){
	switch(type){
	case PlotTypeNone:
		return 0;
	case PlotType1DPoints:
	case PlotType1DColoredPoints:
	case PlotType1DVectors:
	case PlotType1DColoredVectors:
	case PlotType2DPoints:
	case PlotType2DColoredPoints:
	case PlotType2DVectors:
	case PlotType2DColoredVectors:
	case PlotType2DCurve:
	case PlotType2DCurvePoints:
	case PlotType2DGraphLinear:
	case PlotType2DGraphLinearPoints:
	case PlotType2DGraphSplines:
	case PlotType2DGraphSplinesPoints:
	case PlotType2DVariableSizedSpheres:
	case PlotType2DColoredVariableSizedSpheres:
	case PlotType2DFilledCurve:
	case PlotType2DFilledCurves:
	case PlotType2DAnalytic:
		return 1;
	case PlotType3DPoints:
	case PlotType3DColoredPoints:
	case PlotType3DVectors:
	case PlotType3DColoredVectors:
	case PlotType3DCurve:
	case PlotType3DCurvePoints:
		return 2;
	case PlotType1DHistogram:
		return 3;
	case PlotType1DEmpiricalDensity:
		return 4;
	case PlotType2DScatteredValueMap:
		return 5;
	case PlotType2DGridValueMap:
		return 6;
	case PlotType2DScannedValueMap:
		return 7;
	case PlotType2DColorImage:
	case PlotType2DValueImage:
		return 8;
	case PlotType3DScatteredSurface:
		return 9;
	case PlotType3DGridSurface:
	case PlotType3DScannedSurface:
	case PlotType3DValuedScannedSurface:
		return 10;
	default:
		return 11;
	}
}

bool Plot::plotTypesCompatible(PlotType type1, PlotType type2){
	return (compatibilityGroup(type1) == compatibilityGroup(type2));
}



string Plot::vectorHeadStyle2string(VectorHeadStyle style, double headSize){
	switch(style){
	case VectorHeadStyleFilledTriangle: 	return "head filled front size " + STPlot_makeString(headSize) + "*0.6,20";
	case VectorHeadStyleFilledArrow:		return "head filled front size " + STPlot_makeString(headSize) + "*0.6,20,30";
	case VectorHeadStyleEmptyTriangle:		return "head empty front size " + STPlot_makeString(headSize) + "*0.6,20";
	case VectorHeadStyleEmptyArrow:			return "head empty front size " + STPlot_makeString(headSize) + "*0.6,20,30";
	case VectorHeadStyleNone:				return "nohead";
	default: 								return "default";
	}
}



string Plot::evaluateFunctionForColumn(string trafo, unsigned long column){
	string s = "$" + STPlot_makeString(column);
	if(trafo == ""){
		return s;
	}else{
		STPlot_aux_replaceAll(trafo, "#", s.c_str());
		return trafo;
	}
}




string Plot::plotErrorDescription(PlotError error){
	switch(error){
	case PlotErrorNone: 			return "None";
	case PlotErrorInvalidSyntax: 	return "Invalid syntax";
	case PlotErrorException: 		return "Exception";
	case PlotErrorFileReadError:	return "File read error";
	case PlotErrorFileWriteError:	return "File write error";
	case PlotErrorFileNotFoundError:return "File not found";
	case PlotErrorPermissionError:	return "No permission";
	case PlotErrorEOF:				return "EOF";
	case PlotErrorPlotterNotFound:	return "Gnuplot not found";
	case PlotErrorUnknownError:		return "Unknown error";
	case PlotErrorNothingToBeDone:	return "Nothing to be done";
	case PlotErrorParameterError:	return "Parameter error";
	case PlotErrorIncompatiblePlots:return "Incompatiple plots";
	default: 						return "";
	}
}



string Plot::integrationMethod2string(IntegrationMethod method){
	switch(method){
	case IntegrationMethodForwardEuler:			return "Forward Euler";
	case IntegrationMethodClassicalRungeKutta:	return "Runge-Kutta";
	default: 									return "Unknown";
	}
}

unsigned int Plot::highestPlotDimension() const{
	unsigned int hpd = 0;
	for(unsigned int n=0; n<sources.size(); ++n){
		hpd = max(hpd, sources[n].plotDimension());
	}
	return hpd;
}


bool Plot::hasPlotWithEPSBlurProblems() const{
	for(long i=0; i<sources.size(); ++i){
		if(sources[i].isPlotWithEPSBlurProblems()) return true;
	}
	return false;
}



bool Plot::hasPlotWhichShowsColorbox() const{
	for(long i=0; i<sources.size(); ++i){
		if(sources[i].isPlotWhichShowsColorbox()) return true;
	}
	return false;
}


string	Plot::gnuplotDataSeparatorExpression() const{
	if(dataSeparator==""){
		return "set datafile separator whitespace;";
	}else{
		return "set datafile separator \""+dataSeparator+"\";";
	}
}




void Plot::setLineWidthForAllCurrentSources(double width){
	for(long s=0; s<sources.size(); ++s){ sources[s].setLineWidth(width); }
}

void Plot::setLineWidthAutoForAllCurrentSources(){
	for(long s=0; s<sources.size(); ++s){ sources[s].setLineWidthAuto(); }
}

void Plot::setLineTypeForAllCurrentSources(int type){
	for(long s=0; s<sources.size(); ++s){ sources[s].setLineType(type); }
}

void Plot::setLineTypeAutoForAllCurrentSources(){
	for(long s=0; s<sources.size(); ++s){ sources[s].setLineTypeAuto(); }
}

void Plot::setPointSizeForAllCurrentSources(double size){
	for(long s=0; s<sources.size(); ++s){ sources[s].setPointSize(size); }
}

void Plot::setPointSizeAutoForAllCurrentSources(){
	for(long s=0; s<sources.size(); ++s){ sources[s].setPointSizeAuto(); }
}

void Plot::setPointTypeForAllCurrentSources(int type){
	for(long s=0; s<sources.size(); ++s){ sources[s].setPointType(type); }
}

void Plot::setPointTypeAutoForAllCurrentSources(){
	for(long s=0; s<sources.size(); ++s){ sources[s].setPointTypeAuto(); }
}


void Plot::setNiceColorsForAllCurrentSources(){
	vector<Color> colors;
	generateDiscreteColors(sources.size(), true, colors, ColorBlindnessNone);
	for(long s=0; s<sources.size(); ++s){ sources[s].setColor(colors[s]); }
}




PlotError Plot::writeCommandToPipe(	GnuplotPipe 	&pipe, 
									double 			originX, 
									double 			originY,
									double 			availableWidth, 
									double 			availableHeight, 
									double 			resizeNonPlotAreaWidthByFactor,
									double 			resizeNonPlotAreaHeightByFactor,
									double			resizeTicsByFactor,
									const string	&windowFont,
									bool 			interpretSymbols,
									unsigned long 	startSource, 
									const string 	&tag,
									bool 			escapeUnderscoresInTitle,
									bool 			escapeUnderscoresInTicsAndLabels,
									string 			commandJustAfterFirstPlot,
									long 			&nextAvailLabelTag) const{
	PlotError 	error = PlotErrorNone;
	double		xStep, yStep;
	double 		rotationInducedScaleX, 	rotationInducedScaleY;
	double 		lmargin, rmargin, bmargin, tmargin;
	string 		sourceTitle;
	MultiRange	eRange = range; //effective range. Some adjustments might be needed to the original range.
	const int	dim = highestPlotDimension();
	string formattedTitle = STPlot_newLineEscapedVersion((interpretSymbols ? replaceSymbols(title) : title));
	if(escapeUnderscoresInTitle) formattedTitle = STPlot_aux_get_replaceAll(formattedTitle, "_", "\\\\_");

	//skip trivial plots
	for(;startSource<sources.size() && sources[startSource].plotType == PlotTypeNone; ++startSource);
	
	if(startSource >= sources.size()){
		//nothing to be plotted. Display empty box with title
		lmargin	= originX + availableWidth * (ST_DEFAULT_RELATIVE_LMARGIN);
		rmargin = originX + availableWidth * (1.0 - ST_DEFAULT_RELATIVE_RMARGIN);
		bmargin = originY + availableHeight * (ST_DEFAULT_RELATIVE_BMARGIN);
		tmargin	= originY + availableHeight * (1.0 - ST_DEFAULT_RELATIVE_TMARGIN);
		pipe.write("set lmargin at screen %e\nset rmargin at screen %e\n", lmargin, rmargin);
		pipe.write("set bmargin at screen %e\nset tmargin at screen %e\n", bmargin, tmargin);
		
		pipe.write("unset xtics;\nunset xlabel;\nunset ytics;\nunset ylabel;\nunset border;\nset xrange [1:2];\nset yrange [1:2];\n");
		pipe.write("set title \"%s%s\" offset 0, 0;\nunset key;\nunset colorbox;\nunset pm3d;\n", (tag=="" ? "" : ("("+tag+") ").c_str()), formattedTitle.c_str());
		pipe.write("plot 1/0 notitle\n");
		goto CLOSING; //go to end of this function to finish with some technical closing commands.
		return PlotErrorNone;
	}
	
	
	pipe.write("\n");
	if(startSource == 0){
		//color palette
		int redCode, greenCode, blueCode;
		colorPalette2PaletteCodes(colorPalette, redCode, greenCode, blueCode);
		pipe.write("set palette rgbformulae %i,%i,%i;\n", redCode, greenCode, blueCode);
		
		//color conversion formula needed for a bunch of PlotTypes.
		//(r,g,b) to int \in [0, 256^3 - 1], whereas (r,g,b) \in [0,1]^3
		pipe.write("%s(r,g,b) = int(r * 255) * 65536 + int(g * 255) * 256 + int(b * 255) * 1; %smapping (r,g,b) triplets to ints, as encoded by gnuplot\n", ST_GNUPLOT_RGB_FUNCTION_NAME, ST_GNUPLOT_COMMENT_PREFIX);
	}
	
	//plot border
	pipe.write("set border\n");
	
	//data separator
	pipe.write(gnuplotDataSeparatorExpression()+"\n");
	
	//PLOT-TYPE SPECIFIC HEADER
	switch(sources[startSource].plotType){
	case PlotTypeNone:
		rotationInducedScaleX = rotationInducedScaleY = 1.0;
		break;
	case PlotType1DPoints:
	case PlotType1DColoredPoints:
	case PlotType1DVectors:
	case PlotType1DColoredVectors:
		rotationInducedScaleX = rotationInducedScaleY = 1.0;
		//adjust labels, tics & title offsets
		pipe.write("set title offset 0,0\nset xlabel offset 0, 0\nset ylabel offset 0, 0\nset xtic offset 0, 0\nset ytic offset 0, 0\n");
		pipe.write("unset pm3d;\n");
		break;
		
	case PlotType2DPoints:
	case PlotType2DColoredPoints:
	case PlotType2DCurve:
	case PlotType2DCurvePoints:
	case PlotType2DGraphLinear:
	case PlotType2DGraphLinearPoints:
	case PlotType2DGraphSplines:
	case PlotType2DGraphSplinesPoints:
	case PlotType2DVectors:
	case PlotType2DColoredVectors:
	case PlotType2DVariableSizedSpheres:
	case PlotType2DColoredVariableSizedSpheres:
	case PlotType2DFilledCurve:
	case PlotType2DFilledCurves:
	case PlotType2DAnalytic:
		rotationInducedScaleX = rotationInducedScaleY = 1.0;
		//adjust labels, tics & title offsets
		pipe.write("set title offset 0,0\nset xlabel offset 0, 0\nset ylabel offset 0, 0\nset xtic offset 0, 0\nset ytic offset 0, 0\n");
		pipe.write("unset pm3d;\nset autoscale;\n");
		break;
		
	case PlotType3DPoints:
	case PlotType3DColoredPoints:
	case PlotType3DCurve:
	case PlotType3DCurvePoints:
	case PlotType3DVectors:
	case PlotType3DColoredVectors:
		rotationInducedScaleX = 0.73;
		rotationInducedScaleY = 0.7;
		//adjust labels, tics & title offsets
		pipe.write("set title offset 0, 2.6\nset xlabel offset 0, 0\nset ylabel offset 0, 0\nset xtic offset 0, 0\nset ytic offset 0, 0\n");
		pipe.write("unset pm3d;\nset surface;\nunset dgrid3d;\nunset contour;\n");
		pipe.write("set view 60, 30, 1, 1;\nset autoscale;\n");
		break;

	case PlotType1DHistogram:
		rotationInducedScaleX = rotationInducedScaleY = 1.0;
		//adjust labels, tics & title offsets
		pipe.write("set title offset 0,0\nset xlabel offset 0, 0\nset ylabel offset 0, 0\nset xtic offset 0, 0\nset ytic offset 0, 0\n");
		pipe.write("unset pm3d;\nset style data histogram;\nset style histogram clustered gap %e;\nset style fill solid border -1;\n", ST_DEFAULT_ST_GNUPLOT_HISTOGRAM_GAP);
		break;
	
	case PlotType1DEmpiricalDensity:
		rotationInducedScaleX = rotationInducedScaleY = 1.0;
		//adjust labels, tics & title offsets
		pipe.write("set title offset 0,0\nset xlabel offset 0, 0\nset ylabel offset 0, 0\nset xtic offset 0, 0\n");
		pipe.write("unset pm3d;\n");
		break;
		
	case PlotType2DScatteredValueMap:
		rotationInducedScaleX = rotationInducedScaleY = 1.0;
		//adjust labels, tics & title offsets
		pipe.write("set title offset 0,0\nset xlabel offset 0, 0.3\nset ylabel offset -1.5, 0\nset xtic offset 0, 0.3\nset ytic offset 0, 0\n");
		pipe.write("%s\n",contours.gnuplotCommand().c_str());
		pipe.write((showColorMap ? "set surface;\n" : "unset surface;\n"));
		switch(gridColoringType){
		case GridColoringTypeMean:
			pipe.write("set pm3d map corners2color mean;\n");
			break;
		case GridColoringTypeMax:
			pipe.write("set pm3d map corners2color max;\n");
			break;
		case GridColoringTypeMin:
			pipe.write("set pm3d map corners2color min;\n");
			break;
		case GridColoringTypeLL:
			pipe.write("set pm3d map corners2color c1;\n");
			break;
		case GridColoringTypeUR:
			pipe.write("set pm3d map corners2color c4;\n");
			break;
		}	
		//grid interpolating (smoothing)
		switch(gridSmoothingType){
		case GridSmoothingTypeEuclidean:
			pipe.write(	"set dgrid3d %i,%i qnorm %e\n", 
						xGridSize, 
						yGridSize, 
						ST_DEFAULT_ST_GNUPLOT_GRIDNORM);
			break;
		case GridSmoothingTypeGauss:
			pipe.write(	"set dgrid3d %i,%i gauss %e, %e\n", 
						xGridSize, 
						yGridSize, 
						xDistanceScale, 
						yDistanceScale); 
			break;
		case GridSmoothingTypeBox:
			pipe.write(	"set dgrid3d %i,%i box %e, %e\n", 
						xGridSize, 
						yGridSize, 
						xDistanceScale, 
						yDistanceScale); 
			break;
		}
		break;
		
	case PlotType2DGridValueMap:
		rotationInducedScaleX = rotationInducedScaleY = 1.0;
		//pipe.write("set size %e, %e\n", availableWidth * (showColorBox ? 1.1 : 1.23), availableHeight * 1.16);
		//pipe.write("set origin %e, %e\n", originX - 0.11 * availableWidth, originY - 0.073 * availableHeight);
		pipe.write("set title offset 0,0\nset xlabel offset 0, 0.3\nset ylabel offset -1.5, 0\nset xtic offset 0, 0.3\nset ytic offset 0, 0\n");
		pipe.write("unset dgrid3d;\n");
		eRange.scaleRangeByFactor(1.0 + ST_DEFAULT_PLOT_RANGE_TOLLERANCE);
		pipe.write("%s\n", contours.gnuplotCommand().c_str());
		pipe.write((showColorMap ? "set surface;\n" : "unset surface;\n"));
		switch(gridColoringType){
		case GridColoringTypeMean:
			pipe.write("set pm3d map corners2color mean;\n");
			break;
		case GridColoringTypeMax:
			pipe.write("set pm3d map corners2color max;\n");
			break;
		case GridColoringTypeMin:
			pipe.write("set pm3d map corners2color min;\n");
			break;
		case GridColoringTypeLL:
			pipe.write("set pm3d map corners2color c1;\n");
			break;
		case GridColoringTypeUR:
			pipe.write("set pm3d map corners2color c4;\n");
			break;
		}
		break;
		
	case PlotType2DScannedValueMap:
		rotationInducedScaleX = rotationInducedScaleY = 1.0;
		//pipe.write("set size %e, %e\n", availableWidth * (showColorBox ? 1.1 : 1.23), availableHeight * 1.16);
		//pipe.write("set origin %e, %e\n", originX - 0.11 * availableWidth, originY - 0.073 * availableHeight);
		pipe.write("set title offset 0,0\nset xlabel offset 0, 0.3\nset ylabel offset -1.5, 0\nset xtic offset 0, 0.3\nset ytic offset 0, 0\n");
		pipe.write("unset dgrid3d;\n");
		pipe.write("%s\n",contours.gnuplotCommand().c_str());
		pipe.write((showColorMap ? "set surface;\n" : "unset surface;\n"));
		switch(gridColoringType){
		case GridColoringTypeMean:
			pipe.write("set pm3d map corners2color mean;\n");
			break;
		case GridColoringTypeMax:
			pipe.write("set pm3d map corners2color max;\n");
			break;
		case GridColoringTypeMin:
			pipe.write("set pm3d map corners2color min;\n");
			break;
		case GridColoringTypeLL:
			pipe.write("set pm3d map corners2color c1;\n");
			break;
		case GridColoringTypeUR:
			pipe.write("set pm3d map corners2color c4;\n");
			break;
		}
		break;
		
	case PlotType2DColorImage:
	case PlotType2DValueImage:
		rotationInducedScaleX = rotationInducedScaleY = 1.0;
		//color conversion formulas
		pipe.write("div(x,y) = floor(x/y); mod(x,y) = x - floor(x/y) * y;\n");
		//int to (r,g,b) \in [0,1]^3
		pipe.write("%s(v) = div(v, 65536)/255.0; %sred component of inverse of %s()\n", ST_GNUPLOT_INVERSE_RGB_FUNCTION_NAME_RED,ST_GNUPLOT_COMMENT_PREFIX,ST_GNUPLOT_RGB_FUNCTION_NAME);
		pipe.write("%s(v) = div(mod(v, 65536), 256)/255.0; %sgreen component of inverse of %s()\n", ST_GNUPLOT_INVERSE_RGB_FUNCTION_NAME_GREEN,ST_GNUPLOT_COMMENT_PREFIX,ST_GNUPLOT_RGB_FUNCTION_NAME);
		pipe.write("%s(v) = mod(v, 256)/255.0; %sblue component of inverse of %s()\n", ST_GNUPLOT_INVERSE_RGB_FUNCTION_NAME_BLUE,ST_GNUPLOT_COMMENT_PREFIX,ST_GNUPLOT_RGB_FUNCTION_NAME);
		//adjust labels, tics & title offsets
		pipe.write("set title offset 0,0\nset xlabel offset 0, 0\nset ylabel offset 0, 0\nset xtic offset 0, 0\nset ytic offset 0, 0\n");
		//steps and ranges
		xStep = (sources[startSource].xCount<=1 ? 0 : (sources[startSource].maxX - sources[startSource].minX)/(sources[startSource].xCount - 1.0));
		yStep = (sources[startSource].yCount<=1 ? 0 : (sources[startSource].maxY - sources[startSource].minY)/(sources[startSource].yCount - 1.0));
		eRange.mergeWithRange(MultiRange(	sources[startSource].minX - 0.5*xStep, 
											sources[startSource].maxX + 0.5*xStep, 
											sources[startSource].minY - 0.5*yStep, 
											sources[startSource].maxY + 0.5*yStep)); 
		break;
		
	case PlotType3DScatteredSurface:
		rotationInducedScaleX = 0.73;
		rotationInducedScaleY = 0.7;
		//adjust labels, tics & title offsets
		pipe.write("set title offset 0,2.6\nset xlabel offset 0, 0\nset ylabel offset 0, 0\nset xtic offset 0, 0\nset ytic offset 0, 0\n");
		pipe.write("unset hidden3d;\nset surface;\nset view 60, 30, 1, 1;\n");
		pipe.write("%s\n",contours.gnuplotCommand().c_str());
		switch(gridColoringType){
		case GridColoringTypeMean:
			pipe.write("set pm3d at s corners2color mean;\n");
			break;
		case GridColoringTypeMax:
			pipe.write("set pm3d at s corners2color max;\n");
			break;
		case GridColoringTypeMin:
			pipe.write("set pm3d at s corners2color min;\n");
			break;
		case GridColoringTypeLL:
			pipe.write("set pm3d at s corners2color c1;\n");
			break;
		case GridColoringTypeUR:
			pipe.write("set pm3d at s corners2color c4;\n");
			break;
		}	
		//grid interpolating (smoothing)
		switch(gridSmoothingType){
		case GridSmoothingTypeEuclidean:
			pipe.write(	"set dgrid3d %i,%i qnorm %e\n", 
						xGridSize, 
						yGridSize, 
						ST_DEFAULT_ST_GNUPLOT_GRIDNORM);
			break;
		case GridSmoothingTypeGauss:
			pipe.write(	"set dgrid3d %i,%i gauss %e, %e\n", 
						xGridSize, 
						yGridSize, 
						xDistanceScale, 
						yDistanceScale); 
			break;
		case GridSmoothingTypeBox:
			pipe.write(	"set dgrid3d %i,%i box %e, %e\n", 
						xGridSize, 
						yGridSize, 
						xDistanceScale, 
						yDistanceScale); 
			break;
		}
		break;
		
	case PlotType3DGridSurface:	
	case PlotType3DScannedSurface:
	case PlotType3DValuedScannedSurface:
		rotationInducedScaleX = 0.73;
		rotationInducedScaleY = 0.7;
		//adjust labels, tics & title offsets
		pipe.write("set title offset 0,2.6\nset xlabel offset 0, 0\nset ylabel offset 0, 0\nset xtic offset 0, 0\nset ytic offset 0, 0\n");
		pipe.write("unset dgrid3d; unset hidden3d;\nset surface;\nset view 60, 30, 1, 1;\n");
		pipe.write("%s\n",contours.gnuplotCommand().c_str());
		switch(gridColoringType){
		case GridColoringTypeMean:
			pipe.write("set pm3d at s corners2color mean;\n");
			break;
		case GridColoringTypeMax:
			pipe.write("set pm3d at s corners2color max;\n");
			break;
		case GridColoringTypeMin:
			pipe.write("set pm3d at s corners2color min;\n");
			break;
		case GridColoringTypeLL:
			pipe.write("set pm3d at s corners2color c1;\n");
			break;
		case GridColoringTypeUR:
			pipe.write("set pm3d at s corners2color c4;\n");
			break;
		}
		break;
	default:
		rotationInducedScaleX = rotationInducedScaleY = 1.0;
	}
	
	
	//set origin and size of plot by means of margins
	//only write this for the first source
	if(startSource == 0){		
		double colorBoxInducedMarginX = (showColorBox ? ST_DEFAULT_RELATIVE_COLORBOX_WIDTH + ST_DEFAULT_RELATIVE_COLORBOX_DIST + ST_DEFAULT_RELATIVE_COLORBOX_RMARGIN : 0.0);
		double rotationInducedMarginX = (1.0 - ST_DEFAULT_RELATIVE_LMARGIN - ST_DEFAULT_RELATIVE_RMARGIN) * (1.0 - rotationInducedScaleX);
		double rotationInducedMarginY = (1.0 - ST_DEFAULT_RELATIVE_TMARGIN - ST_DEFAULT_RELATIVE_BMARGIN) * (1.0 - rotationInducedScaleY);
		lmargin	= (originX + resizeNonPlotAreaWidthByFactor*availableWidth * (ST_DEFAULT_RELATIVE_LMARGIN  + rotationInducedMarginX/2.0));
		rmargin = (originX + availableWidth * (1.0 - resizeNonPlotAreaWidthByFactor*(ST_DEFAULT_RELATIVE_RMARGIN + rotationInducedMarginX/2.0 + colorBoxInducedMarginX)));
		bmargin = (originY + resizeNonPlotAreaHeightByFactor*availableHeight * (ST_DEFAULT_RELATIVE_BMARGIN + rotationInducedMarginY/2.0));
		tmargin	= (originY + availableHeight * (1.0 - resizeNonPlotAreaHeightByFactor*(ST_DEFAULT_RELATIVE_TMARGIN + rotationInducedMarginY/2.0)));
		pipe.write("set lmargin at screen %e\nset rmargin at screen %e\n", lmargin, rmargin);
		pipe.write("set bmargin at screen %e\nset tmargin at screen %e\n", bmargin, tmargin);
				
		//color box
		if(showColorBox){
			pipe.write("set colorbox vertical user origin %e, %e size %e, %e;\n", 
						originX + availableWidth * (1.0 - resizeNonPlotAreaWidthByFactor*(ST_DEFAULT_RELATIVE_RMARGIN + ST_DEFAULT_RELATIVE_COLORBOX_WIDTH + ST_DEFAULT_RELATIVE_COLORBOX_RMARGIN)), 
						bmargin, 
						resizeNonPlotAreaWidthByFactor * availableWidth * ST_DEFAULT_RELATIVE_COLORBOX_WIDTH, 
						tmargin - bmargin);
		}else{
			pipe.write("unset colorbox;\n");
		}
	}
	
	//Only to be set for the first plot round to avoid duplicate (and thus fuzzy) printings
	if(startSource == 0){
		pipe.write("set title \"%s%s\";\n", (tag=="" ? "" : ("("+tag+") ").c_str()), formattedTitle.c_str());
		pipe.write("%s\n", axisLabels.gnuplotCommand(interpretSymbols, escapeUnderscoresInTicsAndLabels).c_str());	//axis labels
		string formattedCBLabel = STPlot_newLineEscapedVersion((interpretSymbols ? replaceSymbols(CBLabel) : CBLabel));
		if(escapeUnderscoresInTicsAndLabels) STPlot_aux_replaceAll(formattedCBLabel, "_", "\\\\_");
		pipe.write("set cblabel \"%s\";\n", formattedCBLabel.c_str()); //ColorBox label
		pipe.write((showKeys ? 		"set key;\n" 		: "unset key;\n"));		//keys	
		pipe.write("%s\n", CBTics.gnuplotCommand("cbtics", interpretSymbols, escapeUnderscoresInTicsAndLabels, resizeTicsByFactor).c_str());
		
		//custom labels
		for(long l=0; l<labels.size(); ++l){
			commandJustAfterFirstPlot += "; unset label " + STPlot_makeString(nextAvailLabelTag);
			pipe.write(labels[l].gnuplotCommand(nextAvailLabelTag++, interpretSymbols, dim, windowFont, escapeUnderscoresInTicsAndLabels)+"\n");
		}
		
		//background color
		if(!backgroundColor.isWhite()){
			pipe.write("set object %i rectangle from screen %e,%e to screen %e,%e fillcolor rgb '#%s' behind;\n",ST_DEFAULT_BACKGROUND_TAG, lmargin,bmargin,rmargin,tmargin,backgroundColor.getHex().c_str());
		}
			
		//tics. Should come after the offsets, since the latter enables tics.
		pipe.write("%s\n%s\n%s\n", xTics.gnuplotCommand("xtics", interpretSymbols, escapeUnderscoresInTicsAndLabels, resizeTicsByFactor).c_str(), yTics.gnuplotCommand("ytics", interpretSymbols, escapeUnderscoresInTicsAndLabels, resizeTicsByFactor).c_str(), zTics.gnuplotCommand("ztics", interpretSymbols, escapeUnderscoresInTicsAndLabels, resizeTicsByFactor).c_str());	
		
		//better tics format than default
		pipe.write("set format \"%%g\"\n");
	
	}else{
		pipe.write("unset title;\n");
		pipe.write("unset xtics; unset ytics; unset ztics; unset cbtics\n");
		pipe.write("unset xlabel; unset ylabel; unset zlabel; unset cblabel;\n");
		pipe.write("unset key\n"); //only the first compatibility group of sources can display keys, everything additional will create a pile of mess!
	}
	
	//CB Range
	pipe.write("set cbrange %s;\n", CBRange.gnuplotExpression().c_str());
	//Ranges. Should come after the offsets and styles.
	pipe.write("%s\n",eRange.gnuplotCommand().c_str());
	//logarithmic axis scales
	pipe.write("%s log x; %s log y; %s log z;\n", (xLogarithmic ? "set" : "unset"), (yLogarithmic ? "set" : "unset"), (zLogarithmic ? "set" : "unset"));

	//PlotType-specific plot command
	switch(sources[startSource].plotType){
	case PlotTypeNone:
		//do nothing
		break;
	case PlotType1DPoints:
	case PlotType1DColoredPoints:
	case PlotType1DVectors:
	case PlotType1DColoredVectors:
	case PlotType2DPoints:
	case PlotType2DColoredPoints:
	case PlotType2DCurve:
	case PlotType2DCurvePoints:
	case PlotType2DGraphLinear:
	case PlotType2DGraphLinearPoints:
	case PlotType2DGraphSplines:
	case PlotType2DGraphSplinesPoints:
	case PlotType2DVectors:
	case PlotType2DColoredVectors:
	case PlotType2DVariableSizedSpheres:
	case PlotType2DColoredVariableSizedSpheres:
	case PlotType2DFilledCurve:
	case PlotType2DFilledCurves:
	case PlotType2DAnalytic:
		pipe.write("plot\t");
		break;
		
	case PlotType3DPoints:
	case PlotType3DColoredPoints:
	case PlotType3DVectors:
	case PlotType3DColoredVectors:
	case PlotType3DCurve:
	case PlotType3DCurvePoints:
		pipe.write("splot\t");
		break;

	case PlotType1DHistogram:
	case PlotType1DEmpiricalDensity:
		pipe.write("plot\t");
		break;
		
	case PlotType2DScatteredValueMap:
	case PlotType2DGridValueMap:
	case PlotType2DScannedValueMap:
		pipe.write(	"splot\t");
		break;
		
	case PlotType2DColorImage:
	case PlotType2DValueImage:
		pipe.write(	"plot\t");
		break;
		
	case PlotType3DScatteredSurface:
	case PlotType3DGridSurface:
	case PlotType3DScannedSurface:
	case PlotType3DValuedScannedSurface:
		pipe.write("splot\t");
		break;
	}
	
	
	//List plot-sources
	error = PlotErrorNone;
	for(unsigned int s=startSource; s<sources.size(); ++s){
		if(s >= startSource + 1){
			if(!plotTypesCompatible(sources[s].plotType, sources[s-1].plotType)){
				pipe.write("%sChanging plot style within same plot\n", ST_GNUPLOT_COMMENT_PREFIX);
				if(commandJustAfterFirstPlot != ""){
					pipe.write("%s%s\n%s\n", ST_GNUPLOT_COMMENT_PREFIX, "Calling some closing commands after first plot", commandJustAfterFirstPlot.c_str());
				}
				writeCommandToPipe(	pipe, 
									originX, 
									originY, 
									availableWidth, 
									availableHeight, 
									resizeNonPlotAreaWidthByFactor,
									resizeNonPlotAreaHeightByFactor, 
									resizeTicsByFactor,
									windowFont,
									interpretSymbols,
									s, 
									tag,
									escapeUnderscoresInTitle,
									escapeUnderscoresInTicsAndLabels,
									"",
									nextAvailLabelTag);
				break;
			}else{
				pipe.write(ST_DEFAULT_WRITE_COMPOSITE_PLOTS_MULTILINE ? ",\\\n\t\t" : ", ");
			}
		}
		
		if(sources[s].title == ""){
			sourceTitle = "notitle";
		}else{
			sourceTitle = STPlot_newLineEscapedVersion(interpretSymbols ? replaceSymbols(sources[s].title) : sources[s].title);
			if(escapeUnderscoresInTicsAndLabels){ STPlot_aux_replaceAll(sourceTitle,"_","\\\\_"); }
			sourceTitle = "title \"" + sourceTitle + "\"";
		}
		switch(sources[s].plotType){
		case PlotTypeNone:
			//nothing to plot here
			break;
		case PlotType1DPoints:
			pipe.write(	"'%s' index %i using (%s):(0) with points %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotPointTypeExpression().c_str(),
						sources[s].gnuplotPointSizeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType1DColoredPoints:
			pipe.write(	"'%s' index %i using (%s):(0):(%s($%i, $%i, $%i)) with points %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						ST_GNUPLOT_RGB_FUNCTION_NAME, 
						sources[s].columns[1], 
						sources[s].columns[2], 
						sources[s].columns[3],
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotPointTypeExpression().c_str(),
						sources[s].gnuplotPointSizeExpression().c_str(),
						sourceTitle.c_str());
			break;
		
		case PlotType1DVectors:
			pipe.write(	"'%s' index %i using (%s):(0):%i:(0) with vectors %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						sources[s].columns[1], 
						vectorHeadStyle2string(sources[s].vectorHeadStyle, sources[s].vectorHeadSize).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType1DColoredVectors:
			pipe.write(	"'%s' index %i using (%s):(0):%i:(0):(%s($%i, $%i, $%i)) with vectors %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						sources[s].columns[1], 
						ST_GNUPLOT_RGB_FUNCTION_NAME, 
						sources[s].columns[2], 
						sources[s].columns[3], 
						sources[s].columns[4],
						vectorHeadStyle2string(sources[s].vectorHeadStyle, sources[s].vectorHeadSize).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DPoints:
			pipe.write(	"'%s' index %i using (%s):(%s) with points %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotPointTypeExpression().c_str(),
						sources[s].gnuplotPointSizeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DColoredPoints:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s($%i, $%i, $%i)) with points %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						ST_GNUPLOT_RGB_FUNCTION_NAME, 
						sources[s].columns[2], 
						sources[s].columns[3], 
						sources[s].columns[4],
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotPointTypeExpression().c_str(), 
						sources[s].gnuplotPointSizeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DVectors:
			pipe.write(	"'%s' index %i using (%s):(%s):%i:%i with vectors %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].columns[2], 
						sources[s].columns[3],
						vectorHeadStyle2string(sources[s].vectorHeadStyle, sources[s].vectorHeadSize).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DColoredVectors:
			pipe.write(	"'%s' index %i using (%s):(%s):%i:%i:(%s($%i, $%i, $%i)) with vectors %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].columns[2], 
						sources[s].columns[3],
						ST_GNUPLOT_RGB_FUNCTION_NAME, 
						sources[s].columns[4], 
						sources[s].columns[5], 
						sources[s].columns[6], 
						vectorHeadStyle2string(sources[s].vectorHeadStyle, sources[s].vectorHeadSize).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DVariableSizedSpheres:
			pipe.write(	"'%s' index %i using (%s):(%s):%i with circles fillstyle solid %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].columns[2], 
						sources[s].gnuplotColorExpression().c_str(),
						sourceTitle.c_str()); 
			break;
	
		case PlotType2DColoredVariableSizedSpheres:
			pipe.write(	"'%s' index %i using (%s):(%s):%i:(rgb($%i,$%i,$%i)) with circles fillstyle solid %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].columns[2], 
						sources[s].columns[3], 
						sources[s].columns[4], 
						sources[s].columns[5], 
						sources[s].gnuplotColorExpression().c_str(),
						sourceTitle.c_str()); 
			break;
			
		case PlotType2DFilledCurve:
			pipe.write(	"'%s' index %i using (%s):(%s):(%g) with filledcurves %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].fillReference, 
						sources[s].gnuplotFilledCurveStyleExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sourceTitle.c_str()); 
			break;
			
		case PlotType2DFilledCurves:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s) with filledcurves %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						sources[s].gnuplotFilledCurveStyleExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sourceTitle.c_str()); 
			break;
			
		case PlotType2DGraphLinear:
			pipe.write(	"'%s' index %i using (%s):(%s) with lines %s %s %s %s",
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DGraphLinearPoints:
			pipe.write(	"'%s' index %i using (%s):(%s) with linespoints %s %s %s %s %s %s",
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(),
						sources[s].gnuplotPointTypeExpression().c_str(),
						sources[s].gnuplotPointSizeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DCurve:
			pipe.write(	"'%s' index %i using (%s):(%s) with lines %s %s %s %s",
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DCurvePoints:
			pipe.write(	"'%s' index %i using (%s):(%s) with linespoints %s %s %s %s %s %s",
						sources[s].filePath.c_str(), 
						sources[s].index,
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].gnuplotColorExpression().c_str(), 
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(),
						sources[s].gnuplotPointTypeExpression().c_str(),
						sources[s].gnuplotPointSizeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DGraphSplines:
			pipe.write(	"'%s' index %i using (%s):(%s) smooth csplines %s %s %s %s",
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(), 
						sourceTitle.c_str());
			break;
			
		case PlotType2DGraphSplinesPoints:
			pipe.write(	"'%s' index %i using (%s):(%s) smooth csplines %s %s %s notitle, '' index %i using (%s):(%s) with points %s %s %s %s",
						sources[s].filePath.c_str(), 
						sources[s].index,
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotPointTypeExpression().c_str(),
						sources[s].gnuplotPointSizeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DAnalytic:
			pipe.write(	"%s %s %s %s %s",
						sources[s].filePath.c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType3DPoints:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s) with points %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotPointTypeExpression().c_str(),
						sources[s].gnuplotPointSizeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType3DColoredPoints:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s):(%s($%i, $%i, $%i)) with points %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						ST_GNUPLOT_RGB_FUNCTION_NAME, 
						sources[s].columns[3], 
						sources[s].columns[4], 
						sources[s].columns[5], 
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotPointTypeExpression().c_str(),
						sources[s].gnuplotPointSizeExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType3DVectors:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s):%i:i:i with vectors %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						sources[s].columns[3], 
						sources[s].columns[4], 
						sources[s].columns[5], 
						vectorHeadStyle2string(sources[s].vectorHeadStyle, sources[s].vectorHeadSize).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(), 
						sources[s].gnuplotLineTypeExpression().c_str(), 
						sourceTitle.c_str());
			break;
		
		case PlotType3DColoredVectors:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s):%i:%i:%i:(%s($%i, $%i, $%i)) with vectors %s %s %s %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						sources[s].columns[3], 
						sources[s].columns[4], 
						sources[s].columns[5], 
						ST_GNUPLOT_RGB_FUNCTION_NAME, 
						sources[s].columns[6], 
						sources[s].columns[7], 
						sources[s].columns[8], 
						vectorHeadStyle2string(sources[s].vectorHeadStyle, sources[s].vectorHeadSize).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(), 
						sourceTitle.c_str());
			break;
			
		case PlotType3DCurve:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s) with lines %s %s %s %s",
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(),   
						sourceTitle.c_str());
			break;
			
		case PlotType3DCurvePoints:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s) with linespoints %s %s %s %s %s %s",
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						sources[s].gnuplotLineWidthExpression().c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sources[s].gnuplotLineTypeExpression().c_str(), 
						sources[s].gnuplotPointTypeExpression().c_str(),
						sources[s].gnuplotPointSizeExpression().c_str(),
						sourceTitle.c_str());
			break;
	
		case PlotType1DHistogram:
			pipe.write(	"'%s' index %i using (%s):xticlabels(%i) %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(), sources[s].columns[0], 
						sources[s].gnuplotColorExpression().c_str(),
						sourceTitle.c_str());
			break;
		
		case PlotType1DEmpiricalDensity:
			pipe.write(	"'%s' index %i using (%s):(%s) smooth unique with boxes fillstyle solid %s %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						sources[s].gnuplotColorExpression().c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DScatteredValueMap:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s) %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DGridValueMap:
			switch(gridColoringType){
			case GridColoringTypeMean:
				pipe.write(	"'%s' index %i using ($1*%e+%e):($2*%e+%e):(%s) matrix with pm3d %s", 
							sources[s].filePath.c_str(), 
							sources[s].index,
							sources[s].xStep, 
							sources[s].xStart,
							sources[s].yStep, 
							sources[s].yStart, 
							evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
							sourceTitle.c_str());
				break;
				
			case GridColoringTypeMax:
				pipe.write(	"'%s' index %i using ($1*%e+%e):($2*%e+%e):(%s) matrix with pm3d %s", 
							sources[s].filePath.c_str(), 
							sources[s].index,
							sources[s].xStep, 
							sources[s].xStart, 
							sources[s].yStep, 
							sources[s].yStart, 
							evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
							sourceTitle.c_str());
				break;
				
			case GridColoringTypeMin:
				pipe.write(	"'%s' index %i using ($1*%e+%e):($2*%e+%e):(%s) matrix with pm3d %s", 
							sources[s].filePath.c_str(), 
							sources[s].index, 
							sources[s].xStep, 
							sources[s].xStart,
							sources[s].yStep, 
							sources[s].yStart, 
							evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
							sourceTitle.c_str());
				break;
			
			case GridColoringTypeLL:
				pipe.write(	"'%s' index %i using (($1-0.5)*%e+%e):(($2-0.5)*%e+%e):(%s) matrix with pm3d %s", 
							sources[s].filePath.c_str(), 
							sources[s].index,
							sources[s].xStep, 
							sources[s].xStart,
							sources[s].yStep, 
							sources[s].yStart, 
							evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
							sourceTitle.c_str());
				break;
			
			case GridColoringTypeUR:
				pipe.write(	"'%s' index %i using (($1+0.5)*%e+%e):(($2+0.5)*%e+%e):(%s) matrix with pm3d %s", 
							sources[s].filePath.c_str(), 
							sources[s].index,
							sources[s].xStep, 
							sources[s].xStart,
							sources[s].yStep, 
							sources[s].yStart, 
							evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
							sourceTitle.c_str());
				break;
			}
			break;
			
		case PlotType2DScannedValueMap:	
			pipe.write(	"'%s' index %i using (%s):(%s):(%s) %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType2DColorImage:
			if((sources[s].xCount <= 1) || (sources[s].yCount <= 1)){
				error = PlotErrorNothingToBeDone;
				break;
			}
			xStep = (sources[s].maxX - sources[s].minX)/(sources[s].xCount - 1.0);
			yStep = (sources[s].maxY - sources[s].minY)/(sources[s].yCount - 1.0);
			pipe.write(	"'%s' index %i using (%e+%e*$1):(%e+%e*$2):(255 * %s($3)):(255 * %s($3)):(255 * %s($3)) matrix with rgbimage %s\n", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						sources[s].minX, 
						xStep,
						sources[s].minY, 
						yStep,
						ST_GNUPLOT_INVERSE_RGB_FUNCTION_NAME_RED,
						ST_GNUPLOT_INVERSE_RGB_FUNCTION_NAME_GREEN,
						ST_GNUPLOT_INVERSE_RGB_FUNCTION_NAME_BLUE,
						sourceTitle.c_str());
			break;
			
		case PlotType2DValueImage:
			if((sources[s].xCount <= 1) || (sources[s].yCount <= 1)){
				error = PlotErrorNothingToBeDone;
				break;
			}
			xStep = (sources[s].maxX - sources[s].minX)/(sources[s].xCount - 1.0);
			yStep = (sources[s].maxY - sources[s].minY)/(sources[s].yCount - 1.0);
			pipe.write(	"'%s' index %i using (%e+%e*$1):(%e+%e*$2):(%s) matrix with image %s\n", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						sources[s].minX, 
						xStep,
						sources[s].minY, 
						yStep,
						evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType3DScatteredSurface:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s) with pm3d %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						sources[s].columns[0], 
						sources[s].columns[1], 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						sourceTitle.c_str());
			break;
			
		case PlotType3DGridSurface:
			switch(gridColoringType){
			case GridColoringTypeMean:
				pipe.write(	"'%s' index %i using ($1*%e+%e):($2*%e+%e):(%s) matrix with pm3d %s", 
							sources[s].filePath.c_str(), 
							sources[s].index,
							sources[s].xStep, 
							sources[s].xStart,
							sources[s].yStep, 
							sources[s].yStart, 
							evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
							sourceTitle.c_str());
				break;
				
			case GridColoringTypeMax:
				pipe.write(	"'%s' index %i using ($1*%e+%e):($2*%e+%e):(%s) matrix with pm3d %s", 
							sources[s].filePath.c_str(), 
							sources[s].index,
							sources[s].xStep, 
							sources[s].xStart,
							sources[s].yStep, 
							sources[s].yStart, 
							evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
							sourceTitle.c_str());
				break;
				
			case GridColoringTypeMin:
				pipe.write(	"'%s' index %i using ($1*%e+%e):($2*%e+%e):(%s) matrix with pm3d %s", 
							sources[s].filePath.c_str(), 
							sources[s].index, 
							sources[s].xStep, 
							sources[s].xStart, 
							sources[s].yStep, 
							sources[s].yStart, 
							evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
							sourceTitle.c_str());
				break;
			
			case GridColoringTypeLL:
				pipe.write(	"'%s' index %i using (($1-0.5)*%e+%e):(($2-0.5)*%e+%e):(%s) matrix with pm3d %s", 
							sources[s].filePath.c_str(), 
							sources[s].index, 
							sources[s].xStep, 
							sources[s].xStart, 
							sources[s].yStep, 
							sources[s].yStart, 
							evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
							sourceTitle.c_str());
				break;
			
			case GridColoringTypeUR:
				pipe.write(	"'%s' index %i using (($1+0.5)*%e+%e):(($2+0.5)*%e+%e):(%s) matrix with pm3d %s", 
							sources[s].filePath.c_str(), 
							sources[s].index, 
							sources[s].xStep, 
							sources[s].xStart, 
							sources[s].yStep, 
							sources[s].yStart, 
							evaluateFunctionForColumn(sources[s].zTransformationFunction, 3).c_str(),
							sourceTitle.c_str());
				break;
			}
			break;
			
		case PlotType3DScannedSurface:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s) with pm3d %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						sourceTitle.c_str());
			break;
		
		case PlotType3DValuedScannedSurface:
			pipe.write(	"'%s' index %i using (%s):(%s):(%s):(%s) with pm3d %s", 
						sources[s].filePath.c_str(), 
						sources[s].index, 
						evaluateFunctionForColumn(sources[s].xTransformationFunction, sources[s].columns[0]).c_str(),
						evaluateFunctionForColumn(sources[s].yTransformationFunction, sources[s].columns[1]).c_str(),
						evaluateFunctionForColumn(sources[s].zTransformationFunction, sources[s].columns[2]).c_str(),
						evaluateFunctionForColumn(sources[s].wTransformationFunction, sources[s].columns[3]).c_str(),
						sourceTitle.c_str());
			break;
		}
	}
	
	// force appearance of color box
	if(showColorBox && (!hasPlotWhichShowsColorbox())){
		pipe.write(ST_DEFAULT_WRITE_COMPOSITE_PLOTS_MULTILINE ? ",\\\n\t\t" : ", ");
		pipe.write("1/0 lc palette notitle %s Needed to make color box appear", ST_GNUPLOT_COMMENT_PREFIX);
	}
	
	CLOSING:
	pipe.write("\n");
	if(commandJustAfterFirstPlot != ""){
		pipe.write("\n%s%s\n%s\n\n", ST_GNUPLOT_COMMENT_PREFIX, "Calling some closing commands after first plot", commandJustAfterFirstPlot.c_str());
	}
	if(!backgroundColor.isWhite()){
		//If background object is not removed, it will be drawn multiple times with any subsequent plots
		pipe.write("unset object %i;\n",ST_DEFAULT_BACKGROUND_TAG);
	}
	if(error == PlotErrorIncompatiblePlots){
		pipe.write("%sWarning: Incompatible plots requested. Plot incomplete.\n", ST_GNUPLOT_COMMENT_PREFIX);
	}
	pipe.write("\n");
	return error;
}





#pragma mark -
#pragma mark MultiPlot
#pragma mark -


//constructors
MultiPlot::MultiPlot(){
	setDefaults();
	setOutputToTerminal(ST_DEFAULT_TERM,ST_DEFAULT_TERMINAL_WIDTH,ST_DEFAULT_TERMINAL_HEIGHT,"");
	pipe.open();
}



MultiPlot::MultiPlot(const string &gnuplotPath){
	setDefaults();
	setOutputToTerminal(ST_DEFAULT_TERM,ST_DEFAULT_TERMINAL_WIDTH,ST_DEFAULT_TERMINAL_HEIGHT,"");
	if(gnuplotPath==""){ pipe.open(); }
	else{ pipe.open(gnuplotPath); }
}


MultiPlot::MultiPlot(const Plot &_plot, const string &_title){
	setDefaults();
	setOutputToTerminal(ST_DEFAULT_TERM,ST_DEFAULT_TERMINAL_WIDTH,ST_DEFAULT_TERMINAL_HEIGHT, _title);
	plots.push_back(_plot);
	pipe.open();
}


MultiPlot::MultiPlot(const Plot &_plot, const string &_title, const string &gnuplotPath){
	setDefaults();
	setOutputToTerminal(ST_DEFAULT_TERM, ST_DEFAULT_TERMINAL_WIDTH,ST_DEFAULT_TERMINAL_HEIGHT, _title);
	plots.push_back(_plot);
	if(gnuplotPath==""){ pipe.open(); }
	else{ pipe.open(gnuplotPath); }
}


MultiPlot::MultiPlot(const vector<Plot> &_plots, string _title){
	setDefaults();
	setOutputToTerminal(ST_DEFAULT_TERM,ST_DEFAULT_TERMINAL_WIDTH,ST_DEFAULT_TERMINAL_HEIGHT,_title);
	plots = _plots;
	pipe.open();
}
	
	
MultiPlot::MultiPlot(const vector<Plot> &_plots, string _title, string gnuplotPath){
	setDefaults();
	setOutputToTerminal(ST_DEFAULT_TERM,ST_DEFAULT_TERMINAL_WIDTH,ST_DEFAULT_TERMINAL_HEIGHT,_title);
	plots = _plots;
	pipe.open(gnuplotPath);	
}

void MultiPlot::setDefaults(){
	setAspectRatio(ST_DEFAULT_MULTIPLOT_RATIO, ST_DEFAULT_MIN_MULTIPLOT_RATIO, ST_DEFAULT_MAX_MULTIPLOT_RATIO);
	interpretSymbols 	= true;
	tagPlots 			= false;
	onlyResizePlotAreaInResizedPlots = false;
	escapeUnderscoresInFootnote = escapeUnderscoresInTitle = false;
	plotHSpacing 		= ST_DEFAULT_PLOT_H_SPACING;
	plotVSpacing 		= ST_DEFAULT_PLOT_V_SPACING;
}


void MultiPlot::setOutputToTerminal(const string &_terminalName, double width, double height, const string &_title){
	terminalName 		= _terminalName;
	filePath			= "";
	plotWidth			= width;
	plotHeight			= height;
	defaultPlotWidth 	= ST_DEFAULT_TERMINAL_WIDTH;
	defaultPlotHeight 	= ST_DEFAULT_TERMINAL_HEIGHT;
	terminalTitle		= _title;
	outputType			= PlotOutputTypeTerminal;
}

void MultiPlot::setOutputToPDF(const string &_filePath, double width, double height){
	filePath			= _filePath;
	terminalName		= "";
	plotWidth			= width;
	plotHeight			= height;
	defaultPlotWidth 	= ST_DEFAULT_PDF_WIDTH;
	defaultPlotHeight	= ST_DEFAULT_PDF_HEIGHT;
	terminalTitle		= "";
	outputType			= PlotOutputTypeFilePDF;
}


void MultiPlot::setOutputToEPS(const string &_filePath, double width, double height){
	filePath			= _filePath;
	terminalName		= "";
	plotWidth			= width;
	plotHeight			= height;
	defaultPlotWidth 	= ST_DEFAULT_EPS_WIDTH;
	defaultPlotHeight 	= ST_DEFAULT_EPS_HEIGHT;
	terminalTitle		= "";
	outputType			= PlotOutputTypeFileEPS;
}


void MultiPlot::setOutputToSVG(const string &_filePath, double width, double height){
	filePath			= _filePath;
	terminalName		= "";
	plotWidth			= width;
	plotHeight			= height;
	defaultPlotWidth 	= ST_DEFAULT_SVG_WIDTH;
	defaultPlotHeight 	= ST_DEFAULT_SVG_HEIGHT;
	terminalTitle		= "";
	outputType			= PlotOutputTypeFileSVG;
}


void MultiPlot::setPlotSpacing(double Hspacing, double Vspacing){
	plotHSpacing = Hspacing;
	plotVSpacing = Vspacing;
}

bool MultiPlot::program_exists(const string &name){
	string path;
	if(!getProgramPath(name,path)) return false;
	return STPlot_aux_fileExists(path);
}


void MultiPlot::getOutputCommands(	double suggestedWidth, 
									double suggestedHeight, 
									long sessionNumber, 
									string &outputCommand, 
									string &closingShellCommand, 
									string &terminalCommand, 
									string &terminalTitleCommand, 
									double &resizeTicsByFactor) const{
	const bool epstopdf	= program_exists("epstopdf");
	const bool eps2eps	= program_exists("eps2eps");
	const bool ps2pdf	= program_exists("ps2pdf");
	const bool pstopdf	= program_exists("pstopdf");
	ostringstream terminalCommandSS;
	double lineDownscaling, fontDownscaling, width, height;
	
	// make sure total plot dimensions do not exceed certain bounds (rendering problems appear otherwise)
	// if we need to cut down the size, we should also downscale fonts, lines, dashes and tics
	if(outputType==PlotOutputTypeFilePDF || outputType==PlotOutputTypeFileEPS){
		width = min(suggestedWidth, ST_DEFAULT_MAX_EPS_WIDTH);
		height = min(suggestedHeight, ST_DEFAULT_MAX_EPS_HEIGHT);
		// fix aspect ratios
		if(width/suggestedWidth < height/suggestedHeight){
			height = suggestedHeight * (width/suggestedWidth);
		}else{
			width = suggestedWidth * (height/suggestedHeight);
		}
		lineDownscaling = fontDownscaling = min(width/suggestedWidth, height/suggestedHeight);
		resizeTicsByFactor = lineDownscaling;
	}else{
		width = suggestedWidth;
		height = suggestedHeight;
		lineDownscaling = fontDownscaling = resizeTicsByFactor = 1;
	}
		
	if(outputType==PlotOutputTypeTerminal){
		outputCommand = "unset output";
		if(terminalName=="aqua") outputCommand += " fontsize " + STPlot_makeString(ST_DEFAULT_AQUATERM_FONTSIZE);
		closingShellCommand = "";
		terminalCommandSS << terminalName << " " << sessionNumber << " size " << width << ", " << height;
		terminalTitleCommand = "title '" + terminalTitle + "'";
		
	}else if(outputType==PlotOutputTypeFilePDF && (epstopdf || pstopdf || ps2pdf)){
		terminalCommandSS << "postscript eps enhanced color dashed size " << width << "in, " << height << "in fontscale "  << fontDownscaling << " lw " << lineDownscaling << " dashlength " << lineDownscaling;
		terminalTitleCommand = "";
		if(hasPlotWithEPSBlurProblems() && ST_CORRECT_EPS_BLURS && eps2eps){
			string temp_eps1 = filePath + "_temp1.eps";
			string temp_eps2 = filePath + "_temp2.eps";
			closingShellCommand = 	string(ST_GNUPLOT_COMMENT_PREFIX) + "   This command is needed to correct EPS rendering errors with some viewers\n"
									+ "!eps2eps -dNOCACHE '"+temp_eps1+"' '" + temp_eps2 + "'\n"
									ST_GNUPLOT_COMMENT_PREFIX + "   Transform EPS to PDF\n"
									+ (epstopdf ? "!epstopdf --filter < '" + temp_eps2 + "' > " : (ps2pdf ? "!ps2pdf -dEPSCrop '"+temp_eps2+"' " : "!pstopdf '" + temp_eps2 + "' -o ")) + "'" + filePath + "'\n"
									ST_GNUPLOT_COMMENT_PREFIX + "   Remove temporary EPS files\n"
									+ "!rm '" + temp_eps1 + "' '" + temp_eps2 + "'\n";

			outputCommand = "set output '" + temp_eps1 + "'";
		}else{
			outputCommand = string("set output '| ") + (epstopdf ? "epstopdf --filter > " : (ps2pdf ? "ps2pdf -dEPSCrop - " : "pstopdf -i -o ")) +"\"" + filePath + "\"'";
			closingShellCommand = "";
		}
		
	}else if((outputType==PlotOutputTypeFileEPS) || (outputType==PlotOutputTypeFilePDF && (!epstopdf) && (!pstopdf) && (!ps2pdf))){
		string correct_filePath = filePath;
		if(outputType==PlotOutputTypeFilePDF){
			// requested to create PDF, but cannot do it, so change file extension to ps
			STPlot_aux_correctFileExtensionFromTo(correct_filePath,"pdf","ps");
		}
		terminalCommandSS <<"postscript eps enhanced color dashed size " << width << "in, " << height << "in fontscale " << fontDownscaling << " lw " << lineDownscaling << " dashlength " << lineDownscaling;
		terminalTitleCommand = "";
		if(hasPlotWithEPSBlurProblems() && ST_CORRECT_EPS_BLURS && eps2eps){
			string temp_eps = correct_filePath + "_temp.eps";
			closingShellCommand = 	string(ST_GNUPLOT_COMMENT_PREFIX) + "   This command is needed to correct EPS rendering errors with some viewers\n"
									+ "!eps2eps -dNOCACHE '" + temp_eps + "' '" + correct_filePath + "'\n"
									ST_GNUPLOT_COMMENT_PREFIX + "   Remove temporary EPS file\n"
									+ "!rm '" + temp_eps + "'\n";

			outputCommand = "set output '" + temp_eps + "'";
		}else{
			outputCommand = "set output '" + correct_filePath + "'";
			closingShellCommand = "";
		}
		
	}else if(outputType==PlotOutputTypeFileSVG){
		terminalCommandSS 		<< "svg size " << width << ", " << height << " enhanced dashed fsize " + STPlot_makeString(ST_DEFAULT_SVG_FONTSIZE);
		terminalTitleCommand 	= "";
		closingShellCommand 	= "unset output";
		outputCommand 			= "set output '" + filePath + "'";
	}
	terminalCommand = terminalCommandSS.str();
	return;
}

bool MultiPlot::outputIsToFile() const{
	return (outputType==PlotOutputTypeFilePDF) || (outputType==PlotOutputTypeFileEPS) || (outputType==PlotOutputTypeFileSVG);
}



bool MultiPlot::outputIsPS() const{
	return (outputType==PlotOutputTypeFilePDF) || (outputType==PlotOutputTypeFileEPS);
}


template<class SCONTAINER> 
void explodeString(SCONTAINER &parts, const string &haystack, const string &separator, int maxPartsCount){
	parts.clear();
	if(haystack=="") return;
	int previousPos = 0, nextPos, separatorLen = separator.length();
	while(((nextPos = haystack.find(separator, previousPos)) != string::npos) && ((maxPartsCount>1) || (maxPartsCount<0))){
		parts.push_back(haystack.substr(previousPos, nextPos - previousPos));
		previousPos = nextPos + separatorLen;
		-- maxPartsCount;
	}
	parts.push_back(haystack.substr(previousPos));
}


void MultiPlot::renderFootnote(	string &openingCommand, 
								string &closingCommand, 
								double &footnoteScreenHeight,
								double plotHeightWithoutFootnote, 
								double &plotHeightWithFootnote,
								long &nextAvailLabelTag){
	const string formattedFootnote = (escapeUnderscoresInFootnote ? STPlot_aux_get_replaceAll(footnote,"_","\\\\_") : footnote);
	if(formattedFootnote == ""){
		footnoteScreenHeight = 0;
		openingCommand == "";
		closingCommand == "";
		plotHeightWithFootnote = plotHeightWithoutFootnote;
		return;
	}


	vector<string> footnoteLines;
	explodeString(footnoteLines, formattedFootnote, "\n", -1);
	double footnoteHeight = (formattedFootnote == "" ? 0 : getLabelHeight(formattedFootnote));
	footnoteScreenHeight = footnoteHeight / (plotHeightWithoutFootnote + footnoteHeight);
	plotHeightWithFootnote = plotHeightWithoutFootnote + footnoteHeight;
	double screenHeightPerLine = footnoteScreenHeight/footnoteLines.size();
	ostringstream openingCommandStream, closingCommandStream;
	for(unsigned int i=0; i<footnoteLines.size(); ++i){
		openingCommandStream 	<< (i==0 ? "" : "\n") << "set label " << nextAvailLabelTag 
								<< " \"" << (interpretSymbols ? replaceSymbols(footnoteLines[i]) : footnoteLines[i])
								<< "\" at screen 0.01, " << (footnoteScreenHeight - i*screenHeightPerLine);
		closingCommandStream	<< (i==0 ? "" : "\n") << "unset label " << nextAvailLabelTag;
		++nextAvailLabelTag;
	}
	openingCommand = openingCommandStream.str();
	closingCommand = closingCommandStream.str();
}




//plot multiple Plots
PlotError MultiPlot::plot(long sessionNumber){
	unsigned int 	rows, columns;
	unsigned int	row, column;
	double			originX, originY;
	PlotError		error;
	string			commandJustAfterFirstPlot;
	string 			footnoteOpeningCommand, footnoteClosingCommand;
	double 			footnoteScreenHeight, totalPlotHeightWithFootnote, resizeTicsByFactor;
	long			nextAvailLabelTag = 1;
	string			outputCommand, closingShellCommand, terminalCommand, terminalTitleCommand;
	
	if(plots.empty()) 		return PlotErrorNothingToBeDone;
	
	if(outputIsToFile()){
		pipe.transformToScript();
	}else{
		pipe.transformToPipe();
	}	
	if((plots.size() > 1) || outputIsPS()){ //some problems occur with margins in EPS files in non-multiplot mode
		getFactorization(plots.size(), optimalRatio, minRatio, maxRatio, rows, columns);
		const double totalWidth = columns*plotWidth + (columns+1)*plotHSpacing*plotWidth;
		renderFootnote(	footnoteOpeningCommand, 
						footnoteClosingCommand, 
						footnoteScreenHeight,
						plotHeight * rows + (rows+1)*plotVSpacing*plotHeight, 
						totalPlotHeightWithFootnote,
						nextAvailLabelTag);
		getOutputCommands(totalWidth, totalPlotHeightWithFootnote, sessionNumber, outputCommand, closingShellCommand, terminalCommand, terminalTitleCommand, resizeTicsByFactor);
						
		//prepare plot window and footnote
		pipe.write(	"set term %s %s %s;\nset size 1,1;\n%s;\n", 
					terminalCommand.c_str(),
					terminalTitleCommand.c_str(),
					(font == "" ? "" : (" fname '" + font + "'").c_str()),
					outputCommand.c_str());
		pipe.write("set multiplot\n%s\n",(outputType==PlotOutputTypeFileSVG ? "" : "clear\n")); // for EPS this needs to come after multiplot, but for SVGs that generates an error
		//pipe.write("set lmargin %i; set rmargin %i\n\n", ST_DEFAULT_MULTIPLOT_LMARGIN, ST_DEFAULT_MULTIPLOT_RMARGIN);
		if(footnote != ""){
			pipe.write("%s%s\n%s\n\n", ST_GNUPLOT_COMMENT_PREFIX, "Footnote", footnoteOpeningCommand.c_str());
			commandJustAfterFirstPlot = footnoteClosingCommand;
		}else{
			footnoteScreenHeight = 0;
		}
		
		const double relativePlotHSpacing 	= plotHSpacing*plotWidth/totalWidth;
		const double relativePlotVSpacing 	= plotVSpacing*plotHeight/totalPlotHeightWithFootnote;
		const double relativePlotWidth 		= plotWidth/totalWidth;
		const double relativePlotHeight 	= plotHeight/totalPlotHeightWithFootnote;
		
		for(unsigned int i=0; i<plots.size(); ++i){
			column 	= i % columns;
			row 	= (i - column)/columns;
			originX = column * relativePlotWidth + (column+1)*relativePlotHSpacing;
			originY	= 1 - (row+1) * (relativePlotHeight + relativePlotVSpacing);
			//pipe.write("set size %e, %e\nset origin %e, %e\n", relativePlotWidth, relativePlotHeight, originX, originY);
			error = plots[i].writeCommandToPipe(pipe, 
												originX, 
												originY, 
												relativePlotWidth, 
												relativePlotHeight,
												(onlyResizePlotAreaInResizedPlots ? defaultPlotWidth/plotWidth : 1.0),
												(onlyResizePlotAreaInResizedPlots ? defaultPlotHeight/plotHeight : 1.0), 
												resizeTicsByFactor,
												font,
												interpretSymbols,
												0, 
												(tagPlots && (plots.size()>1) ? generatePlotTag(i) : ""),
												escapeUnderscoresInTitle,
												escapeUnderscoresInTicsAndLabels,
												(i == 0 ? commandJustAfterFirstPlot : ""),
												nextAvailLabelTag);
			if(error != PlotErrorNone){
				pipe.transformToPipe();
				return error;
			}
		}
		
		pipe.write("\nunset multiplot;\n");
		pipe.write("reset;\n\n\n");
	
	}else{
		//Only a single plot. Do not use multiplot modus, since it needlessly prevents terminal-user interaction
		const double totalWidth = plotWidth + 2*plotHSpacing*plotWidth;
		renderFootnote(	footnoteOpeningCommand, 
						footnoteClosingCommand, 
						footnoteScreenHeight,
						plotHeight+2*plotVSpacing*plotHeight, 
						totalPlotHeightWithFootnote,
						nextAvailLabelTag);
		getOutputCommands(totalWidth, totalPlotHeightWithFootnote, sessionNumber, outputCommand, closingShellCommand, terminalCommand, terminalTitleCommand, resizeTicsByFactor);
		
		pipe.write(	"set term %s %s %s;\n%s;\n%s\n", 
					terminalCommand.c_str(), 
					terminalTitleCommand.c_str(),
					(font == "" ? "" : (" fname '" + font + "'").c_str()),
					outputCommand.c_str(),
					(outputType==PlotOutputTypeFileSVG ? "" : "clear\n"));

		//footnote
		if(footnote != ""){
			pipe.write("%s%s\n%s\n\n", ST_GNUPLOT_COMMENT_PREFIX, "Footnote", footnoteOpeningCommand.c_str());
			commandJustAfterFirstPlot = footnoteClosingCommand;
		}else{
			footnoteScreenHeight = 0;
		}
		

		error = plots[0].writeCommandToPipe(pipe, 
											plotHSpacing*plotWidth/totalWidth, 
											footnoteScreenHeight + plotVSpacing*plotHeight/totalPlotHeightWithFootnote, 
											plotWidth/totalWidth, 
											plotHeight/totalPlotHeightWithFootnote,
											(onlyResizePlotAreaInResizedPlots ? defaultPlotWidth/plotWidth : 1.0), 
											(onlyResizePlotAreaInResizedPlots ? defaultPlotHeight/plotHeight : 1.0), 
											resizeTicsByFactor,
											font,
											interpretSymbols,
											0,
											"",
											escapeUnderscoresInTitle,
											escapeUnderscoresInTicsAndLabels,
											commandJustAfterFirstPlot,
											nextAvailLabelTag);
		pipe.write("\nreset;\n\n\n");
		if(error != PlotErrorNone){
			pipe.transformToPipe();
			return error;
		}
	}
	
	if(closingShellCommand!=""){
		pipe.write("\n%s Closing shell commands\n%s\n", ST_GNUPLOT_COMMENT_PREFIX, closingShellCommand.c_str());
	}
	
	if(outputIsToFile()){
		bool er = !pipe.executeScript();
		pipe.transformToPipe();
		if(er){ return (foundGnuplot() ? PlotErrorUnknownError : PlotErrorPlotterNotFound); }
	}
	return PlotErrorNone;
}


PlotError MultiPlot::plot(const Plot &plot, const string &title){
	MultiPlot mp(plot, title);
	return mp.plot();	
}


PlotError MultiPlot::plot(const Plot &plot, const string &title, const string &gnuplotPath){
	MultiPlot mp(plot, title, gnuplotPath);
	return mp.plot();	
}

	

void MultiPlot::setLogStream(ostream *logStream){
	pipe.setLogStream(logStream);
}


void MultiPlot::setAspectRatio(double _optimalRatio){
	optimalRatio 	= _optimalRatio;
	minRatio 		= 0;
	maxRatio 		= numeric_limits<double>::infinity();
}

void MultiPlot::setAspectRatio(double _optimalRatio, double _minRatio, double _maxRatio){
	optimalRatio	= _optimalRatio;
	minRatio		= min(_optimalRatio * (1.0 - ST_DEFAULT_ASPECT_RATIO_RELATIVE_TOLLERANCE), _minRatio);
	maxRatio		= max(_optimalRatio * (1.0 + ST_DEFAULT_ASPECT_RATIO_RELATIVE_TOLLERANCE), _maxRatio);
}




void MultiPlot::getFactorization(unsigned int value, double ratio, unsigned int &r, unsigned int &c){
	if(value == 0){
		r = c = 0;
		return;
	}
	double r0 = sqrt(value/ratio);
	double c0 = sqrt(value * ratio);
	
	unsigned int rt=1, ct=value;
	for(int factor = (int)max(1.0,floor(r0)); factor >= 1; --factor){
		if((value % factor) == 0){
			r = factor;
			c = value / factor;
			break;
		}
	}
	double dist = fabs(r - r0) + fabs(c - c0);
	
	for(int factor = (int)max(1.0,floor(c)); factor >= 1; --factor){
		if((value % factor) == 0){
			rt = factor;
			ct = value / factor;
			break;
		}
	}
	if((fabs(rt - r0) + fabs(ct - c0)) < dist){
		r = rt; 
		c = ct;
	}
}


void MultiPlot::getFactorization(unsigned int value, double optimalRatio, double minRatio, double maxRatio, unsigned int &r, unsigned int &c){
	do{
		getFactorization(value, optimalRatio, r, c);
		++ value;
	}while((c > r * maxRatio) || (c < r * minRatio));
}


unsigned int MultiPlot::countNeedle(const string &haystack, const string &needle){
	int pos = 0, count = 0; 
	int needleLen = needle.length();
	while((pos = haystack.find(needle, pos))!= string::npos){
		pos += needleLen;
		++count;
	}
	return count;
}


double MultiPlot::getLabelHeight(const string &label){
	return (1 + countNeedle(label, "\n")) * (outputIsPS() ? ST_DEFAULT_LABEL_HEIGHT_EPS : (outputType==PlotOutputTypeFileSVG ? ST_DEFAULT_LABEL_HEIGHT_SVG : ST_DEFAULT_LABEL_HEIGHT_TERM));
}



/*
void MultiPlot::copyTo(MultiPlot &MP) const{
	MP.plots 				= plots;
	MP.footnote 			= footnote;
	MP.interpretSymbols= interpretSymbols;
	MP.optimalRatio 		= optimalRatio;
	MP.minRatio 			= minRatio;
	MP.maxRatio 			= maxRatio;
	MP.filePath			 	= filePath;
	MP.terminalName 		= terminalName;
	MP.outputType 			= outputType;
	MP.terminalTitle 		= terminalTitle;
	MP.plotWidth 			= plotWidth;
	MP.plotHeight 			= plotHeight;
	MP.plotHSpacing			= plotHSpacing;
	MP.plotVSpacing			= plotVSpacing;
}
*/

unsigned int MultiPlot::highestPlotDimension() const{
	unsigned int hpd = 0;
	for(unsigned int n=0; n<plots.size(); ++n){
		hpd = max(hpd, plots[n].highestPlotDimension());
	}
	return hpd;
}



set<string> MultiPlot::sourceFiles() const{
	set<string> files;
	long p,s;
	for(p=0; p<plots.size(); ++p){
		for(s=0; s<plots[p].sources.size(); ++s){
			if(plots[p].sources[s].usesSourceFile()){ files.insert(plots[p].sources[s].filePath); }
		}
	}
	return files;
}


string	MultiPlot::generatePlotTag(long number) const{
	if(number==0) return "a";
	string tag="";
	const int period = (1+'z'-'a');
	while(number>0){
		tag = string(1,'a' + number % period) + tag;
		number = (number/period)  -1;
	}
	if(number==0) tag = "a"+tag;
	return tag;
}


bool MultiPlot::hasPlotWithEPSBlurProblems() const{
	for(long i=0; i<plots.size(); ++i){
		if(plots[i].hasPlotWithEPSBlurProblems()) return true;
	}
	return false;
}




void MultiPlot::setLineWidthForAllCurrentPlots(double width){
	for(long p=0; p<plots.size(); ++p){
		plots[p].setLineWidthForAllCurrentSources(width);
	}
}


void MultiPlot::setLineWidthAutoForAllCurrentPlots(){
	for(long p=0; p<plots.size(); ++p){
		plots[p].setLineWidthAutoForAllCurrentSources();
	}
}


void MultiPlot::setLineTypeForAllCurrentPlots(int type){
	for(long p=0; p<plots.size(); ++p){
		plots[p].setLineTypeForAllCurrentSources(type);
	}
}


void MultiPlot::setLineTypeAutoForAllCurrentPlots(){
	for(long p=0; p<plots.size(); ++p){
		plots[p].setLineTypeAutoForAllCurrentSources();
	}
}


void MultiPlot::setPointSizeForAllCurrentPlots(double size){
	for(long p=0; p<plots.size(); ++p){
		plots[p].setPointSizeForAllCurrentSources(size);
	}
}


void MultiPlot::setPointSizeAutoForAllCurrentPlots(){
	for(long p=0; p<plots.size(); ++p){
		plots[p].setPointSizeAutoForAllCurrentSources();
	}
}


void MultiPlot::setPointTypeForAllCurrentPlots(int type){
	for(long p=0; p<plots.size(); ++p){
		plots[p].setPointTypeForAllCurrentSources(type);
	}
}


void MultiPlot::setPointTypeAutoForAllCurrentPlots(){
	for(long p=0; p<plots.size(); ++p){
		plots[p].setPointTypeAutoForAllCurrentSources();
	}
}



void MultiPlot::setColorForAllCurrentPlots(const Color &color){
	for(long p=0; p<plots.size(); ++p){
		for(long s=0; s<plots[p].sources.size(); ++s){
			plots[p].sources[s].setColor(color);
		}
	}
}



#pragma mark -
#pragma mark High-level plot creation
#pragma mark -


PlotSource vectorField2PlotSource(	bool (*vectorField)(double, double, double&, double&),
									string filePath, 
									string title, 
									unsigned int resolutionX, 
									unsigned int resolutionY, 
									double minX, 
									double maxX, 
									double minY, 
									double maxY, 
									VectorColoring coloring,
									double meanVectorLength,
									PlotError &error)
{
	ofstream fout;
	Color color;
	double r,g,b;
	double x,y,vx,vy,xStep,yStep, norm, vectorScale, relativeNorm;
	double maxVx, maxVy, minVx, minVy, minNorm, maxNorm, avWNorm;
	fout.open(filePath.c_str(), ios::out);
	if(!fout.good()){
		error = PlotErrorFileWriteError;
		return PlotSource(PlotTypeNone, "", "", 0);
	}
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Vector field sampled " << (coloring == VectorColoringNone ? "" : " and colored") << " on a " << resolutionX << " x " << resolutionY << " grid, within the domain [" << minX << ", " << maxX << "] x [" << minY << ", " << maxY << "]\n";
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Data columns are: x y vx vy" << (coloring == VectorColoringNone ? "" : " r g b") << endl;
	error = PlotErrorNone;
	double Dx = (maxX - minX);
	double Dy = (maxY - minY);
	xStep = Dx / resolutionX;
	yStep = Dy / resolutionY;
	
	//get extreme values for correct scaling and coloring
	vectorField(minX, minY, minVx, minVy);
	maxVx = minVx; maxVy = minVy;
	minNorm = maxNorm = pow(minVx,2) + pow(minVy,2);
	avWNorm = sqrt(pow(minVx/Dx,2) + pow(minVy/Dy,2));
	long count=1;
	for(x = minX; x <= maxX; x += xStep){
		for(y = minY; y <= maxY; y += yStep){
			if(!vectorField(x, y, vx, vy)) continue;
			norm  = pow(vx,2) + pow(vy,2);
			minVx = min(minVx, vx);
			minVy = min(minVy, vy);
			maxVx = max(maxVx, vx);
			maxVy = max(maxVy, vy);
			minNorm = min(minNorm, norm);
			maxNorm = max(maxNorm, norm);
			avWNorm = avWNorm + sqrt(pow(vx/Dx,2) + pow(vy/Dy,2));
			++count;
		 }
	}
	minNorm = sqrt(minNorm); maxNorm = sqrt(maxNorm); avWNorm /= count;
	vectorScale = (meanVectorLength<0 ? 1.0: meanVectorLength/avWNorm);
	
	if(coloring == VectorColoringNone){
		//no coloring, just write sampled vectorField to file
		for(x = minX; x <= maxX; x += xStep){
			for(y = minY; y <= maxY; y += yStep){
				 vectorField(x, y, vx, vy);
				 fout << x << " " << y << " " << vectorScale * vx << " " << vectorScale * vy << "\n";
			 }
		}
		fout << endl;
		fout.close();
		return PlotSource(PlotType2DVectors, title, filePath, 0, 1, 2, 3, 4);
	
	}else{
		//coloring enabled, write color-info with vectorField to file
		
		//write sampled vectorField + colors to file
		for(x = minX; x <= maxX; x += xStep){
			for(y = minY; y <= maxY; y += yStep){
				if(!vectorField(x, y, vx, vy)) continue;
				norm = sqrt(pow(vx, 2) + pow(vy,2));
				color = vector2RGBColor(	(maxVx == minVx ? 1.0 : (vx - minVx)/(maxVx - minVx)), 
											(maxVy == minVy ? 1.0 : (vy - minVy)/(maxVy - minVy)), 
											0.0, 
											(maxNorm == minNorm ? 1.0 : (norm - minNorm)/(maxNorm - minNorm)),
											coloring);
				color.getRGB(r,g,b);
				fout << x << " " << y << " " << vectorScale * vx << " " << vectorScale * vy << " " << r << " " << g << " " << b << "\n";
			 }
		}
		fout << endl;
		fout.close();
		
		return PlotSource(PlotType2DColoredVectors, title, filePath, 0, 1, 2, 3, 4, 5, 6, 7);
	}
	
}





PlotSource vectorField2PlotSource(	bool (*vectorField)(double, double, double, double&, double&, double&),
									string filePath, 
									string title, 
									unsigned int resolutionX, 
									unsigned int resolutionY, 
									unsigned int resolutionZ,
									double minX, 
									double maxX, 
									double minY, 
									double maxY,
									double minZ,
									double maxZ,
									VectorColoring coloring,
									double meanVectorLength,
									PlotError &error)
{
	ofstream fout;
	Color color;
	double r,g,b;
	double x,y,z,vx,vy,vz,xStep,yStep,zStep, norm, vectorScale, relativeNorm;
	double maxVx, maxVy, maxVz, minVx, minVy, minVz, minNorm, maxNorm, avWNorm;
	fout.open(filePath.c_str(), ios::out);
	if(!fout.good()){
		error = PlotErrorFileWriteError;
		return PlotSource(PlotTypeNone, "", "", 0);
	}
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Vector field sampled " << (coloring == VectorColoringNone ? "" : " and colored") << " on a " << resolutionX << " x " << resolutionY << " x " << resolutionZ << " grid, within the domain [" << minX << ", " << maxX << "] x [" << minY << ", " << maxY << "] x [" << minZ << ", " << maxZ << "]\n";
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Data columns are: x y z vx vy vz" << (coloring == VectorColoringNone ? "" : " r g b") << endl;
	error = PlotErrorNone;
	
	double Dx = (maxX - minX);
	double Dy = (maxY - minY);
	double Dz = (maxZ - minZ);
	xStep = Dx / resolutionX;
	yStep = Dy / resolutionY;
	zStep = Dz / resolutionZ;
	
	//get extreme values for correct scaling and coloring
	vectorField(minX, minY, minZ, minVx, minVy, minVz);
	maxVx = minVx; maxVy = minVy; maxVz = minVz; 
	minNorm = maxNorm = pow(minVx, 2) + pow(minVy, 2) + pow(minVz, 2);
	avWNorm = sqrt(pow(minVx/Dx, 2) + pow(minVy/Dy, 2) + pow(minVz/Dz, 2));
	long count=1;
	for(x = minX; x <= maxX; x += xStep){
		for(y = minY; y <= maxY; y += yStep){
			for(z = minZ; z <= maxZ; z += zStep){
				if(!vectorField(x, y, z, vx, vy, vz)) continue;
				norm = pow(vx, 2) + pow(vy, 2) + pow(vz, 2);
				minVx = min(minVx, vx);
				minVy = min(minVy, vy);
				minVz = min(minVz, vz);
				maxVx = max(maxVx, vx);
				maxVy = max(maxVy, vy);
				maxVz = max(maxVz, vz);
				minNorm = min(minNorm, norm);
				maxNorm = max(maxNorm, norm);
				avWNorm = avWNorm + sqrt(pow(vx/Dx,2) + pow(vy/Dy,2) + pow(vz/Dz,2));
				++count;
			}
		 }
	}
	minNorm = sqrt(minNorm); maxNorm = sqrt(maxNorm); avWNorm /= count;
	vectorScale = (meanVectorLength<0 ? 1.0 : meanVectorLength/avWNorm);
	
	if(coloring == VectorColoringNone){
		//no coloring, just write sampled vectorField to file
		for(x = minX; x <= maxX; x += xStep){
			for(y = minY; y <= maxY; y += yStep){
				for(z = minZ; z <= maxZ; z += zStep){
				 	if(!vectorField(x, y, z, vx, vy, vz)) continue;
				 	fout << x << " " << y << " " << z << " " << vectorScale * vx << " " << vectorScale * vy << " " << vectorScale * vz << "\n";
			 	}
			 }
		}
		fout << endl;
		fout.close();
		return PlotSource(PlotType3DVectors, title, filePath, 0, 1, 2, 3, 4, 5, 6);
	
	}else{
		//coloring enabled, write color-info with vectorField to file
		
		for(x = minX; x <= maxX; x += xStep){
			for(y = minY; y <= maxY; y += yStep){
				for(z = minZ; z <= maxZ; z += zStep){
				 	if(!vectorField(x, y, z, vx, vy, vz)) continue;
				 	norm = pow(vx, 2) + pow(vy, 2) + pow(vz, 2);
					color = vector2RGBColor(	(maxVx == minVx ? 1.0 : (vx - minVx)/(maxVx - minVx)), 
												(maxVy == minVy ? 1.0 : (vy - minVy)/(maxVy - minVy)), 
												(maxVz == minVz ? 1.0 : (vz - minVz)/(maxVz - minVz)), 
												(maxNorm == minNorm ? 1.0 : (norm - minNorm)/(maxNorm - minNorm)), 
												coloring);
					color.getRGB(r,g,b);
					fout << x << " " << y << " " << z << " " << vectorScale * vx << " " << vectorScale * vy << " " << vectorScale * vz << " " << r << " " << g << " " << b << "\n";
				}
			}
		}
		fout << endl;
		fout.close();
		return PlotSource(PlotType3DColoredVectors, title, filePath, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9);
	}
}



PlotSource vectorFieldFlowLines2PlotSource( bool (*vectorField)(double, double, double&, double&),
											string filePath, 
											string title, 
											unsigned int resolutionX, 
											unsigned int resolutionY,
											double timeStep,
											double flowTime,
											IntegrationMethod integrationMethod,
											bool followFlowsOutsideBounds,
											double minX, 
											double maxX, 
											double minY, 
											double maxY,
											PlotError &error)
{
	ofstream fout;
	unsigned int i;
	double x,y,xStep,yStep;
	double t, point[2], V[2], tempPointA[2], tempPointB[2], tempPointC[2], tempVA[2], tempVB[2], tempVC[2];
	fout.open(filePath.c_str(), ios::out);
	if(!fout.good()){
		error = PlotErrorFileWriteError;
		return PlotSource(PlotTypeNone, "", "", 0);
	}
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Vector field flow-lines (trajectories) integrated up to time " << flowTime << ", starting on a " << resolutionX << " x " << resolutionY << " grid, within the domain [" << minX << ", " << maxX << "] x [" << minY << ", " << maxY << "]\n";
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Integration method was " << Plot::integrationMethod2string(integrationMethod) << " using a time-step of " << timeStep << "\n";
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Data columns are: x y\n";
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Different flow-lines are separated by one empty line" << endl;
	error = PlotErrorNone;
	
	xStep = (maxX - minX) / resolutionX;
	yStep = (maxY - minY) / resolutionY;
	
	for(x = (resolutionX == 1 ? (maxX + minX)/2.0 : minX); x <= maxX; x += xStep){
		for(y = (resolutionY == 1 ? (maxY + minY)/2.0 : minY); y <= maxY; y += yStep){
			//follow flow line starting at point (x,y) up to time flowTime
			point[0] = x; point[1] = y;
			switch(integrationMethod){
			case IntegrationMethodForwardEuler:
				for(t = 0.0; t <= flowTime; t += timeStep){
					//write point to file
					fout << point[0] << " " << point[1] << " " << point[2] << "\n";
					//check if point is within bounds
					if(!followFlowsOutsideBounds){
						if((point[0] > maxX) || (point[0] < minX) || (point[1] > maxY) || (point[1] < minY)) break;
					}
					//calculate new point
					if(!vectorField(point[0], point[1], V[0], V[1])) break;
					for(i=0; i<2; ++i) point[i] = point[i] + timeStep * V[i];
				}
				break;
				
			case IntegrationMethodClassicalRungeKutta:
				for(t = 0.0; t <= flowTime; t += timeStep){
					//write point to file
					fout << point[0] << " " << point[1] << " " << point[2] << "\n";
					//check if point is within bounds
					if(!followFlowsOutsideBounds){
						if((point[0] > maxX) || (point[0] < minX) || (point[1] > maxY) || (point[1] < minY)) break;
					}
					//calculate new point
					if(!vectorField(point[0], point[1], V[0], V[1])) break;
					for(i=0; i<2; ++i) tempPointA[i] = point[i] + (timeStep/2.0) * V[i];
					if(!vectorField(tempPointA[0], tempPointA[1], tempVA[0], tempVA[1])) break;
					for(i=0; i<2; ++i) tempPointB[i] = point[i] + (timeStep/2.0) * tempVA[i];
					if(!vectorField(tempPointB[0], tempPointB[1], tempVB[0], tempVB[1])) break;
					for(i=0; i<2; ++i) tempPointC[i] = point[i] + timeStep * tempVB[i];
					if(!vectorField(tempPointC[0], tempPointC[1], tempVC[0], tempVC[1])) break;
					for(i=0; i<2; ++i) point[i] = point[i] + (timeStep/6.0) * (V[i] + 2 * (tempVA[i] + tempVB[i]) + tempVC[i]);
				}
				break;
			}
			fout << endl;
		
		 }
	}
	fout << endl;
	fout.close();
	return PlotSource(PlotType2DCurve, title, filePath, 0, 1, 2);

}



PlotSource vectorFieldFlowLines2PlotSource( bool (*vectorField)(double, double, double, double&, double&, double&),
											string filePath, 
											string title, 
											unsigned int resolutionX, 
											unsigned int resolutionY,
											unsigned int resolutionZ,
											double timeStep,
											double flowTime,
											IntegrationMethod integrationMethod,
											bool followFlowsOutsideBounds,
											double minX, 
											double maxX, 
											double minY, 
											double maxY, 
											double minZ,
											double maxZ,
											PlotError &error)
{
	ofstream fout;
	unsigned int i;
	double x,y,z,xStep,yStep,zStep;
	double t, point[3], V[3], tempPointA[3], tempPointB[3], tempPointC[3], tempVA[3], tempVB[3], tempVC[3];
	fout.open(filePath.c_str(), ios::out);
	if(!fout.good()){
		error = PlotErrorFileWriteError;
		return PlotSource(PlotTypeNone, "", "", 0);
	}
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Vector field flow-lines (trajectories) integrated up to time " << flowTime << ", starting on a " << resolutionX << " x " << resolutionY << " x " << resolutionZ << " grid, within the domain [" << minX << ", " << maxX << "] x [" << minY << ", " << maxY << "] x [" << minZ << ", " << maxZ << "]\n";
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Integration method was " << Plot::integrationMethod2string(integrationMethod) << " using a time-step of " << timeStep << "\n";
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Data columns are: x y z\n";
	fout << ST_GNUPLOT_COMMENT_PREFIX << "Different flow-lines are separated by one empty line" << endl;
	error = PlotErrorNone;
	
	xStep = (maxX - minX) / resolutionX;
	yStep = (maxY - minY) / resolutionY;
	zStep = (maxZ - minZ) / resolutionZ;
	
	for(x = (resolutionX == 1 ? (maxX + minX)/2.0 : minX); x <= maxX; x += xStep){
		for(y = (resolutionY == 1 ? (maxY + minY)/2.0 : minY); y <= maxY; y += yStep){
			for(z = (resolutionZ == 1 ? (maxZ + minZ)/2.0 : minZ); z <= maxZ; z += zStep){
				//follow flow line starting at point (x,y,z) up to time flowTime
				point[0] = x; point[1] = y; point[2] = z;
				switch(integrationMethod){
				case IntegrationMethodForwardEuler:
					for(t = 0.0; t <= flowTime; t += timeStep){
						//write point to file
						fout << point[0] << " " << point[1] << " " << point[2] << "\n";
						//check if point is within bounds
						if(!followFlowsOutsideBounds){
							if((point[0] > maxX) || (point[0] < minX) || (point[1] > maxY) || (point[1] < minY) || (point[2] > maxZ) || (point[2] < minZ)) break;
						}
						//calculate new point
						if(!vectorField(point[0], point[1], point[2], V[0], V[1], V[2])) break;
						for(i=0; i<3; ++i) point[i] = point[i] + timeStep * V[i];
					}
					break;
					
				case IntegrationMethodClassicalRungeKutta:
					for(t = 0.0; t <= flowTime; t += timeStep){
						//write point to file
						fout << point[0] << " " << point[1] << " " << point[2] << "\n";
						//check if point is within bounds
						if(!followFlowsOutsideBounds){
							if((point[0] > maxX) || (point[0] < minX) || (point[1] > maxY) || (point[1] < minY) || (point[2] > maxZ) || (point[2] < minZ)) break;
						}
						//calculate new point
						if(!vectorField(point[0], point[1], point[2], V[0], V[1], V[2])) break;
						for(i=0; i<3; ++i) tempPointA[i] = point[i] + (timeStep/2.0) * V[i];
						if(!vectorField(tempPointA[0], tempPointA[1], tempPointA[2], tempVA[0], tempVA[1], tempVA[2])) break;
						for(i=0; i<3; ++i) tempPointB[i] = point[i] + (timeStep/2.0) * tempVA[i];
						if(!vectorField(tempPointB[0], tempPointB[1], tempPointB[2], tempVB[0], tempVB[1], tempVB[2])) break;
						for(i=0; i<3; ++i) tempPointC[i] = point[i] + timeStep * tempVB[i];
						if(!vectorField(tempPointC[0], tempPointC[1], tempPointC[2], tempVC[0], tempVC[1], tempVC[2])) break;
						for(i=0; i<3; ++i) point[i] = point[i] + (timeStep/6.0) * (V[i] + 2 * (tempVA[i] + tempVB[i]) + tempVC[i]);
					}
					break;
				}
				fout << endl;
				
			}
		 }
	}
	fout << endl;
	fout.close();
	return PlotSource(PlotType3DCurve, title, filePath, 0, 1, 2, 3);

}


PlotSource complexMapping2PlotSource( void (*mapping)(double, double, double&, double&), 
										string filePath,
										string title,
										unsigned int resolutionR,
										unsigned int resolutionI,
										double minR, 
										double maxR, 
										double minI, 
										double maxI,
										double zeroRadius,
										double inftyRadius,
										PlotError &error)
{
 	ofstream fout;
	double r,g,b;
	double xr,xi,yr,yi,rStep,iStep;
	fout.open(filePath.c_str(), ios::out);
	if(!fout.good()){
		error = PlotErrorFileWriteError;
		return PlotSource(PlotTypeNone, "", "", 0);
	}
	fout 	<< ST_GNUPLOT_COMMENT_PREFIX << "Complex mapping C -> C sampled on a " << resolutionR << " x " << resolutionI << " grid, within the domain [" << minR << ", " << maxR << "] x [" << minI << ", " << maxI << "]\n"
	 		<< ST_GNUPLOT_COMMENT_PREFIX << "Sampling grid given as matrix, whose entries are color numbers, corresponding to the complex value of the function at that grid point.\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << "Complex to color transformation is done using the HSV-cylinder color model. 0 is white and INFTY is black.\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << "Upper left matrix entry is point (" << minR << ", " << minI << ")\n";
	error = PlotErrorNone;
	rStep = (maxR - minR) / resolutionR;
	iStep = (maxI - minI) / resolutionI;
		
	for(xi = minI; xi <= maxI; xi += iStep){
		for(xr = minR; xr <= maxR; xr += rStep){
			 mapping(xr,xi,yr,yi);
			 complex2RGB(yr,yi,zeroRadius,inftyRadius,r,g,b);
			 fout << RGB2ColorNumber(r,g,b) << " ";
		 }
		 fout << "\n";
	}
	fout << endl;
	fout.close();
	
	//create plot source
	PlotSource source = PlotSource(PlotType2DColorImage, title, filePath, 0);
	source.minX = minR;
	source.maxX = maxR;
	source.minY = minI;
	source.maxY = maxI;
	source.xCount = resolutionR;
	source.yCount = resolutionI;
	return source;
}





Plot HeatMapHClustered2Plot(	const vector<double>	&data,		// (INPUT) data matrix in row-major format
								long					Nrows,
								long					Ncolumns,
								const string			&Rlegend,
								const string			&Clegend,
								const string			&ValueLegend,
								const vector<string>	&Rnames,
								const vector<string>	&Cnames,
								const string			&filePath,
								ofstream				&fout,	
								long					&currentDataIndex,
								const string			&title,
								const string			&fileComment,
								double					columnNamesRotation,
								HClusterWhat			clusterWhat,
								HClusterOrdering		Rordering,
								HClusterOrdering		Cordering,
								HClusterMetric			metric,
								HClusterLinkage			linkage,
								bool					uniformize,
								double					&recommendedPlotWidthCorrection,
								double					&recommendedPlotHeightCorrection,
								double					&recommendedRelativePlotHSpacing,
								double					&recommendedRelativePlotVSpacing,
								PlotError 				&error){
	long r,r2,c,c2,n;
	bool foutWasOpen = fout.is_open();
	if(!foutWasOpen){
		fout.open(filePath.c_str(), ios::app);
	}
	if(!fout.good()){
		error = PlotErrorFileWriteError;
		return Plot();
	}
	error = PlotErrorNone;
	string nnFileComment = STPlot_aux_get_replaceAll(fileComment, "\n", (string(ST_GNUPLOT_COMMENT_PREFIX) + "\n").c_str());
	
	const bool clusterRows 		= (clusterWhat==HClusterRows || clusterWhat==HClusterRowsAndColumns || clusterWhat==HClusterRowsAndLockColumns2Rows);
	const bool clusterColumns 	= (clusterWhat==HClusterColumns || clusterWhat==HClusterRowsAndColumns || clusterWhat==HClusterColumnsAndLockRows2Columns);

	long Nmax = max((clusterRows ? Nrows : 0), (clusterColumns ? Ncolumns : 0));
	long *orderR = new long[Nrows];
	long *orderC = new long[Ncolumns];
	long *merge1 = new long[Nmax];
	long *merge2 = new long[Nmax];
	double *distances = new double[Nmax*Nmax];
	bool needDistancesBackup = (((Rordering == HClusterOrderingByDistance) && clusterRows) || ((Cordering == HClusterOrderingByDistance) && clusterColumns));
	double *distances2 = NULL;
	if(needDistancesBackup){ distances2 = new double[Nmax*Nmax]; }
	
	//cluster rows if needed
	if(clusterRows){
		//calculate distance matrix (upper-triangular)
		for(r=0; r<Nrows; ++r){
			for(r2=r+1; r2<Nrows; ++r2){
				distances[r*Nrows + r2] = STPlot_aux_distanceBetweenDataRows(&data[0], r, r2, Ncolumns, metric);
			}
		}
		if(Rordering == HClusterOrderingByDistance){
			//make backup of distance matrix because it will be modified by HClusterAgglomerative()
			for(r=0; r<Nrows; ++r){
				for(r2=r+1; r2<Nrows; ++r2){ 
					distances2[r*Nrows + r2] = distances[r*Nrows + r2];
				}
			}
		}
		HClusterAgglomerative(Nrows, distances, linkage, merge1, merge2, (Rordering == HClusterOrderingArbitrary ? orderR : NULL));
		if(Rordering == HClusterOrderingCentroid){
			orderHClusterRowsCentroid(&data[0], Nrows, Ncolumns, false, metric, merge1, merge2, orderR);
		}else if(Rordering == HClusterOrderingAscending){
			orderHClusterAscending(Nrows, merge1, merge2, orderR);
		}else if(Rordering == HClusterOrderingByDistance){
			orderHClusterByDistance(distances2, Nrows, merge1, merge2, orderR);
		}
	}else if(clusterWhat!=HClusterColumnsAndLockRows2Columns){
		for(r=0; r<Nrows; ++r){ orderR[r] = r; }
	}
	
	//cluster columns if needed
	if(clusterColumns){
		//calculate distance matrix (upper-triangular)
		for(c=0; c<Ncolumns; ++c){
			for(c2=c+1; c2<Ncolumns; ++c2){
				distances[c*Ncolumns + c2] = STPlot_aux_distanceBetweenDataColumns(&data[0], c, c2, Nrows, Ncolumns, metric);
			}
		}
		if(Cordering == HClusterOrderingByDistance){
			//make backup of distance matrix because it will be modified by HClusterAgglomerative()
			for(c=0; c<Ncolumns; ++c){
				for(c2=c+1; c2<Ncolumns; ++c2){ 
					distances2[c*Ncolumns + c2] = distances[c*Ncolumns + c2];
				}
			}
		}
		HClusterAgglomerative(Ncolumns, distances, linkage, merge1, merge2, (Cordering == HClusterOrderingArbitrary ? orderC : NULL));
		if(Cordering == HClusterOrderingCentroid){
			orderHClusterRowsCentroid(&data[0], Nrows, Ncolumns, true, metric, merge1, merge2, orderC);
		}else if(Cordering == HClusterOrderingAscending){
			orderHClusterAscending(Ncolumns, merge1, merge2, orderC);
		}else if(Cordering == HClusterOrderingByDistance){
			orderHClusterByDistance(distances2, Ncolumns, merge1, merge2, orderC);
		}
	}else if(clusterWhat!=HClusterRowsAndLockColumns2Rows){
		for(c=0; c<Ncolumns; ++c){ orderC[c] = c; }
	}
	
	// lock order of rows/columns to each other if necessary
	if(clusterWhat==HClusterColumnsAndLockRows2Columns){
		for(r=0; r<min(Nrows,Ncolumns); ++r){ orderR[r] = orderC[r]; }
		for(r=Ncolumns; r<Nrows; ++r){ orderR[r] = r; }
	}
	if(clusterWhat==HClusterRowsAndLockColumns2Rows){
		for(c=0; c<min(Ncolumns,Nrows); ++c){ orderC[c] = orderR[c]; }
		for(c=Nrows; c<Ncolumns; ++c){ orderC[c] = c; }
	}

	//write file header
	fout 	<< (currentDataIndex>0 ? "\n\n" : "") << ST_GNUPLOT_COMMENT_PREFIX << " Data index: " << currentDataIndex << "\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << " " << nnFileComment << "\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << (clusterWhat==HClusterNone ? " Heatmap" : " Hierarchically clustered heatmap") << " (" << Nrows << " rows x " << Ncolumns << " columns)\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << " Rows" << (Rlegend=="" ? "" : " ("+Rlegend+")") << " are " << (clusterRows ? "" : " not") << " clustered\n";
	if(clusterRows){
		fout	<< ST_GNUPLOT_COMMENT_PREFIX << " Plot_row  original_row\n";
		for(r=0; r<Nrows; ++r){ fout << ST_GNUPLOT_COMMENT_PREFIX << " " << r << " " << orderR[r] << (orderR[r]<Rnames.size() ? " ("+Rnames[orderR[r]]+")" : "") << "\n"; }
		fout 	<< ST_GNUPLOT_COMMENT_PREFIX << "\n";
	}
	fout	<< ST_GNUPLOT_COMMENT_PREFIX << " Columns" << (Clegend=="" ? "" : " ("+Clegend+")") << " are " << (clusterColumns ? "" : " not") << " clustered\n";
	if(clusterColumns){
		fout	<< ST_GNUPLOT_COMMENT_PREFIX << " Plot_column  original_column\n";
		for(c=0; c<Ncolumns; ++c){ fout << ST_GNUPLOT_COMMENT_PREFIX << " " << c << " " << orderC[c] << (orderC[c]<Cnames.size() ? " ("+Cnames[orderC[c]]+")" : "") << "\n"; }
		fout 	<< ST_GNUPLOT_COMMENT_PREFIX << "\n";
	}
	//write data as matrix
	//if matrix is single-columned and/or single-rowed, add a 2nd duplicate column and/or row
	fout 	<< ST_GNUPLOT_COMMENT_PREFIX << " Data: top-left entry corresponds to bottom-left plot pixel\n";
	if(Nrows==1){ fout << ST_GNUPLOT_COMMENT_PREFIX << " Note: 2nd row is a dummy (duplicate) row for technical purposes and should be ignored\n"; }
	if(Ncolumns==1){ fout << ST_GNUPLOT_COMMENT_PREFIX << " Note: 2nd column is a dummy (duplicate) column for technical purposes and should be ignored\n"; }
	for(r=max(2l,Nrows)-1; r>=0; --r){
		for(c=0; c<max(2l,Ncolumns); ++c){
			fout << data[orderR[min(r,Nrows-1)]*Ncolumns + orderC[min(c,Ncolumns-1)]] << " ";
		}
		fout << "\n";
	}
	fout << endl;

	Tics adjustedTics;
	vector<double> transformedData;
	if(uniformize){
		++currentDataIndex;
		STPlot_transformDataToUniformDistribution(	data,
													0,
													Nrows*Ncolumns-1,
													transformedData,
													adjustedTics);
		//write transformed data as matrix, in similar way as original data
		fout 	<< "\n\n"
				<< ST_GNUPLOT_COMMENT_PREFIX << " Data index: " << currentDataIndex << "\n"
				<< ST_GNUPLOT_COMMENT_PREFIX << " Transformed data for plotting purposes\n"
				<< ST_GNUPLOT_COMMENT_PREFIX << " The order of the values is preserved, but values have been non-linearly transformed to become uniformly distributed on the interval [0,1]\n"
			 	<< ST_GNUPLOT_COMMENT_PREFIX << " Data: top-left entry corresponds to bottom-left plot pixel\n";
		if(Nrows==1){ fout << ST_GNUPLOT_COMMENT_PREFIX << " Note: 2nd row is a dummy (duplicate) row for technical purposes and should be ignored\n"; }
		if(Ncolumns==1){ fout << ST_GNUPLOT_COMMENT_PREFIX << " Note: 2nd column is a dummy (duplicate) column for technical purposes and should be ignored\n"; }
		for(r=max(2l,Nrows)-1; r>=0; --r){
			for(c=0; c<max(2l,Ncolumns); ++c){
				fout << transformedData[orderR[min(r,Nrows-1)]*Ncolumns + orderC[min(c,Ncolumns-1)]] << " ";
			}
			fout << "\n";
		}
		fout << endl;
	}
	
	//create plot
	PlotSource source 	= PlotSource(PlotType2DValueImage, "", filePath, currentDataIndex);
	source.minX 		= 0;
	source.maxX 		= max(2l,Ncolumns)-1;
	source.minY 		= 0;
	source.maxY 		= max(2l,Nrows)-1;
	source.xCount 		= max(2l,Ncolumns);
	source.yCount 		= max(2l,Nrows);
	Plot plot(source, MultiLabel(Clegend,Rlegend), MultiRange(Range(),Range()), "");
	if(uniformize){
		plot.CBTics 	= adjustedTics;
		plot.CBRange 	= Range(0,1);
	}
	++ currentDataIndex;

	
	//print row & column category names (if given)
	plot.xTics = plot.yTics = Tics(0);
	for(r=0; r<Nrows; ++r){
		if(orderR[r]<Rnames.size()){ plot.yTics.add(Rnames[orderR[r]], (Nrows==1 ? 0.5 : Nrows-r-1)); }
	}
	for(c=0; c<Ncolumns; ++c){
		if(orderC[c]<Cnames.size()){ plot.xTics.add(Cnames[orderC[c]], (Ncolumns==1 ? 0.5 : c)); }
	}
	plot.xTics.setRotation(columnNamesRotation);
	plot.title = title;
	plot.CBLabel = ValueLegend;
	if(ValueLegend!="") plot.showColorBox = true;
	
	//recommended plot dimensions (relative to defaults)
	recommendedPlotWidthCorrection 	= max(1.0, Ncolumns*ST_DEFAULT_HEATMAP_CORRECTION_PER_COLUMN);
	recommendedPlotHeightCorrection = max(1.0, Nrows*ST_DEFAULT_HEATMAP_CORRECTION_PER_ROW);
	recommendedRelativePlotHSpacing	= max(ST_DEFAULT_PLOT_H_SPACING, max(10.0, 10+abs(cos(columnNamesRotation*180/3.14))*STPlot_aux_getMaxStringLength(Rnames))*ST_DEFAULT_HEATMAP_CORRECTION_PER_ROW_NAME_CHAR - ST_DEFAULT_RELATIVE_LMARGIN)/recommendedPlotWidthCorrection;
	recommendedRelativePlotVSpacing = max(ST_DEFAULT_PLOT_V_SPACING, max(10.0, 10+abs(sin(columnNamesRotation*180/3.14))*STPlot_aux_getMaxStringLength(Cnames))*ST_DEFAULT_HEATMAP_CORRECTION_PER_COLUMN_NAME_CHAR - ST_DEFAULT_RELATIVE_BMARGIN)/recommendedPlotHeightCorrection;

	delete[] orderR, orderC, merge1, merge2, distances;
	if(needDistancesBackup){ delete[] distances2; }
	if(!foutWasOpen) fout.close();
	return plot;
}






Plot Histogram2Plot(	const vector<string> 	&names,
						const vector<double> 	&values,
						const string			&nameLegend,
						const string			&valueLegend,
						const string			&filePath,	
						long					&dataIndex,	
						const string			&title,
						const string			&fileComment,
						const string			&separator,
						double					namesRotation,
						SortHow					sortHow,
						bool					escapeUnderscoresInNames,
						double					&recommendedPlotWidthCorrection,
						double					&recommendedPlotHeightCorrection,
						double					&recommendedRelativePlotHSpacing,
						double					&recommendedRelativePlotVSpacing,
						long					&histogramIndex,
						PlotError 				&error){
	const long N=names.size();
	if(N==0){ error = PlotErrorNothingToBeDone; return Plot(); }
	ofstream fout;
	fout.open(filePath.c_str(), (dataIndex<0 ? ios::out : ios::app));
	if(!fout.good()){ error = PlotErrorFileWriteError; return Plot(); }
	if(dataIndex>=0) fout << "\n\n" << ST_GNUPLOT_COMMENT_PREFIX << " Data index: " << dataIndex << "\n";
	error = PlotErrorNone;
	string nnFileComment = STPlot_aux_get_replaceAll(fileComment, "\n", (string(ST_GNUPLOT_COMMENT_PREFIX) + "\n").c_str());
	double maxValue=values[0], minValue=values[0];

	// figure out category order
	vector<long> order(N);
	if(sortHow!=SortNone){
		STPlot_qsortIndices(values, order, (sortHow==SortAscending ? true : false));
	}else{
		for(long i=0; i<N; ++i){ order[i] = i; }
	}

	//write file header
	fout 	<< ST_GNUPLOT_COMMENT_PREFIX << " " << nnFileComment << "\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << " Histogram (" << N << " categories)\n";
	if(sortHow!=SortNone) fout << ST_GNUPLOT_COMMENT_PREFIX << " Categories are sorted by " << (sortHow==SortDescending ? "decreasing" : "increasing") << " value\n";		
	
	//write data
	for(long i=0; i<N; ++i){
		fout << (escapeUnderscoresInNames ? STPlot_aux_get_replaceAll(names[order[i]],"_","\\\\_") : names[order[i]]) << (separator=="" ? " " : separator) << values[order[i]] << "\n";
		minValue = min(minValue, values[i]);
		maxValue = max(maxValue, values[i]);
	}
	fout.close();

	//recommended plot dimensions (relative to defaults)
	recommendedPlotWidthCorrection 	= 1.0;
	recommendedPlotHeightCorrection	= max(1.0, N*ST_DEFAULT_HEATMAP_CORRECTION_PER_ROW);
	recommendedRelativePlotHSpacing	= ST_DEFAULT_PLOT_H_SPACING/recommendedPlotWidthCorrection;
	recommendedRelativePlotVSpacing = max(ST_DEFAULT_PLOT_V_SPACING, max(10.0, 10+abs(sin(namesRotation*180/3.14))*STPlot_aux_getMaxStringLength(names))*ST_DEFAULT_HEATMAP_CORRECTION_PER_COLUMN_NAME_CHAR/recommendedPlotHeightCorrection)/recommendedPlotHeightCorrection;
	
	//create empty plot
	Plot plot;
	plot.range = MultiRange(Range(),Range());
	plot.axisLabels = MultiLabel(nameLegend,valueLegend);
	plot.title = title;
	
	// add histogram plot
	plot.sources.push_back(PlotSource(PlotType1DHistogram, "", filePath, (dataIndex<0 ? 0 : dataIndex), 1, 2));	
	if(minValue>=0){ plot.range.yRange = Range(0,maxValue*1.2); }
	else if(maxValue<=0){ plot.range.yRange = Range(minValue*1.2,0); }
	else{ plot.range.yRange = Range(1.2*minValue, 1.2*maxValue); }
	plot.dataSeparator = separator;
	plot.xTics.setRotation(namesRotation);
	plot.sources.back().setColor(Color::getNicePlotColor(0));
	histogramIndex = plot.sources.size()-1;
	
	// if bars are both negative & positive, include a horizontal line at zero
	// this must come after the histogram plot so that xtics are correct
	if((minValue<0) && (maxValue>0)){
		plot.sources.push_back(PlotSource(PlotType2DAnalytic, "", "0", 0));
		plot.sources.back().setColor(Color(0.6,0.6,0.6));	
	}
		
	++dataIndex;
	return plot;
}







Plot StackedFilledCurves2Plot(	const double			data[],
								unsigned long			xcount,
								unsigned long			ycount,
								const vector<Color>		&colors,
								const vector<string>	&names,
								const string			&filePath,
								long					dataIndex,	
								const string			&title,
								const string			&fileComment,
								const string			&xlabel,
								const string			&ylabel,
								bool					relative,
								ColorBlindness			correctForColorBlindness,
								PlotError				&error){
	long i,j,n;
	double cumulative, rescaleAtZero;
	Color color;
	if(ycount==0){
		error = PlotErrorNothingToBeDone;
		return Plot();
	}
	
	//prepare output file
	ofstream fout;
	fout.open(filePath.c_str(), (dataIndex<0 ? ios::out : ios::app));
	if(!fout.good()){ error = PlotErrorFileWriteError; return Plot(); }
	if(dataIndex>=0) fout << "\n\n" << ST_GNUPLOT_COMMENT_PREFIX << " Data index: " << dataIndex << "\n";
	error = PlotErrorNone;	
	string nnFileComment = fileComment; 
	STPlot_aux_replaceAll(nnFileComment, "\n", (string(ST_GNUPLOT_COMMENT_PREFIX) + "\n").c_str());

	
	fout 	<< ST_GNUPLOT_COMMENT_PREFIX << " " << nnFileComment << "\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << " Data columns: 1:x" << (xlabel=="" ? "" : "("+xlabel+")"); for(j=0; j<ycount; ++j){ fout << "  " << (j+2) << ":cumulative_until_y" << (j+1) << (names.size()>j ? "("+names[j]+")" : ""); } fout << "\n";
	for(i=0; i<xcount; ++i){
		fout << data[i*(1+ycount)];
		for(j=0, cumulative=0; j<ycount; ++j){
			cumulative += data[i*(1+ycount) + j+1];
			fout << " " << cumulative;
		}
		fout << "\n";
	}
	fout.close();
	
	//define curve plots
	Plot plot;
	plot.title = title;
	plot.axisLabels = MultiLabel(xlabel, ylabel);
	plot.range = MultiRange(Range(), Range());
	plot.sources.push_back(PlotSource(PlotType2DFilledCurve, (0<names.size() ? names[0] : ""), filePath, (dataIndex<0 ? 0 : dataIndex), 1, 2));
	plot.sources[0].fillReference = 0;
	for(j=1; j<ycount; ++j){
		plot.sources.push_back(PlotSource(PlotType2DFilledCurves, (j<names.size() ? names[j] : ""), filePath, (dataIndex<0 ? 0 : dataIndex), 1, 2+j, 1+j));
	}
	if(relative){
		for(j=0, rescaleAtZero=0; j<ycount; ++j){ 
			plot.sources[j].yTransformationFunction = plot.sources[j].zTransformationFunction = "#/$"+STPlot_makeString(1+ycount);
			rescaleAtZero += data[j+1];
		}
	}else{
		rescaleAtZero = 1;
	}
	
	//set fill colors
	vector<Color> extraColors;
	generateDiscreteColors(max(0l,(long)ycount - (long)colors.size()), true, extraColors, ColorBlindnessNone);
	for(j=0; j<ycount; ++j){
		color = (j<colors.size() ? colors[j] : extraColors[j-colors.size()]);
		color.daltonize(correctForColorBlindness);
		plot.sources[j].setColor(color);
		plot.sources[j].filledCurveStyle = FilledCurveStyleAbove;
	}
	
	//set labels
	const long NL = min(ycount, (unsigned long)names.size());
	const double labelSize = (NL<10 ? -1 : max(1.0,12.0 * (10.0/NL)));
	double below;
	for(j=0, below=0; j<NL; ++j){
		plot.labels.push_back(Label(names[j], true, data[0], (0.5*(below+data[j+1])+0.5*below)/rescaleAtZero, LabelAlignmentLeft));
		plot.labels.back().fontSize = labelSize;
		below += data[j+1];
	}
	
	return plot;
}





Plot TimeSeriesDistribution2Plot(	const vector<double> 			&times,				// (INPUT) common times for all sample time series, strictly increasing
									const vector<double>			&centiles,			// (INPUT) centiles (e.g. 0.5, 0.95) at which to show contours and change coloring
									const vector<vector<double> > 	&centileMinValues, 	// (INPUT) values[t][c] is min value of c-th centile at time times[t] 
									const vector<vector<double> > 	&centileMaxValues, 	// (INPUT) values[t][c] is max value of c-th centile at time times[t] . That is, at time t, time series will be within [centileMinValues[t][c], centileMaxValues[t][c]] with probablity centiles[c]
									const vector<double>			&highlightValues,		// (INPUT) values[t] is value of median at time times[t]. Can also be empty (omit).
									const string					&highlightName,
									const vector<double>			&meanValues,
									double							lineWidth,
									ofstream						&fout,
									const string					&filePath,
									long							&currentDataIndex,
									const string					&title,
									const string					&fileComment,
									const string					&xlabel,
									const string					&ylabel,
									const string					&dataName,
									double							epsilon,
									ColorPalette					colorPalette,
									bool 							useKeysInsteadOfColorBar, // if true, centile color scale is shown as plot keys instead of a color bar
									PlotError						&error){
	const long NC = centiles.size();
	const long NT = times.size();
	if(NT<1){ error = PlotErrorNothingToBeDone; return Plot(); }
	if((NC==0) && highlightValues.empty()){ error = PlotErrorNothingToBeDone; return Plot(); }
	
	const bool foutWasOpen = fout.is_open();
	if(!foutWasOpen) fout.open(filePath.c_str());
	if(!fout.is_open()){ error = PlotErrorFileWriteError; return Plot(); }
	fout << (currentDataIndex==0 ? "" : "\n\n") << ST_GNUPLOT_COMMENT_PREFIX << " Data index: " << currentDataIndex << "\n";
	error = PlotErrorNone;
	string nnFileComment = STPlot_aux_get_replaceAll(fileComment, "\n", (string(ST_GNUPLOT_COMMENT_PREFIX) + "\n").c_str());

	// sort centiles just in case
	vector<long> sortedCentileIndices;
	STPlot_qsortIndices(centiles, sortedCentileIndices, true);

	//write file header
	fout 	<< ST_GNUPLOT_COMMENT_PREFIX << " " << nnFileComment << "\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << " Distribution of " << dataName << " time series (resolved at " << NC << " centiles)\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << " 1:(" << xlabel << ")";
	
	// write raw centile data
	long columns=1;
	for(long c=0; c<NC; ++c){
		fout << "\t" << (++columns) << ":" << centiles[sortedCentileIndices[c]] << "_centile_min\t" << (++columns) << ":" << centiles[sortedCentileIndices[c]] << "_centile_max";
	}
	if(!highlightValues.empty()){ fout << "\t" << (++columns) << ":" << highlightName; }
	if(!meanValues.empty()){ fout << "\t" << (++columns) << ":mean"; }
	fout << "\n";
	for(long t=0; t<NT; ++t){
		fout << times[t];
		for(long c=0; c<NC; ++c){
			fout << "\t" << centileMinValues[t][sortedCentileIndices[c]] << "\t" << centileMaxValues[t][sortedCentileIndices[c]];
		}
		if(!highlightValues.empty()){ fout << "\t" << highlightValues[t]; }
		if(!meanValues.empty()){ fout << "\t" << meanValues[t]; }
		fout << "\n";
	}
	
	// get data extremes
	double minValue = (NC==0 ? highlightValues[0] : centileMinValues[0][0]);
	double maxValue = (NC==0 ? highlightValues[0] : centileMaxValues[0][0]);
	for(long t=0; t<NT; ++t){
		for(long c=0; c<NC; ++c){
			minValue = min(minValue, centileMinValues[t][c]);
			maxValue = max(maxValue, centileMaxValues[t][c]);
		}
		if(!highlightValues.empty()){
			minValue = min(minValue, highlightValues[t]);
			maxValue = max(maxValue, highlightValues[t]);
		}
		if(!meanValues.empty()){
			minValue = min(minValue, meanValues[t]);
			maxValue = max(maxValue, meanValues[t]);
		}
	}

	// create plots
	Plot plot;
	plot.title = title;
	plot.range = MultiRange(Range(times.front(), times.back()), Range::getNicePlotRange(minValue, maxValue, false, false, epsilon));
	plot.axisLabels = MultiLabel(xlabel, ylabel);
	for(long c=NC-1; c>=0; --c){
		plot.sources.push_back(PlotSource(PlotType2DFilledCurves,(useKeysInsteadOfColorBar ? STPlot_makeString(100*centiles[sortedCentileIndices[c]])+"%" : ""),filePath,currentDataIndex,1,2+2*c,3+2*c));
		plot.sources.back().setColorAsValue(centiles[sortedCentileIndices[c]]);
		plot.sources.back().filledCurveStyle = FilledCurveStyleAll;
	}
	columns = 1+2*NC;
	if(!highlightValues.empty()){
		// show highlighted time series as 'double curve' (i.e. two alternating colors)
		plot.sources.push_back(PlotSource(PlotType2DCurve,"",filePath,currentDataIndex,1,++columns));
		plot.sources.back().setColorLineTypeLineWidth(Color::white(),1,lineWidth);
		plot.sources.push_back(plot.sources.back());
		plot.sources.back().setColorLineTypeLineWidth(Color::black(),2,lineWidth);
		plot.sources.back().title = highlightName;
	}
	if(!meanValues.empty()){
		// show mean time series as 'double curve' (i.e. two alternating colors)
		plot.sources.push_back(PlotSource(PlotType2DCurve,"",filePath,currentDataIndex,1,++columns));
		plot.sources.back().setColorLineTypeLineWidth(Color::white(),1,lineWidth);
		plot.sources.push_back(plot.sources.back());
		plot.sources.back().setColorLineTypeLineWidth(Color::black(),3,lineWidth);
		plot.sources.back().title = "mean";
	}
	++currentDataIndex;
	
	// add color bar or keys
	if((!highlightValues.empty()) || (!meanValues.empty()) || useKeysInsteadOfColorBar){
		plot.showKeys = true;
	}
	if(!useKeysInsteadOfColorBar){
		plot.showColorBox = true;
		plot.CBRange = Range(0,centiles[sortedCentileIndices.back()]);
		plot.colorPalette = colorPalette;
		plot.CBTics.disableAutomaticTics();
		plot.CBLabel = "centile";
		for(long c=0; c<NC; ++c){ plot.CBTics.add(STPlot_makeString(100*centiles[c],3)+"%%", centiles[c]); }
	}
	
	return plot;
}




Plot irregular2DGridValueMap2Plot(	const vector<double>	&gridValues, 	// (INPUT) data grid values as matrix in row-major format
									long					NR,				// (INPUT) number of grid rows
									long					NC,				// (INPUT) number of grid columns
									const vector<double>	&rowValues,		// (INPUT) coordinates at row centres, in ascenting order. Of size at least NR
									const vector<double>	&columnValues,	// (INPUT) coordinates at column centres, in ascending order. Of size at least NC
									ofstream				&fout,
									const string			&filePath,
									long					&currentDataIndex,		
									const string			&title,
									const string			&fileComment,
									const string			&xlabel,
									const string			&ylabel,
									const string			&dataName,
									bool					alsoSaveOriginalData,
									PlotError				&error){
	if(NR*NC==0){
		error = PlotErrorNothingToBeDone;
		return Plot();
	}
	
	//prepare output file
	if(!fout.is_open()){ error = PlotErrorFileWriteError; return Plot(); }
	long r,c;
	error = PlotErrorNone;	
	string nnFileComment = fileComment; 
	STPlot_aux_replaceAll(nnFileComment, "\n", (string(ST_GNUPLOT_COMMENT_PREFIX) + "\n").c_str());
	
	// write original data
	if(alsoSaveOriginalData){
		fout 	<< ST_GNUPLOT_COMMENT_PREFIX << " Data index: " << currentDataIndex << "\n"
			 	<< ST_GNUPLOT_COMMENT_PREFIX << " " << nnFileComment << "\n"
				<< ST_GNUPLOT_COMMENT_PREFIX << " Original data\n"
				<< ST_GNUPLOT_COMMENT_PREFIX << " Data columns: 1:(" << xlabel << ")\t2:(" << ylabel << ")\t3:(" << dataName << ")\n";	
		for(r=0; r<NR; ++r){
			for(c=0; c<NC; ++c){
				fout << columnValues[c] << "\t" << rowValues[r] << "\t" << gridValues[r*NC + c] << "\n";	
			}
			fout << "\n";
		}
		fout << "\n\n";
		++currentDataIndex;
	}

	// write data for plotting
	fout 	<< ST_GNUPLOT_COMMENT_PREFIX << " Data index: " << currentDataIndex << "\n"
		 	<< ST_GNUPLOT_COMMENT_PREFIX << " " << nnFileComment << "\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << " Data modified to generate correct plot in gnuplot. Use with caution\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << " Data columns: 1:(" << xlabel << ", modified)\t2:(" << ylabel << ", modified)\t3:(" << dataName << ")\n";	
	for(r=-1; r<NR; ++r){
		for(c=-1; c<NC; ++c){
			fout << (c==-1 ? columnValues[0] : (c==NC-1 ? columnValues[NC-1] : (columnValues[c]+columnValues[c+1])/2)) << "\t" << (r==-1 ? rowValues[0] : (r==NR-1 ? rowValues[NR-1] : (rowValues[r]+rowValues[r+1])/2)) << "\t" << gridValues[max(0l,r)*NC + max(0l,c)] << "\n";	
		}
		fout << "\n";
	}
	fout << "\n\n";
	
	// create plot
	Plot plot(PlotSource(	PlotType2DScannedValueMap, "", filePath, currentDataIndex, 1, 2, 3),
							MultiLabel(xlabel, ylabel),
							MultiRange(Range(columnValues[0], columnValues[NC-1]), Range(rowValues[0], rowValues[NR-1])),
							title);
	plot.gridColoringType = GridColoringTypeUR;
	++currentDataIndex;
	return plot;
}
									
	
	

Plot circularNetwork2Plot(	const double 			adjacencies[],
							long					N,
							const vector<double>	&sizeParam,			// (INPUT) object size >0 (~ node area) parameter (only relative values matter). Should either contain N entries or be empty.
							const vector<double>	&colorParam,		// (INPUT) object color parameter (only relative values matter). Should either contain N entries or be empty.
							const string			&filePath,
							const string			&title,
							const string			&fileComment,
							ColorPalette			nodeColorPalette,	// (INPUT) color palette from which to choose object colors
							const Color				&colorNeg,
							const Color				&colorPos,
							double					thresholdNeg,		// (INPUT) negative correlation threshold (negative adjacencies greater than this are not plotted)
							double					thresholdPos,		// (INPUT) positive correlation threshold (positive adjacencies smaller than this are not plotted)
							bool					cluster,
							PlotError 				&error){
	long i,j,n;
	double temp,r,g,b,x1,y1,x2,y2;
	ofstream fout;
	fout.open(filePath.c_str());
	if(!fout.good()){
		error = PlotErrorFileWriteError;
		return Plot();
	}
	error 					= PlotErrorNone;
	bool varSizes 			= (N<=sizeParam.size());
	bool varColors 			= (N<=colorParam.size());
	string nnFileComment	= fileComment; 
	STPlot_aux_replaceAll(nnFileComment, "\n", (string(ST_GNUPLOT_COMMENT_PREFIX) + "\n").c_str());
	
	//prepare geometry
	const double circleGap	= 0.2*PI;
	const double domega 	= -(2.0*PI-circleGap)/N; //move clockwise
	const double omegao 	= PI/2 - circleGap/2;
	double omega, sizeScale, maxSize, colorScale, minColorParam, maxColorParam;
	const double defaultNodeRadius = sqrt(ST_DEFAULT_NETWORK_NODE_COVERAGE/(PI*N));
	if(varSizes){
		for(i=0, temp=0, maxSize=sizeParam[0]; i<N; ++i){ 
			temp += abs(sizeParam[i]); 
			maxSize = max(maxSize, sizeParam[i]);
		}
		if(temp<EPSILON){
			varSizes 	= false;
		}else{
			sizeScale 	= ST_DEFAULT_NETWORK_NODE_COVERAGE/temp;
			maxSize		= maxSize*sizeScale;
		}
	}
	if(varColors){
		for(i=0, minColorParam=maxColorParam=colorParam[0]; i<N; ++i){
			minColorParam = min(minColorParam, colorParam[i]);
			maxColorParam = max(maxColorParam, colorParam[i]);
		}
		if(maxColorParam < minColorParam + EPSILON){
			varColors = false;
		}else{
			colorScale = 1/(maxColorParam - minColorParam);
		}
	}
	double const R = max(0.1,0.49 - (varSizes ? sqrt(maxSize/PI) : defaultNodeRadius));	//radius of circular network
	
	//figure out node order on circular network 
	long *order = new long[N];
	if(cluster){
		//hierarchical clustering
		double *distances 	= new double[N*N];
		double *distances2 	= new double[N*N];
		long *merge1 		= new long[N-1];
		long *merge2 		= new long[N-1];
		double minAdj, maxAdj;
		
		//get min & max adjacencies in data set
		for(i=0, minAdj=maxAdj=adjacencies[0]; i<N; ++i){ for(j=i+1; j<N; ++j){ minAdj = min(minAdj, adjacencies[i*N+j]); maxAdj = max(maxAdj, adjacencies[i*N+j]); } }
		double adjSpan = maxAdj-minAdj;
		
		if(adjSpan<EPSILON){	
			//all adjacencies are identical, so distances will be as well. So no reason for clustering
			for(i=0; i<N; ++i){ order[i] = i; }
		}else{			
			//heuristically transform adjacencies into "distances"
			for(i=0; i<N; ++i){
				for(j=i+1; j<N; ++j){
					distances2[i*N+j] = distances[i*N+j] = pow(100.0, 2*(maxAdj+EPSILON - adjacencies[i*N+j])/adjSpan);
				}
			}
			//cluster hierarchically according to distances
			HClusterAgglomerative(N, distances, HClusterLinkageUPGMA, merge1, merge2, NULL);
			//order nodes according to tree structure
			orderHClusterByDistance(distances2, N, merge1, merge2, order);
		}
		delete[] distances, distances2, merge1, merge2;
	}else{
		//no clustering, order stays unchanged
		for(i=0; i<N; ++i){ order[i] = i; }
	}
	
	//write line segments
	fout 	<< ST_GNUPLOT_COMMENT_PREFIX << nnFileComment << "\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << "Circular network for " << N << " nodes" << (cluster ? " (cluster)" : "") << ", cutoff thresholds " << thresholdNeg << " (neg) & " << thresholdPos << " (pos)\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << "Data index 0: Edges (adjacencies) as line segments:\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << "Data columns: node1  node2  x  y  dx  dy  r  g  b  corr\n";
	for(i=0; i<N; ++i){
		x1 = 0.5 + R*cos(omegao + domega*i); 
		y1 = 0.5 + R*sin(omegao + domega*i);
		for(j=i+1; j<N; ++j){
			x2 = 0.5 + R*cos(omegao + domega*j); 
			y2 = 0.5 + R*sin(omegao + domega*j);
			n = (order[i]<order[j] ? order[i]*N + order[j] : order[j]*N + order[i]); // adjacencies[] is only well-defined for upper-triangular part
			if((adjacencies[n]>thresholdNeg) && (adjacencies[n]<thresholdPos)) continue; //correlation too weak
			fout << order[i] << " " << order[j] << " " << x1 << " " << y1 << " " << (x2-x1) << " " << (y2-y1) << " ";
			if(adjacencies[n]<0){
				colorNeg.getRGB(r,g,b);
			}else{
				colorPos.getRGB(r,g,b);
			}
			fout << r << " " << g << " " << b << " " << adjacencies[n] << "\n";
			
		}
	}
	
	//write network nodes (circles)
	fout 	<< "\n\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << "Data index 1: Nodes (objects) as circles\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << "Data columns: node x y radius r g b\n";
	for(i=0; i<N; ++i){
		omega = omegao + domega * i;
		n = order[i];
		r = (varSizes ? sqrt(sizeParam[n]*sizeScale/PI) : defaultNodeRadius);
		fout << order[i] << " " << (0.5 + R*cos(omega)) << " " << (0.5 + R*sin(omega)) << " " << r;
		if(varColors){
			value2RGBViaPalette((colorParam[n]-minColorParam)*colorScale, nodeColorPalette, r, g, b);
			fout << " " << r << " " << g << " " << b << "\n";
		}else{
			fout << " 0 0 0\n"; //no node colors specified, make all nodes black
		}
	}
	
	//write adjacency matrix (not needed for the plot)
	fout 	<< "\n\n"
			<< ST_GNUPLOT_COMMENT_PREFIX << "Data index 2: Adjacency matrix (original order)\n";
	for(i=0; i<N; ++i){
		for(j=0; j<N; ++j){
			fout << (j==0 ? "" : " ") << (i==j ? 0 : adjacencies[(i<j ? i*N+j : j*N+i)]);
		}
		fout << "\n";
	}
	fout.close();
	
	//prepare plot
	Plot plot;
	plot.sources.push_back(PlotSource(PlotType2DColoredVectors, "", filePath, 0, 3, 4, 5, 6, 7, 8, 9));
	plot.sources.push_back(PlotSource(PlotType2DColoredVariableSizedSpheres, "", filePath, 1, 2, 3, 4, 5, 6, 7));
	plot.sources[0].vectorHeadStyle = VectorHeadStyleNone;
	plot.sources[0].setLineWidth((N<100 ? 2 : 1));
	plot.axisLabels = MultiLabel("", "");
	plot.title = title;
	plot.range = MultiRange(Range(0,1), Range(0,1));
	plot.showKeys = false;
	plot.xTics = plot.yTics = Tics(0);
	
	delete[] order;
	return plot;

}





#pragma mark -
#pragma mark Hierarchical Clustering
#pragma mark -



void aux_placeHCNode(const long merge1[], const long merge2[], long node, long order[], long &nextSlot){
	if(merge1[node]<=0){ order[nextSlot] = -merge1[node]; ++nextSlot; }
	else{ aux_placeHCNode(merge1, merge2, merge1[node]-1, order, nextSlot); }
	
	if(merge2[node]<=0){ order[nextSlot] = -merge2[node]; ++nextSlot; }
	else{ aux_placeHCNode(merge1, merge2, merge2[node]-1, order, nextSlot); }
}






void HClusterAgglomerative(	long			N,				// (INPUT) sample size
							double 			distances[],	// (INPUT) distance matrix in row-major format (only strictly super-diagonal part needs to be filled). Might be modified upon function return so make a copy of it if you still need it.
							HClusterLinkage	linkage,
							long 			merge1[], 		// (OUTPUT) needs to be pre-allocated of size N-1, or NULL
							long			merge2[],		// (OUTPUT) needs to be pre-allocated of size N-1, or NULL
							long			order[]){		// (OUTPUT) needs to be pre-allocated of size N, or NULL
	long Nclust = N;
	if(N==1){
		if(order!=NULL) order[0] = 0;
		return;
	}
	long n,m,clmin1,clmin2,r,c,c1,c2;
	double td, dmin, d1, d2;
	long mass1, mass2, newmass;
	bool trackMass = (linkage == HClusterLinkageUPGMA);
	long *clusters 	= new long[N];
	long *clID 		= new long[N];
	long *masses; if(trackMass){ masses = new long[N]; }
	long Nmerges 	= 0;
	bool deleteMerge1=false, deleteMerge2=false;
	if(order!=NULL){
		//we need merge1[] and merge2[] to figure out leaf order
		if(merge1==NULL){ merge1 = new long[N-1]; deleteMerge1=true; }
		if(merge2==NULL){ merge2 = new long[N-1]; deleteMerge2=true; }
	}
	//setup cluster data positions (within distances matrix) and IDs (negative for leafs, strictly positive for composite nodes)
	for(n=0; n<N; ++n){ clusters[n] = n; clID[n] = -n; if(trackMass) masses[n] = 1; }
	
	while(Nclust>1){
		//find cluster pair with least distance
		//min-distance pair will be clusters clmin1 & clmin2
		dmin = distances[STPlot_aux_array_rowColumn2Index(clusters[clmin1=0],clusters[clmin2=1],N)]; 
		for(r=0; r<Nclust; ++r){
			for(c=r+1; c<Nclust; ++c){
				td = distances[STPlot_aux_arrayUT_rowColumn2Index(clusters[r],clusters[c],N)];
				if(td<dmin){
					clmin1 = r; clmin2 = c; dmin=td;
				}
			}
		}
		c1 		= clusters[clmin1];
		c2 		= clusters[clmin2];
		if(trackMass){
			mass1 	= masses[c1];
			mass2 	= masses[c2];
			newmass = mass1 + mass2;
		}
		
		//merge this min-distance pair into one cluster (take place of clmin1)
		if(merge1!=NULL){ merge1[Nmerges] = clID[c1]; }
		if(merge2!=NULL){ merge2[Nmerges] = clID[c2]; }
		++Nmerges;
		
		//update cluster distance matrix: Distance to merged clusters depends on distances to the two sub-clusters, as described by the linkage criterion
		for(c=clmin1+1; c<Nclust; ++c){
			d1 	= distances[STPlot_aux_arrayUT_rowColumn2Index(c1, clusters[c], N)];
			d2 	= distances[STPlot_aux_arrayUT_rowColumn2Index(c2, clusters[c], N)];
			n	= STPlot_aux_array_rowColumn2Index(c1, clusters[c], N);
			switch(linkage){
			case HClusterLinkageSingle: 	distances[n]=min(d1, d2); break;
			case HClusterLinkageComplete: 	distances[n]=max(d1, d2); break;
			case HClusterLinkageUPGMA: 		distances[n]=(d1*mass1 + d2*mass2)/newmass; break;
			}
		}
		for(r=0; r<clmin1; ++r){
			d1 	= distances[STPlot_aux_arrayUT_rowColumn2Index(clusters[r], c1, N)];
			d2 	= distances[STPlot_aux_arrayUT_rowColumn2Index(clusters[r], c2, N)];
			n	= STPlot_aux_array_rowColumn2Index(clusters[r], c1, N);
			switch(linkage){
			case HClusterLinkageSingle: 	distances[n]=min(d1, d2); break;
			case HClusterLinkageComplete: 	distances[n]=max(d1, d2); break;
			case HClusterLinkageUPGMA: 		distances[n]=(d1*mass1 + d2*mass2)/newmass; break;
			}
		}
		
		//merge clmin1 & clmin2 and delete old cluster clmin2	
		if(trackMass){ masses[c1] = newmass; }	
		clID[c1] = Nmerges; 
		clusters[clmin2] = clusters[Nclust-1];
		--Nclust;
	}
		
	if(order != NULL){
		//order leafs so that tree branches do not cross
		//simple ordering algorithm: recursively traverse binary tree in left-child first order
		long nextSlot=0;
		aux_placeHCNode(merge1, merge2, Nmerges-1, order, nextSlot);
	}
	
//	delete[] clID, clusters, distances;
	if(deleteMerge1) delete[] merge1;
	if(deleteMerge2) delete[] merge2;
	if(trackMass) delete[] masses;
}



//auxiliary function for setting subtree ranges and reordering subtrees
void aux_HCsplitRange(	long p, 
						long rangeMin, 
						long rangeMax, 
						long ndown, 
						long nup, 
						long leafs[],
						long PRangeL[], 
						long PRangeU[], 
						long merge1[], 
						long merge2[],
						const long Pmass[]){
	if(merge1!=NULL) merge1[p] = ndown; 
	if(merge2!=NULL) merge2[p] = nup;
	if(ndown<=0){ leafs[-ndown] = rangeMin; }
	else{ PRangeL[ndown-1] = rangeMin; PRangeU[ndown-1] = rangeMin+Pmass[ndown-1]-1; } 
	
	if(nup<=0){ leafs[-nup] = rangeMax; }
	else{ PRangeL[nup-1] = rangeMax-Pmass[nup-1]+1; PRangeU[nup-1] = rangeMax; }
}





void orderHClusterRowsCentroid(	const double	data[],		//(INPUT) data matrix in row-major format
								long			Nrows,
								long			Ncolumns,
								bool			transpose,
								HClusterMetric	metric,
								long 			merge1[], 	//(INPUT) (might be modified)
								long 			merge2[], 	//(INPUT) (might be modified)
								long			order[]){ 	//(OUTPUT) improved row order (if transpose==false) or column order (if transpose==true)
	long n,m,p,c1,c2,gchild1,gchild2,sibling,child,ci,si;
	double d1,d2;
	long mass1,mass2;
	bool keepOrder;
	long N 			= (transpose ? Ncolumns : Nrows);
	long dim		= (transpose ? Nrows : Ncolumns);
	long Nmerges 	= N-1;

	// calculate mass & centroids for internal nodes
	double 	*Pcentroid 	= new double[Nmerges*dim]; 	// centroids of parent nodes (to be calculated). Matrix in row-major format. Each row of this matrix is one centroid of a different node.
	long 	*Pmass 		= new long[Nmerges];		// mass of parent nodes (to be calculated)
	for(p=0; p<Nmerges; ++p){
		c1 			= merge1[p];
		c2 			= merge2[p];
		mass1 		= (c1<=0 ? 1 : Pmass[c1-1]);
		mass2 		= (c2<=0 ? 1 : Pmass[c2-1]);
		Pmass[p]	= mass1 + mass2;
		for(n=0; n<dim; ++n){
			Pcentroid[p*dim+n] = ((c1>0 ? Pcentroid[(c1-1)*dim+n] : data[STPlot_aux_array_rowColumn2Index(Ncolumns, transpose, -c1, n)])*mass1 + (c2>0 ? Pcentroid[(c2-1)*dim+n] : data[STPlot_aux_array_rowColumn2Index(Ncolumns, transpose, -c2, n)])*mass2)/Pmass[p];
		}
	}	
	
	// order nodes according to the proximity of their centroids to the centroid of their parent’s sibling
	long *PRangeL 	= new long[Nmerges];
	long *PRangeU 	= new long[Nmerges];	// PRangeL[i] - PRangeU[i] is the range within which the leafs descending from node i are to be placed
	long *leafs 	= new long[N];		// leafs[i] is the new position of the i-th object
	
	
	
	// root leaf range is entire index space (1,..,N)
	// root children order does not matter
	PRangeL[Nmerges-1] = 0; PRangeU[Nmerges-1] = N-1; 
	aux_HCsplitRange(	Nmerges-1, 
						0, 
						N-1, 
						merge1[Nmerges-1], 
						merge2[Nmerges-1],
						leafs,
						PRangeL,
						PRangeU,
						merge1,
						merge2,
						Pmass);
												
	for(p=Nmerges-1; p>=0; --p){
		for(ci=0; ci<=1; ++ci){ // reorder children of each child ci (i.e. all grandchildren of p)
			child = (ci==0 ? merge1[p] : merge2[p]);
			if(child<=0) continue; // child is a leaf, so has no grandchildren
			child 	-= 1;
			gchild1 = merge1[child];
			gchild2 = merge2[child];
			
			// get distance of grandchildren to sibling's centroid
			si = 1-ci; sibling = (si==0 ? merge1[p] : merge2[p]); 
			d1 = STPlot_aux_distanceBetweenDataRecords((sibling<=0 ? data 		: Pcentroid),
												(sibling<=0 ? Ncolumns 	: dim),
												(sibling<=0 ? -sibling 	: sibling-1),
												(sibling<=0 ? !transpose: true),
												(gchild1<=0 ? data 		: Pcentroid),
												(gchild1<=0 ? Ncolumns 	: dim),
												(gchild1<=0 ? -gchild1 	: gchild1-1),
												(gchild1<=0 ? !transpose: true),
												dim,
												metric);
			d2 = STPlot_aux_distanceBetweenDataRecords((sibling<=0 ? data 		: Pcentroid),
												(sibling<=0 ? Ncolumns 	: dim),
												(sibling<=0 ? -sibling 	: sibling-1),
												(sibling<=0 ? !transpose: true),
												(gchild2<=0 ? data 		: Pcentroid),
												(gchild2<=0 ? Ncolumns 	: dim),
												(gchild2<=0 ? -gchild2 	: gchild2-1),
												(gchild2<=0 ? !transpose: true),
												dim,
												metric);
							
			// determine appropriate order of grandchildren
			keepOrder = ((d1<d2)&(si<ci))|((d1>d2)&(si>ci));
			
			// set appropriate order and leaf ranges for grandchildren
			aux_HCsplitRange(	child, 
								PRangeL[child], 
								PRangeU[child], 
								(keepOrder ? gchild1 : gchild2), 
								(keepOrder ? gchild2 : gchild1),
								leafs,
								PRangeL, 
								PRangeU, 
								merge1, 
								merge2,
								Pmass);

		}		
	}

	//order[] is inverse permutation of leafs[]
	for(p=0; p<N; ++p){ order[leafs[p]] = p; }	
	delete[] Pcentroid, Pmass, PRangeL, PRangeU, leafs;	
}






double aux_distanceBetweenHClusters(	long N,
										const double leafDistances[], //distance matrix in row-major format. Only strictly upper-triangular part needed
										const long merge1[],
										const long merge2[],
										const long Pmass[],
										long cl1,		// <=0 if leaf, >0 if composite cluster
										long cl2,		// <=0 if leaf, >0 if composite cluster
										double clDistances[]){ //Distance matrix between all clusters, in row-major format. Only strictly upper-triangular part needed. Already weighted by Pmass[]. Negative if distance has not been calculated yet for a given cluster pair.
	long i,j,n,m;
	if(cl2>0){
		//cl2 is a composite cluster
		long child1 = merge1[cl2-1], child2 = merge2[cl2-1];
		i = N-1+cl1; 
		j = N-1+child1; 
		n = (i<j ? i*N+j : j*N+i);
		if(clDistances[n]<0){ 
			clDistances[n] = (child1<=0 ? 1 : Pmass[child1-1])*aux_distanceBetweenHClusters(N, leafDistances, merge1, merge2, Pmass, cl1, child1, clDistances)/Pmass[cl2-1];
		}
		j = N-1+child2; 
		m = (i<j ? i*N+j : j*N+i);
		if(clDistances[m]<0){ 
			clDistances[m] = (child2<=0 ? 1 : Pmass[child2-1])*aux_distanceBetweenHClusters(N, leafDistances, merge1, merge2, Pmass, cl1, child2, clDistances)/Pmass[cl2-1];
		}
		return (clDistances[n] + clDistances[m]);
	}
	
	if(cl1>0){
		//cl1 is a composite cluster but cl2 is a leaf
		long child1 = merge1[cl1-1], child2 = merge2[cl1-1];
		i = N-1+cl2; 
		j = N-1+child1; 
		n = (i<j ? i*N+j : j*N+i);
		if(clDistances[n]<0){ 
			clDistances[n] = (child1<=0 ? 1 : Pmass[child1-1])*aux_distanceBetweenHClusters(N, leafDistances, merge1, merge2, Pmass, cl2, child1, clDistances)/Pmass[cl1-1];
		}
		j = N-1+child2; 
		m = (i<j ? i*N+j : j*N+i);
		if(clDistances[m]<0){ 
			clDistances[m] = (child2<=0 ? 1 : Pmass[child2-1])*aux_distanceBetweenHClusters(N, leafDistances, merge1, merge2, Pmass, cl2, child2, clDistances)/Pmass[cl1-1];
		}
		return clDistances[n] + clDistances[m];
	}
	
	//both cl1 & cl2 are leafs
	cl1 = -cl1; cl2 = -cl2;
	return leafDistances[(cl1<cl2 ? cl1*N+cl2 : cl2*N+cl1)];
}





void orderHClusterByDistance(	const double	distances[], 
								long			N, 
								long			merge1[], 
								long			merge2[], 
								long			order[]){
	long n,m,p,c1,c2,gchild1,gchild2,sibling,child,ci,si;
	double d1,d2;
	long mass1,mass2;
	bool keepOrder;
	long Nmerges = N-1;

	// order nodes according to their average distance to their parent’s sibling
	long *PRangeL 	= new long[Nmerges];
	long *PRangeU 	= new long[Nmerges];	// PRangeL[i] - PRangeU[i] is the range within which the leafs descending from node i are to be placed
	long *leafs 	= new long[N];			// leafs[i] is the new position of the i-th object
	
	//calculate cluster mass
	long 	*Pmass 	= new long[Nmerges];
	for(p=0; p<Nmerges; ++p){
		c1 			= merge1[p];
		c2 			= merge2[p];
		mass1 		= (c1<=0 ? 1 : Pmass[c1-1]);
		mass2 		= (c2<=0 ? 1 : Pmass[c2-1]);
		Pmass[p]	= mass1 + mass2;
	}
	
	//recursively calculate distances between clusters only on demand
	double *clDistances = new double[(N+Nmerges)*(N+Nmerges)];
	for(n=0; n<N; ++n){
		for(m=0; m<N; ++m){
			c1 = N-1-n; c2 = N-1-m;
			clDistances[(c1<c2 ? c1*N+c2 : c2*N+c1)] = distances[(n<m ? n*N+m : m*N+n)];
		}
	}
	
	// root leaf range is entire index space (1,..,N)
	// root children order does not matter
	PRangeL[Nmerges-1] = 0; PRangeU[Nmerges-1] = N-1; 
	aux_HCsplitRange(	Nmerges-1, 
						0, 
						N-1, 
						merge1[Nmerges-1], 
						merge2[Nmerges-1],
						leafs,
						PRangeL,
						PRangeU,
						merge1,
						merge2,
						Pmass);
												
	for(p=Nmerges-1; p>=0; --p){
		for(ci=0; ci<=1; ++ci){ // reorder children of each child ci (i.e. all grandchildren of p)
			child = (ci==0 ? merge1[p] : merge2[p]);
			if(child<=0) continue; // child is a leaf, so has no grandchildren
			child 	-= 1;
			gchild1 = merge1[child];
			gchild2 = merge2[child];
			
			// get distance of grandchildren to sibling's centroid
			si = 1-ci; sibling = (si==0 ? merge1[p] : merge2[p]); 
			d1 = aux_distanceBetweenHClusters(	N,
												distances,
												merge1,
												merge2,
												Pmass,
												sibling,
												gchild1,
												clDistances);
			d2 = aux_distanceBetweenHClusters(	N,
												distances,
												merge1,
												merge2,
												Pmass,
												sibling,
												gchild2,
												clDistances);
							
			// determine appropriate order of grandchildren
			keepOrder = ((d1<d2)&(si<ci))|((d1>d2)&(si>ci));
			
			// set appropriate order and leaf ranges for grandchildren
			aux_HCsplitRange(	child, 
								PRangeL[child], 
								PRangeU[child], 
								(keepOrder ? gchild1 : gchild2), 
								(keepOrder ? gchild2 : gchild1),
								leafs,
								PRangeL, 
								PRangeU, 
								merge1, 
								merge2,
								Pmass);

		}		
	}

	//order[] is inverse permutation of leafs[]
	for(p=0; p<N; ++p){ order[leafs[p]] = p; }	
	delete[] Pmass, PRangeL, PRangeU, leafs, clDistances;	
}
								
								


void orderHClusterAscending(long		N,
							const long	merge1[],
							const long	merge2[],
							long		order[]){
	long p,c1,c2;
	double ID1, ID2;
	long mass1, mass2;
	long Nmerges = N-1;
	// calculate average original orders (ID) for internal nodes
	float 	*PID 	= new float[Nmerges]; 	// average IDs of parent nodes (to be calculated). In sync with merge tables.
	long 	*Pmass 	= new long[Nmerges];	// mass of parent nodes (to be calculated)
	for(p=0; p<Nmerges; ++p){
		c1 			= merge1[p];
		c2 			= merge2[p];
		mass1 		= (c1<=0 ? 1 : Pmass[c1-1]);
		mass2 		= (c2<=0 ? 1 : Pmass[c2-1]);
		Pmass[p]	= mass1 + mass2;
		PID[p] 		= (((float)mass1)*(c1<=0 ? (-c1) : PID[c1-1]) + ((float)mass2)*(c2<=0 ? (-c2) : PID[c2-1]))/((float)Pmass[p]);
	}
	
	// order nodes according to the proximity of their centroids to the centroid of their parent’s sibling
	long *PRangeL 	= new long[Nmerges];
	long *PRangeU 	= new long[Nmerges];	// PRangeL[i] - PRangeU[i] is the range within which the leafs descending from node i are to be placed
	long *leafs 	= new long[N];			// leafs[i] is the new position of the i-th M-row	
		
	// root range is entire index space (1,..,N)
	PRangeL[Nmerges-1] = 0; PRangeU[Nmerges-1] = N-1; 

	for(p=Nmerges-1; p>=0; --p){
		c1 	= merge1[p];
		c2 	= merge2[p];
		ID1	= (c1<=0 ? (-c1) : PID[c1-1]);
		ID2	= (c2<=0 ? (-c2) : PID[c2-1]);
		
		// set appropriate order and leaf ranges for children
		aux_HCsplitRange(	p, 
							PRangeL[p], 
							PRangeU[p], 
							(ID1<ID2 ? c1 : c2), 
							(ID1<ID2 ? c2 : c1),
							leafs,
							PRangeL, 
							PRangeU,
							NULL,
							NULL,
							Pmass);
	}

	//order[] is inverse permutation of leafs[]
	for(p=0; p<N; ++p){ order[leafs[p]] = p; }
	delete[] PID, Pmass, PRangeL, PRangeU, leafs;
}







#pragma mark - 
#pragma mark Plot handlers
#pragma mark


bool plotWithLog(	MultiPlot 		&MP,
					const string 	&logFile,		//set to empty for no log (script) file
					const string 	&plotPDFFile,	//set to empty for no PDF plotting
					const string 	&plotSVGFile,	//set to empty for no SVG plotting
					const string 	&plotName,		//descriptive name for this plot (only relevant if streamOut != NULL)
					ostream* 		streamOut, 		//feedback stream. Set to NULL for no feedback 
					bool 			plotTerm,
					const string	&shellScript,
					string 			&errorMessage){
	return plotWithLog(	MP,
						logFile,	
						plotPDFFile,
						plotSVGFile,
						plotName,	
						streamOut, 	
						plotTerm,	
						shellScript,
						ST_DEFAULT_TERMINAL_WIDTH,
						ST_DEFAULT_TERMINAL_HEIGHT,
						ST_DEFAULT_PDF_WIDTH,
						ST_DEFAULT_PDF_HEIGHT,
						ST_DEFAULT_SVG_WIDTH,
						ST_DEFAULT_SVG_HEIGHT,
						ST_DEFAULT_PLOT_H_SPACING,
						ST_DEFAULT_PLOT_V_SPACING,
						errorMessage);
}




bool plotWithLog(	MultiPlot 		&MP,
					const string 	&logFile,	
					const string 	&plotPDFFile,
					const string 	&plotSVGFile,
					const string 	&plotName,	
					ostream* 		streamOut, 	
					bool 			plotTerm,	
					const string	&shellScript,
					double			plotWidthTerm,
					double			plotHeightTerm,
					double			plotWidthPDF,
					double			plotHeightPDF,
					double			plotWidthSVG,
					double			plotHeightSVG,
					double			plotHSpacing,
					double			plotVSpacing,
					string 			&errorMessage){
	string plotOutput, plotFilePath;
	ofstream log;
	PlotError error;
	const bool plotPDF = (plotPDFFile != "");
	const bool plotSVG = (plotSVGFile != "");
	if((!plotPDF) && (!plotTerm) && (!plotSVG)) return true; //nothing to do
	if(logFile != ""){ log.open(logFile.c_str()); }
	bool termOK=false, PDFOK=false, logOK=false, SVGOK=false;
	MP.setAspectRatio(1, ST_DEFAULT_MIN_MULTIPLOT_RATIO, ST_DEFAULT_MAX_MULTIPLOT_RATIO);
	MP.setPlotSpacing(plotHSpacing, plotVSpacing);
	errorMessage = "";
	
	if(plotTerm){
		//Plot everything to terminal
		if(log.is_open()){ MP.setLogStream(&log); }
		MP.setOutputToTerminal(ST_DEFAULT_TERM, plotWidthTerm, plotHeightTerm, plotName);
		if((error = MP.plot()) != PlotErrorNone){
			errorMessage = Plot::plotErrorDescription(error);
			if(streamOut != NULL) (*streamOut) << "ERROR: " << errorMessage << endl;
		}else{
			termOK = true;
		}
		
	}
	if(plotPDF){
		//Plot everything to PDF
		if((!plotTerm) && log.is_open()){
			MP.setLogStream(&log);
		}else{
			MP.setLogStream(NULL);
		}
		MP.setOutputToPDF(plotPDFFile, plotWidthPDF, plotHeightPDF);
		if((error = MP.plot()) != PlotErrorNone){
			errorMessage = Plot::plotErrorDescription(error);
			if(streamOut != NULL) (*streamOut) << "ERROR: " << errorMessage << endl;
		}else{
			PDFOK = true;
		}
	}
	if(plotSVG){
		//Plot everything to SVG
		if((!plotTerm) && (!plotPDF) && log.is_open()){
			MP.setLogStream(&log);
		}else{
			MP.setLogStream(NULL);
		}
		MP.setOutputToSVG(plotSVGFile, plotWidthSVG, plotHeightSVG);
		if((error = MP.plot()) != PlotErrorNone){
			errorMessage = Plot::plotErrorDescription(error);
			if(streamOut != NULL) (*streamOut) << "ERROR: " << errorMessage << endl;
		}else{
			SVGOK = true;
		}
	}
	if(log.is_open()){
		logOK = !log.fail();
		log.close();
		MP.setLogStream(NULL);
	}
	
	//add gnuplot command to shell script
	if((shellScript!="") && (logFile!="")){
		ofstream shellFout;
		bool newSS = STPlot_aux_file_missing_or_empty(shellScript);
		shellFout.open(shellScript.c_str(), ios_base::app);
		if(shellFout.is_open()){
			if(newSS){ 
				shellFout << "#!/bin/bash\n"; 
				if(system(NULL)){ system(("chmod u+x '"+shellScript+"'").c_str()); }
//				chmod(shellScript.c_str(), S_IRWXU);
			}
			set<string> sources = MP.sourceFiles();
			shellFout << "\n# Plot" << (plotName=="" ? "" : " " + plotName) << ":\n";
			if(plotTerm){ shellFout << "# Will plot to the " << ST_DEFAULT_TERM << " terminal\n"; }
			else if(plotPDF){ shellFout << "# Will generate the PDF plot: " << plotPDFFile << "\n"; }
			else if(plotSVG){ shellFout << "# Will generate the SVG plot: " << plotSVGFile << "\n"; }
			if(!sources.empty()){ shellFout << "# Required source files:\n"; }
			for(set<string>::const_iterator i=sources.begin(); i!=sources.end(); ++i){
				shellFout << "#   " << *i << "\n";
			}
			if(logOK){
				shellFout << "# Required gnuplot script file: " << logFile << "\n";
				shellFout << (MP.foundGnuplot() ? MP.gnuplotPath() : "gnuplot") << " '" << logFile << "'\n";
			}else{
				shellFout << "# Error: Failed to write to file '" << logFile << "'\n";
			}
		}
	}

	//some feedback
	if(streamOut != NULL){
		if(PDFOK) 	plotOutput = "PDF";
		if(SVGOK) 	plotOutput = plotOutput + (plotOutput=="" ? "" : " & ") + "SVG";
		if(termOK) 	plotOutput = plotOutput + (plotOutput=="" ? "" : " & ") + "terminal";
		(*streamOut) << "Successfully plotted " << plotName << " to " << plotOutput << endl;
	}
	
	return ((PDFOK || (!plotPDF)) && (SVGOK || (!plotSVG)) && (termOK || (!plotTerm)));
}










#pragma mark -
#pragma mark Auxiliary functions
#pragma mark -





void STPlot_transformDataToUniformDistribution(	const vector<double> 	&original,
												long 					start,
												long					end,
												vector<double> 			&transformed,
												Tics 					&adjustedTics){
	const long N = end-start+1;
	transformed.resize(N);
	vector<long> indices;
	adjustedTics.disableAutomaticTics();
	if(N<=0) return;
		
	// map data to ascending order values (equivalent to transforming using empirical CDF)
	
	indices.resize(N); for(long i=0; i<N; ++i){ indices[i] = start+i; }
	STPlot_aux_qsortIndices(original, indices, start, end, true);
	long Ndifferent=0;
	for(long i=0; i<N; ++i){
		if(i==0 || original[indices[i]]!=original[indices[i-1]]) ++Ndifferent;
	}
	long order=-1;
	vector<double> orderedDifferentValues(Ndifferent);
	for(long i=0; i<N; ++i){
		if(i==0 || original[indices[i]]!=original[indices[i-1]]) ++order;
		transformed[indices[i]] 		= order/(1.0*max(1l,Ndifferent-1));
		orderedDifferentValues[order] 	= original[indices[i]];
	}

	// create tics placed at regular transformed values but labeled with original values
	vector<long> gridI;
	STPlot_getRegularIntegerGrid(0,Ndifferent-1,gridI);
	vector<double> ticLabels(gridI.size()); for(long i=0; i<gridI.size(); ++i){ ticLabels[i] = orderedDifferentValues[gridI[i]]; }
	const long labelPrecision = STPlot_aux_getPrecisionToDistinguishValues(ticLabels);
	for(long i=0; i<gridI.size(); ++i){
		adjustedTics.add(STPlot_makeString(ticLabels[i], labelPrecision), gridI[i]/(1.0*max(1l,Ndifferent-1)));
	}
}






void STPlot_getRegularIntegerGrid(long minI, long maxI, vector<long> &grid){
	grid.clear();
	if(maxI<minI) return;
	grid.push_back(minI);
	if(minI==maxI) return;
	
	// figure out a preferred grid step
	long step=-1;
	const long Ncount = 13;
	const long preferredN[Ncount] = {5, 4, 6, 7, 3, 8, 9, 10, 11, 13, 17, 19, 21};
	for(int i=0; i<Ncount; ++i){
		if((maxI-minI) % (preferredN[i]-1) == 0){ step = (maxI-minI)/(preferredN[i]-1); break; }
	}
	if(step<=0){
		for(int i=0; i<Ncount; ++i){
			if((maxI-minI-1) % (preferredN[i]-1) == 0){ step = (maxI-minI-1)/(preferredN[i]-1); break; }
		}
	}
	if(step<=0) step = maxI-minI;
	
	// add points at preferred regular steps
	grid.reserve(long((maxI-minI)/step + 2));
	for(int j=1; minI+j*step<=maxI; ++j){
		grid.push_back(minI + j*step);
	}
}


double STPlot_aux_distanceBetweenDataRows(	const double 	data[],	// (INPUT) data matrix in row-major format
											long			row1,
											long			row2,
											long			Ncolumns,
											HClusterMetric	metric){
	double D = 0;
	if(metric == HClusterMetricEuclidean){
		for(long c=0; c<Ncolumns; ++c){
			D += STPlot_aux_SQ(data[row1*Ncolumns + c] - data[row2*Ncolumns + c]);
		}
		D = sqrt(D);
	}else if(metric == HClusterMetricManhattan){
		for(long c=0; c<Ncolumns; ++c){
			D += abs(data[row1*Ncolumns + c] - data[row2*Ncolumns + c]);
		}
	}else{ //HClusterMetricMaximum
		for(long c=0; c<Ncolumns; ++c){
			D = max(D, abs(data[row1*Ncolumns + c] - data[row2*Ncolumns + c]));
		}
	}
	return D;
}




double STPlot_aux_distanceBetweenDataColumns(	const double 	data[],	// (INPUT) data matrix in row-major format
												long			column1,
												long			column2,
												long			Nrows,
												long			Ncolumns,
												HClusterMetric	metric){
	double D = 0;
	if(metric == HClusterMetricEuclidean){
		for(long r=0; r<Nrows; ++r){
			D += STPlot_aux_SQ(data[r*Ncolumns + column1] - data[r*Ncolumns + column2]);
		}
		D = sqrt(D);
	}else if(metric == HClusterMetricManhattan){
		for(long r=0; r<Nrows; ++r){
			D += abs(data[r*Ncolumns + column1] - data[r*Ncolumns + column2]);
		}
	}else{ //HClusterMetricMaximum
		for(long r=0; r<Nrows; ++r){
			D = max(D, abs(data[r*Ncolumns + column1] - data[r*Ncolumns + column2]));
		}
	}
		return D;
}



double STPlot_aux_distanceBetweenDataRecords(	const double 	data1[],		// (INPUT) data matrix 1 in row-major format
												long			Ncolumns1,
												long			recordID1,
												bool			asRow1,
												const double	data2[],		// (INPUT) data matrix 2 in row-major format
												long			Ncolumns2,
												long			recordID2,
												bool			asRow2,
												long			recordDim,		// (INPUT)
												HClusterMetric	metric){
	double td,D = 0;
	for(long n=0; n<recordDim; ++n){
		td = data1[asRow1 ? (recordID1*Ncolumns1 + n) : (n*Ncolumns1+recordID1)] - data2[asRow2 ? (recordID2*Ncolumns2 + n) : (n*Ncolumns2 + recordID2)];
		switch(metric){
		case HClusterMetricEuclidean: D += STPlot_aux_SQ(td); break;
		case HClusterMetricManhattan: D += abs(td); break;
		default: D = max(D, abs(td)); break;
		}
	}
	if(metric == HClusterMetricEuclidean){ D = sqrt(D); }
	return D;
}




bool STPlot_aux_file_missing_or_empty(string file){
	ifstream fin;
	fin.open(file.c_str());
	return fin.fail();
}


int STPlot_aux_replaceAll(string& haystack, const char *needle, const char *replacement){
	int pos = 0, count = 0; 
	int needleLen = strlen(needle), replacementLen = strlen(replacement);
	while((pos = haystack.find(needle, pos))!= string::npos){
		haystack.replace(pos, needleLen, replacement);
		pos += replacementLen;
		++count;
	}
	return count;
}




string STPlot_aux_get_replaceAll(const string& haystack, const char *needle, const char *replacement){
	string s=haystack;
	STPlot_aux_replaceAll(s, needle, replacement);
	return s;
}


double STPlot_aux_modulo(double numerator, double denominator){
	if((numerator>=0) && (denominator >= 0)){
		return numerator - denominator * floor(numerator/denominator);
	}else if((numerator < 0) && (denominator < 0)){
		return numerator - denominator * floor(numerator/denominator);
	}else if((numerator < 0) && (denominator > 0)){
		return numerator + denominator * ceil(abs(numerator/denominator));
	}else if((numerator >0) && (denominator < 0)){
		return numerator + denominator * ceil(numerator/abs(denominator));
	}
	return 0; //should not have arrived here
}




bool STPlot_aux_fileExists(const string &path){ 
	struct stat fileInfo; 
	if(stat(path.c_str(),&fileInfo) == 0){ 
		return true; 
	}else{ 
		return false; 
	}
}





void STPlot_qsortIndices(const vector<double> &values, vector<long> &sortedIndices, bool ascending){
	sortedIndices.resize(values.size());
	for(long n=0; n<sortedIndices.size(); ++n) sortedIndices[n] = n;
	STPlot_aux_qsortIndices(values, sortedIndices, 0, sortedIndices.size()-1, ascending);
}


long STPlot_aux_qsortIndices_partition(const vector<double> &values, vector<long> &indices, long start, long end, bool ascending){
	const long pi = start+(end-start)/2; //middle
	const double &pv = values[indices[pi]];
	long temp;
	
	// swap
	temp = indices[pi]; indices[pi] = indices[end]; indices[end] = temp;
	
	long si = start;
	for(long i=start; i<end; ++i){
		if(ascending ? values[indices[i]] < pv : values[indices[i]] > pv){
			// swap
			temp = indices[i]; indices[i] = indices[si]; indices[si] = temp;
			++si;
		}
	}
	// swap
	temp = indices[si]; indices[si] = indices[end]; indices[end] = temp;
	
	return si;
}


void STPlot_aux_qsortIndices(const vector<double> &values, vector<long> &indices, long start, long end, bool ascending){
	if(start<end){
		long p = STPlot_aux_qsortIndices_partition(values, indices, start, end, ascending);
		STPlot_aux_qsortIndices(values, indices, start, p-1, ascending);
		STPlot_aux_qsortIndices(values, indices, p+1, end, ascending);
	}
}


void STPlot_aux_correctFileExtensionFromTo(string &path, const string &from, const string &to){
	const long dotpos = path.rfind('.');
	if(dotpos==string::npos) return;
	const long slashpos = path.rfind('/');
	if((slashpos!=string::npos) && (dotpos<slashpos)) return;
	if(STPlot_aux_lowercase(path.substr(dotpos+1))==STPlot_aux_lowercase(from)) path.replace(dotpos+1,path.length()-dotpos-1,to);
}




string STPlot_aux_getUniqueTempFilePath(){
	const string basePath = string(P_tmpdir)+"/STPLOT";
	string path = basePath;
	long counter=0;
	while(STPlot_aux_fileExists(path)){
		path = basePath + STPlot_makeString(counter);
		++counter;
	}
	return path;
}


long STPlot_aux_getPrecisionToDistinguishValues(const vector<double> &values){
	double minDifference = -1;
	for(long i=1; i<values.size(); ++i){
		if((values[i]!=values[i-1]) && (values[i]!=0 || values[i-1]!=0)){
			minDifference = min(minDifference, abs(values[i]-values[i-1])/max(abs(values[i]), abs(values[i-1])));
		}
	}
	if(minDifference<=0) return 3;
	return max(0.0,ceil(-log10(minDifference)))+2;
}







}//END OF NAMESPACE ST

#endif