function [fittedp pcaanalysis calcconcs calcabs]=newmain(filename,ncomp,method,cutoff,totconc,baselinecorrect,maxcut) 
%Needs: filename(string), ncomp, method(string), cutoff, totconc, baselinecorrect(bool)
% calls: getdata, trimdata, getpca, efa, opa, rtl, model, modres, oparef, mcr, normalise, altlsq, plotresults

%need octave-odepkg, octave-optim

%edited 26/11/2012
% filename='myfile.dat'
% ncomp=
% method='efa','pka','rtl','opa','ore'
% cutoff=210
% totconc=1
% baselinecorrect=1,0
% maxcut=500
%
%
%
tic
% parallel execution
	noprocs=5; %define number of cores/instances to use
%
%	tim1=time; another thing matlab can't do
%%% UGLY
disp('execution begun')
toc
tic
if (or(method=='pka',method=='rtl'))
		do_mcrals=false;
	else
		do_mcrals=true;
end
%do_mcrals=false; % set to false if pka or rtl
%disp('init mcrals done')
%toc
%tic

% if you use 'rtl' you must supply these variables
	xr=	[1,0,0,0;
		0,1,0,0;
		0,0,1,0;];

	xp=	[0,1,0,0; % same number of rows and columns as xr
		0,0,1,0;
		0,0,0,1];

	c0=[1,0,0,0].*abs(totconc); % same number of columns as in xr
	pguess=[0.07;0.08;0.09]; % number of rows must be equal to number of rows in xr
	fittedp=0;
% end of rtl parameters
disp('init rtl done')
%toc
%tic

% if you use 'pka' you must supply these variables
	kguess=[7,9]
	kguess=10.^(-kguess);
	fittedk=0;
% end of pka parameters
disp('init pka done')
%toc
tic

%%% Begin program %%%

% read data
	[spectra wavevector timevector]=getdata(filename);
disp('read data done')
toc
tic

% apply baseline
	if (baselinecorrect)
		wavelength='last';
		spectra=baseline(spectra,wavelength,wavevector);
	end
disp('baseline correct done')
toc
tic

% trim data
	[obsdata wavevector]=trimdata(spectra,cutoff,wavevector,maxcut);
disp('data trim done')
toc
tic

tmp_obsdata=size(obsdata); %another hack forced by matlab

% simple pca
	if and((tmp_obsdata(2)<20000),(tmp_obsdata(1)<20000))
		[pcaanalysis]=getpca(obsdata);
	else
		pcaanalysis=[1:10' 1:10' 1:10' 1:10'];
	end
disp('basic pca done')
toc
tic

% choose and apply method

	if (method=='efa') % evolving factor analysis
		[calcforward calcbackward calcconcs calcabs]=efa(obsdata,ncomp,noprocs);
		firstdimension='col';
		resids=obsdata-calcconcs*calcabs;

	elseif (method=='pka') % fit pka
		phvector=10.^(-timevector);
		[calcabs calcconcs fittedk]=pka(obsdata, phvector, kguess);
		firstdimension='col';
		resids=obsdata-calcconcs*calcabs;

	elseif (method=='opa') % orthogonal projection approach
		[calcabs calcconcs]=opa(obsdata, ncomp);
		firstdimension='row';
		resids=obsdata-calcconcs*calcabs;

	elseif (method=='rtl') % apply kinetic model
		[calcabs calcconcs fittedp]=rtl(xr,xp,obsdata,c0,timevector,pguess);
		firstdimension='col';
		resids=obsdata-calcconcs*calcabs;

%	elseif (method=='ore') % referenced opa
%		xref=load('-ascii','ref.out');
%		if (columns(spectra)~=columns(xref))
%			disp('There are different numbers of wavelengths in your experimental data')
%			disp('and your reference spectra! Your results will be nonsense')
%		end
%		resids=obsdata-calcconcs*calcabs;
%		xref=xref(rows(xref),columns(xref)-columns(wavevector)+2:columns(xref));
%		[calcabs calcconcs]=oparef(obsdata,ncomp,xref');
%		firstdimension='row';

	end

	disp('Sizes of concentration and absorbance matrices')
	[size(calcconcs) size(calcabs)]
disp('method done')
toc
%tic

% do mcr
% is it wise to do mcr for rtl?
%	if (and((method~='rtl'),(method~='pka')))
%		[calcabs calcconcs rms finalrms]=mcr(obsdata,calcconcs,calcabs,firstdimension);
%		resids=obsdata-calcconcs*calcabs;
%	end		
%disp('mcr done')
%toc
tic

%normalise
	if (and((method~='rtl'),(method~='pka')))
		if (totconc~=-1)
			[calcconcs calcabs]=normalise(calcconcs,obsdata,totconc);
		end
	end
disp('normalise done')
toc
tic

% do alternating least squares
	if ( or( (and((method~='rtl'),(method~='pka'))), (do_mcrals==true) ) )
		[calcabs calcconcs resids]=altlsq(calcconcs,calcabs,obsdata,ncomp,method);
	end
disp('altlsq done')
toc
tic

% plot results
	plotresults(obsdata, timevector,wavevector, calcconcs',calcabs, pcaanalysis,resids)
disp('plotted results')
toc
tic

% save results

	concs=[timevector calcconcs]; %t vs conc
	cmpds=[wavevector' calcabs']; % wavelen vs abs for cmpds/ncomp
	calcspecs=[wavevector' (calcconcs*calcabs)'];
	obsspecs=[wavevector' obsdata'];
	resids=[wavevector' resids'];
%	zparams=[filename; ncomp;method;cutoff;totconc;baselinecorrect];
	disp('prepared results')
toc
tic
	save ('-ascii','zparams.out','filename','ncomp','method','cutoff','totconc','baselinecorrect')
	save 'zconcs.out' concs
	save 'zcmpds.out' cmpds 
	save 'zcalcspecs.out' calcspecs
	save 'zresids.out' resids
	save 'zobsspec.out' obsspecs
	save 'zprcomp.out' pcaanalysis
	
	disp('And now we are done done.')
%	exectime=time-tim1; another thing not supported by matlab
	disp('sum of residuals')
	totalsumofresiduals=sum(sum(abs(resids)))
	if ((fittedp==0)&(method=='pka'))
		fittedp=fittedk;
	end
disp('saved results')
toc

end
