diff --git a/src/archive/newsrc.tar b/src/archive/newsrc.tar
deleted file mode 100644
index cbc313c..0000000
Binary files a/src/archive/newsrc.tar and /dev/null differ
diff --git a/src/archive/oldsrc.tar b/src/archive/oldsrc.tar
deleted file mode 100644
index 8f1eca2..0000000
Binary files a/src/archive/oldsrc.tar and /dev/null differ
diff --git a/src/inc/def.hpp b/src/inc/def.hpp
deleted file mode 100644
index 84fa165..0000000
--- a/src/inc/def.hpp
+++ /dev/null
@@ -1,36 +0,0 @@
-/*
- * def.hpp
- *
- *  Created on: May 24, 2013
- *      Author: abduz
- */
-
-#ifndef DEF_HPP_
-#define DEF_HPP_
-
-#include <alglib/ap.h>
-#include <iostream>
-#include <iomanip>
-
-using namespace alglib;
-using namespace std;
-
-typedef real_1d_array vec;
-typedef real_2d_array vec2;
-typedef integer_1d_array ivec;
-
-
-class lcurve {
-public:
-	int len; double dt; vec lc,lce,t;
-	lcurve(int n=1){t.setlength(n); lc.setlength(n); lce.setlength(n) ;len=n;dt=1.0;}
-	lcurve(vec& t_,vec& lc_,vec& lce_,double dt_){t=t_;lc=lc_;lce=lce_;len=t.length();dt=dt_;}
-	double demean(){double m=0;int i;for(i=0;i<len;i++){m+=lc[i];} m/=len;for(i=0;i<len;i++){lc[i]-=m;}; return m;}
-	friend ostream& operator<<(ostream&o,lcurve&l){ o << setprecision(12);
-			for(int i=0;i<l.len;i++){o << l.t[i] << " " << l.lc[i] << " " << l.lce[i] << endl;}
-			return(o);}
-};
-
-
-
-#endif /* DEF_HPP_ */
diff --git a/src/inc/lag.hpp b/src/inc/lag.hpp
deleted file mode 100644
index e82d4a6..0000000
--- a/src/inc/lag.hpp
+++ /dev/null
@@ -1,56 +0,0 @@
-/*
- * lag.hpp
- *
- *  Created on: Jun 1, 2013
- *      Author: azoghbi
- */
-
-#ifndef LAG_HPP_
-#define LAG_HPP_
-
-#include "mod.hpp"
-
-class lag : public mod {
-
-	int		n1;
-	vec		p1,p2;
-	ivec	icx,iphi;
-
-	void	_C( vec );
-	void	_dC( vec , int );
-	void	_dCx( vec , int );
-	void	_dPhi( vec , int );
-
-
-public:
-	lag( lcurve , lcurve , vec , vec );
-	virtual ~lag();
-
-	void	step_pars( int , vec& , vec& );
-	void	print_pars( vec& , vec& );
-};
-
-// ----------------------------------- //
-
-class lag10 : public mod {
-
-	int		n1;
-	vec		p1,p2;
-	ivec	icx,iphi;
-
-	void	_C( vec );
-	void	_dC( vec , int );
-	void	_dCx( vec , int );
-	void	_dPhi( vec , int );
-
-
-public:
-	lag10( lcurve , lcurve , vec , vec );
-	virtual ~lag10();
-
-	void	step_pars( int , vec& , vec& );
-	void	print_pars( vec& , vec& );
-	void	what_pars( int& , int& );
-};
-
-#endif /* LAG_HPP_ */
diff --git a/src/inc/lag_rms.hpp b/src/inc/lag_rms.hpp
deleted file mode 100644
index 5c933e1..0000000
--- a/src/inc/lag_rms.hpp
+++ /dev/null
@@ -1,58 +0,0 @@
-/*
- * lag.hpp
- *
- *  Created on: Jun 1, 2013
- *      Author: azoghbi
- */
-
-#ifndef LAG_RMS_HPP_
-#define LAG_RMS_HPP_
-
-#include "mod.hpp"
-
-class lag_rms : public mod {
-
-	int		n1;
-	vec		p1,p2;
-	ivec	icx,iphi;
-	double	mean1,mean2,mean1sq,mean2sq,mean12;
-
-	void	_C( vec );
-	void	_dC( vec , int );
-	void	_dCx( vec , int );
-	void	_dPhi( vec , int );
-
-
-public:
-	lag_rms( lcurve , lcurve , vec , vec );
-	virtual ~lag_rms();
-
-	void	step_pars( int , vec& , vec& );
-	void	print_pars( vec& , vec& );
-};
-
-// ----------------------------------- //
-
-class lag10_rms : public mod {
-
-	int		n1;
-	vec		p1,p2;
-	ivec	icx,iphi;
-	double	mean1,mean2,mean1sq,mean2sq,mean12;
-
-	void	_C( vec );
-	void	_dC( vec , int );
-	void	_dCx( vec , int );
-	void	_dPhi( vec , int );
-
-
-public:
-	lag10_rms( lcurve , lcurve , vec , vec );
-	virtual ~lag10_rms();
-
-	void	step_pars( int , vec& , vec& );
-	void	print_pars( vec& , vec& );
-	void	what_pars( int& , int& );
-};
-
-#endif /* LAG_HPP_ */
diff --git a/src/inc/main.hpp b/src/inc/main.hpp
deleted file mode 100644
index 72571d0..0000000
--- a/src/inc/main.hpp
+++ /dev/null
@@ -1,58 +0,0 @@
-/*
- * main.hpp
- *
- *  Created on: May 31, 2013
- *      Author: azoghbi
- */
-
-#ifndef MAIN_HPP_
-#define MAIN_HPP_
-
-
-#include <fstream>
-#include <cstring>
-#include <sstream>
-#include <vector>
-
-#include "def.hpp"
-#include "psd.hpp"
-#include "psd_rms.hpp"
-#include "lag.hpp"
-#include "lag_rms.hpp"
-#include "psdlag.hpp"
-#include "psdlag_rms.hpp"
-#include "mcmc.hpp"
-
-
-void usage();
-void do_work( char* );
-void readLC( vector<vector<lcurve> >&, string , int , int , int , bool );
-
-
-double mcmc_lag10( vec& x ,void*ptr ){
-	int i,hnp,np = x.length();hnp=np/2; double phi,twopi=2*M_PI;
-	for(i=0;i<hnp;i++){
-		if(x[i]>7){x[i]=7;}if(x[i]<-7){x[i]=-7;}
-		phi = x[i+hnp]+M_PI;
-		while(phi>twopi){phi -= twopi;}
-		while(phi<0){phi += twopi;}
-		x[i+hnp] = phi-M_PI;
-	}
-	Mod<lag10_rms> *mod = (Mod<lag10_rms>*) ptr; double logl=mod->loglikelihood(x);
-	return logl;
-}
-
-double mcmc_psdlag10( vec& x ,void*ptr ){
-	int i,hnp,np = x.length();hnp=np/2; double phi,twopi=2*M_PI;
-	for(i=0;i<hnp;i++){
-		if(x[i]>7){x[i]=7;}if(x[i]<-7){x[i]=-7;}
-		phi = x[i+hnp]+M_PI;
-		while(phi>twopi){phi -= twopi;}
-		while(phi<0){phi += twopi;}
-		x[i+hnp] = phi-M_PI;
-	}
-	Mod<psdlag10_rms> *mod = (Mod<psdlag10_rms>*) ptr; double logl=mod->loglikelihood(x);
-	return logl;
-}
-
-#endif /* MAIN_HPP_ */
diff --git a/src/inc/mcmc.hpp b/src/inc/mcmc.hpp
deleted file mode 100644
index 7cfb0e4..0000000
--- a/src/inc/mcmc.hpp
+++ /dev/null
@@ -1,35 +0,0 @@
-/*
- * mcmc.hpp
- *
- *  Created on: May 14, 2013
- *      Author: abduz
- */
-
-#ifndef MCMC_HPP_
-#define MCMC_HPP_
-
-#include <alglib/ap.h>
-#include <alglib/alglibmisc.h>
-#include <fstream>
-#include <iostream>
-#include <iomanip>
-using namespace alglib;
-using namespace std;
-
-typedef real_1d_array vec;
-typedef real_2d_array vec2;
-typedef integer_1d_array ivec;
-typedef integer_2d_array ivec2;
-
-class mcmc {
-	int npar; double(*loglikelihood)(vec&,void*); void*ptr;
-	hqrndstate rnd;
-public:
-	int nrun,nburn,nwk,ncheck; double avalue;
-	mcmc(int np,double (*f)(vec&,void*),void *p);
-	void setseed(int s1,int s2 ){hqrndseed(s1,s2,rnd);}
-	void run( vec& , vec& , const char*fname="mcmc.dat");
-	virtual ~mcmc();
-};
-
-#endif /* MCMC_HPP_ */
diff --git a/src/inc/mod.hpp b/src/inc/mod.hpp
deleted file mode 100644
index fa24844..0000000
--- a/src/inc/mod.hpp
+++ /dev/null
@@ -1,289 +0,0 @@
-/*
- * mod.hpp
- *
- *  Created on: May 31, 2013
- *      Author: azoghbi
- */
-
-#ifndef MOD_HPP_
-#define MOD_HPP_
-
-#include "def.hpp"
-#include <alglib/specialfunctions.h>
-#include <alglib/solvers.h>
-
-class mod {
-
-	void	_CSfq( vec , vector<vec2>& , vector<vec2>& );
-
-public:
-
-	int				n,nfq,npar;
-	double			dt,f1,f2;
-	vec				lc,lce,t,FqL;
-	vector<vec2>	Cfq,Sfq,Cv,Cfq2,Sfq2;
-
-	vec2			C,Ci,I,yyT,yyTmC,dC,C2,C3;
-	vec				Cilc;
-
-	ae_int_t		info; densesolverreport rep;
-
-	mod();
-	virtual ~mod();
-
-	virtual void _C( vec );
-	virtual void _dC( vec , int );
-	virtual void step_pars( int , vec& , vec& );
-	virtual void print_pars( vec& , vec& );
-	virtual void print_pars( vec& , vec& , vec& );
-	virtual void what_pars( int&, int& );
-
-	void	setlc(){ lc.setlength(n); lce.setlength(n); t.setlength(n);}
-	void	init( vec , int );
-	double	loglikelihood( vec );
-	void	dlikelihood( vec , double& , vec& , vec2& );
-	void	optimize( vec& , vec& );
-};
-
-// ----------------------------------------- //
-// ---------- mod container ---------------- //
-// ----------------------------------------- //
-
-template <class M>
-class Mod {
-
-	int			nmod,nfq;
-	vector<M>	mods;
-
-public:
-	int			npar;
-	Mod( vector<lcurve> inlc , vec fqL ){
-		int		i;
-		nmod	=	inlc.size();
-		for( i=0 ; i<nmod ; i++ ){ mods.push_back( M( inlc[i] , fqL ) ); }
-		npar	=	mods[0].npar;
-		nfq		=	mods[0].nfq;
-	}
-	Mod( vector<lcurve> lc1, vector<lcurve> lc2 , vec fqL , vec pars ){
-		int		i;
-		nmod	=	lc1.size();
-		for( i=0 ; i<nmod ; i++ ){ mods.push_back( M( lc1[i] , lc2[i] , fqL , pars ) ); }
-		npar	=	mods[0].npar;
-		nfq		=	mods[0].nfq;
-	}
-	Mod( vector<lcurve> lc1, vector<lcurve> lc2 , vec fqL ){
-		int		i;
-		nmod	=	lc1.size();
-		for( i=0 ; i<nmod ; i++ ){ mods.push_back( M( lc1[i] , lc2[i] , fqL ) ); }
-		npar	=	mods[0].npar;
-		nfq		=	mods[0].nfq;
-	}
-	virtual ~Mod(){}
-
-	double	loglikelihood( vec p ){
-		double l=0; for( int i=0 ; i<nmod ; i++ ){ l += mods[i].loglikelihood(p); }
-		return l;
-	}
-
-	void	dlikelihood( vec p , double& loglike , vec& grad , vec2& hess ){
-		int i,j,nm;vec g; vec2 h; double l;g.setlength(npar);h.setlength(npar,npar);
-		loglike = 0; for( i=0 ; i<npar ; i++ ){ grad[i]=0; for(j=0;j<npar;j++){ hess(i,j) = 0;}}
-
-		for( nm=0 ; nm<nmod ; nm++ ){
-			mods[nm].dlikelihood( p , l , g , h );
-			loglike	+=	l;
-			for( i=0 ; i<npar ; i++ ){ grad[i] += g[i]; for(j=0;j<npar;j++){ hess(i,j) += h(i,j);}}
-		}
-	}
-
-	void	step_pars( int n, vec& dpar , vec& pars ){ mods[0].step_pars( n , dpar , pars );}
-	void	print_pars( vec& pars , vec& errs ){ mods[0].print_pars( pars , errs );}
-	void	print_pars( vec& pars , vec& errs1 , vec& errs2 ){ mods[0].print_pars( pars , errs1 , errs2 );}
-	void	what_pars( int& ip1 , int& ip2 ){ mods[0].what_pars( ip1 , ip2 );}
-
-	double	optimize( vec& pars , vec& errs ){
-		int		nmax = 400;
-		double	tol = 1e-3;
-
-		ae_int_t info; densesolverreport rep;
-		int		i,j,n,sing=0,np;
-		double	loglike,dpmax,dl,prevl=-1e20;
-		vec		tmppars,dpar,grad,prevp; vec2 hess,hessi,ii;
-		tmppars.setlength(npar); dpar.setlength(npar); grad.setlength(npar); prevp.setlength(npar);
-		hess.setlength( npar , npar );ii.setlength(npar,npar);hessi.setlength(npar,npar);
-		for( i=0 ; i<npar ; i++ ){
-			tmppars[i] = pars[i]; dpar[i] = 0.01*pars[i]; ii(i,i)=1;for(j=0;j<i;j++){ii(i,j)=0;ii(j,i)=0;}
-			prevp[i] = pars[i];
-		}
-
-		for( n=1 ; n<= nmax ; n++ ){
-
-			dlikelihood( tmppars , loglike , grad , hess );
-			// ---------- //
-			if( loglike != loglike ){ for( i=0 ; i<npar ; i++ ){ dpar[i]/=10;tmppars[i] = prevp[i] + dpar[i]; } continue;}
-			dl = loglike - prevl;
-			// ---------- //
-			vector<int> ic,iv; sing = 0;
-			for( i=0; i <npar ; i++ ){ if( fabs(grad[i]) < 1e-5 ){ sing=1;ic.push_back(i); }else { iv.push_back(i);}}
-			//if( sing == 1 and iv.size()!=0 ){
-			if( sing == 1 ){
-				np = iv.size();
-				vec	g; vec2 h,hi,iii; g.setlength(np); h.setlength(np,np); hi.setlength(np,np);iii.setlength(np,np);
-				for( i=0 ; i<np ; i++ ){
-					iii(i,i) = 1; h(i,i) = hess(iv[i],iv[i]);
-					g[i] = grad[iv[i]]; for(j=0;j<i;j++){ h(i,j) = hess(iv[i],iv[j]);h(j,i)=h(i,j);iii(i,j)=0;iii(j,i)=0;}
-				}
-				spdmatrixcholesky( h , np , false );
-				spdmatrixcholeskysolvem( h , np , false , iii , np , info , rep , hi );
-
-				dpmax = -1e20;
-				for( i=0 ; i<np ; i++ ){
-					dpar[iv[i]] = 0; for( j=0 ; j<np ; j++ ){ dpar[iv[i]] += g[j] * hi(i,j);}
-					dpmax	=	max( dpmax , fabs(dpar[iv[i]]) );
-				}
-				for( i=0 ; i<int(ic.size()); i++ ){ dpar[ic[i]] = 0;}
-			}else {
-			// --------- //
-
-				spdmatrixcholesky( hess , npar , false );
-				spdmatrixcholeskysolvem( hess , npar , false , ii , npar , info , rep , hessi );
-
-				dpmax = -1e20;
-				for( i=0 ; i<npar ; i++ ){
-					dpar[i] = 0; for( j=0 ; j<npar ; j++ ){ dpar[i] += grad[j] * hessi(i,j);}
-					dpmax	=	max( dpmax , fabs(dpar[i]) );
-				}
-			}
-			for(i=0;i<npar;i++){prevp[i] = tmppars[i];}
-			prevl = loglike;
-			if ( dl < 2 and n>10 ){for(i=0;i<npar;i++){dpar[i]/=2;}}
-			step_pars( n , dpar , tmppars );
-
-			cout << n << " " << dpmax << " " << loglike << " "; for(j=0;j<npar;j++){cout << tmppars[j] << " ";} cout << endl;
-			if ( dpmax<tol ) break;
-		}
-		for( i=0 ; i<npar ; i++ ){
-			pars[i] = tmppars[i];
-			errs[i]	= sqrt(hessi(i,i));
-		}
-		cout << setfill('-') << setw(40) << "\n";
-		for(i=0;i<npar;i++){cout << pars[i] << " ";} cout << endl;
-		for(i=0;i<npar;i++){cout << errs[i] << " ";} cout << endl;
-		for(i=0;i<npar;i++){cout << grad[i] << " ";} cout << endl;
-		cout << setfill('-') << setw(40) << "\n" << setfill(' ');
-		print_pars( pars, errs );
-		return loglike;
-	}
-
-	double	opti( vec pars , vec errs ,int k ){
-		int		nmax = 400;
-		double	tol = 1e-3;
-
-		ae_int_t info; densesolverreport rep;
-		int		i,j,n,np;
-		double	loglike,dpmax,dl,prevl=-1e20;
-		vec		tmppars,dpar,grad,prevp; vec2 hess,hessi,ii;
-		tmppars.setlength(npar); dpar.setlength(npar); grad.setlength(npar); prevp.setlength(npar);
-		hess.setlength( npar , npar );ii.setlength(npar,npar);hessi.setlength(npar,npar);
-		vector<int> ic,iv;
-		vec	g; vec2 h,hi,iii;
-
-		dlikelihood( pars , loglike , grad , hess );
-		if( fabs(grad[k])<1e-5 ) { return loglike;}
-
-		for( i=0 ; i<npar ; i++ ){
-			tmppars[i] = pars[i]; dpar[i] = 0.01*pars[i]; ii(i,i)=1;for(j=0;j<i;j++){ii(i,j)=0;ii(j,i)=0;}
-			prevp[i] = pars[i];
-			if(i==k or fabs(grad[i])<1e-5 ){ic.push_back(i);}else{ iv.push_back(i); }
-		}
-		np = npar - ic.size();
-		g.setlength(np); h.setlength(np,np); hi.setlength(np,np);iii.setlength(np,np);
-
-
-		for( n=1 ; n<= nmax ; n++ ){
-
-			dlikelihood( tmppars , loglike , grad , hess );
-			// ---------- //
-			if( loglike != loglike ){ for( i=0 ; i<npar ; i++ ){ tmppars[i] = prevp[i] + dpar[i]/10; } continue;}
-			dl = loglike - prevl;
-			// ---------- //
-
-
-			for( i=0 ; i<np ; i++ ){
-				iii(i,i) = 1; h(i,i) = hess(iv[i],iv[i]);
-				g[i] = grad[iv[i]]; for(j=0;j<i;j++){ h(i,j) = hess(iv[i],iv[j]);h(j,i)=h(i,j);iii(i,j)=0;iii(j,i)=0;}
-			}
-			spdmatrixcholesky( h , np , false );
-			spdmatrixcholeskysolvem( h , np , false , iii , np , info , rep , hi );
-
-			dpmax = -1e20;
-			for( i=0 ; i<np ; i++ ){
-				dpar[iv[i]] = 0; for( j=0 ; j<np ; j++ ){ dpar[iv[i]] += g[j] * hi(i,j);}
-				dpmax	=	max( dpmax , fabs(dpar[iv[i]]) );
-			}
-			for( i=0 ; i<int(ic.size()); i++ ){ dpar[ic[i]] = 0;}
-
-			for(i=0;i<npar;i++){prevp[i] = tmppars[i];}
-			prevl = loglike;
-			if ( dl < 2 and n>10 ){for(i=0;i<npar;i++){dpar[i]/=2;}}
-			step_pars( n , dpar , tmppars );
-
-			cout << n << " " << dpmax << " " << loglike << " "; for(j=0;j<npar;j++){cout << tmppars[j] << " ";} cout << endl;
-			if ( dpmax<tol ) break;
-		}
-		return loglike;
-	}
-
-	void	errors_avg( vec& pars , vec& errs ){
-		vec errs1,errs2;errs1.setlength(npar);errs2.setlength(npar);
-		errors( pars , errs1 , 1 );
-		errors( pars , errs2 , -1 );
-		for(int i=0;i<npar;i++){errs[i] = (errs1[i]+errs2[i])/2;}
-	}
-
-	void	errors( vec& pars , vec& errs1 , vec& errs2 ){
-		errors( pars , errs1 , 1 );
-		errors( pars , errs2 , -1 );
-	}
-
-	void	errors( vec& pars , vec& errs , int sign=1){
-		double		tol = 1e-2 , delchi = 1;
-		int			i,ip,ip1,ip2;
-		double		pu,pd,phalf=0,bestlike,tmplike,dl;
-		vec			pars0,errs0,tmpp; pars0.setlength(npar);errs0.setlength(npar);tmpp.setlength(npar);
-		bestlike	=	optimize( pars , errs );
-		for(i=0;i<npar;i++){pars0[i] = pars[i];errs0[i]=errs[i];}
-
-		what_pars( ip1 , ip2 );
-
-		for( ip=ip1 ; ip<ip2 ; ip++ ){
-			cout << endl << "******* par " << ip << " *******" << endl << endl;
-			if( errs0[ip]>10 ) { continue ;}
-			pd	=	pars[ip];
-			for( i=0 ;i<npar ; i++ ){tmpp[i] = pars0[i]; }
-
-			i=2; dl = 0;
-			while( dl < delchi ){
-				pu			=	pars[ip] + sign*i*errs0[ip];
-				tmpp[ip]	=	pu;
-				tmplike		=	opti( tmpp , errs , ip );
-				dl			=	2*(bestlike - tmplike);
-				i++;
-			}
-			while( fabs(dl-delchi)>tol ){
-				phalf		=	(pu+pd)/2.0;
-				tmpp[ip]	=	phalf;
-				tmplike		=	opti( tmpp , errs , ip );
-				dl			=	2*(bestlike - tmplike);
-				cout << "+++ " << bestlike << " " << tmplike << " " << pars[ip] << " "<< phalf << " " << dl << endl;
-				if( dl<delchi ){ pd = phalf;}else{ pu = phalf;}
-			}
-			errs0[ip]		=	sign* ( phalf - pars0[ip] );
-		}
-		for( i=0 ; i<npar ; i++ ){ pars[i] = pars0[i]; errs[i] = errs0[i];}
-		print_pars( pars0 , errs0 );
-		//optimize( pars0 , errs );
-	}
-};
-
-#endif /* MOD_HPP_ */
diff --git a/src/inc/psd.hpp b/src/inc/psd.hpp
deleted file mode 100644
index 9a3a0f0..0000000
--- a/src/inc/psd.hpp
+++ /dev/null
@@ -1,37 +0,0 @@
-/*
- * psd.hpp
- *
- *  Created on: May 31, 2013
- *      Author: azoghbi
- */
-
-#ifndef PSD_HPP_
-#define PSD_HPP_
-
-#include "mod.hpp"
-
-class psd: public mod {
-
-	void	_C( vec );
-	void	_dC( vec , int );
-
-public:
-	psd( lcurve , vec );
-	virtual ~psd();
-};
-
-// ---------------------- //
-
-class psd10: public mod {
-
-	void	_C( vec );
-	void	_dC( vec , int );
-
-public:
-	psd10( lcurve , vec );
-	virtual ~psd10();
-
-	void	step_pars( int , vec& , vec& );
-};
-
-#endif /* PSD_HPP_ */
diff --git a/src/inc/psd_rms.hpp b/src/inc/psd_rms.hpp
deleted file mode 100644
index dce3e40..0000000
--- a/src/inc/psd_rms.hpp
+++ /dev/null
@@ -1,41 +0,0 @@
-/*
- * psd.hpp
- *
- *  Created on: May 31, 2013
- *      Author: azoghbi
- */
-
-#ifndef PSD_RMS_HPP_
-#define PSD_RMS_HPP_
-
-#include "mod.hpp"
-
-class psd_rms: public mod {
-
-	double	mean,mean2;
-
-	void	_C( vec );
-	void	_dC( vec , int );
-
-public:
-	psd_rms( lcurve , vec );
-	virtual ~psd_rms();
-};
-
-// ---------------------- //
-
-class psd10_rms: public mod {
-
-	double	mean,mean2;
-
-	void	_C( vec );
-	void	_dC( vec , int );
-
-public:
-	psd10_rms( lcurve , vec );
-	virtual ~psd10_rms();
-
-	void	step_pars( int , vec& , vec& );
-};
-
-#endif /* PSD_RMS_HPP_ */
diff --git a/src/inc/psdlag.hpp b/src/inc/psdlag.hpp
deleted file mode 100644
index b486a45..0000000
--- a/src/inc/psdlag.hpp
+++ /dev/null
@@ -1,54 +0,0 @@
-/*
- * psdlag.hpp
- *
- *  Created on: Jun 1, 2013
- *      Author: azoghbi
- */
-
-#include "mod.hpp"
-
-#ifndef PSDLAG_HPP_
-#define PSDLAG_HPP_
-
-class psdlag : public mod {
-
-	int		n1;
-	ivec	ip1,ip2,icx,iphi;
-
-	void	_C( vec );
-	void	_dC( vec , int );
-	void	_dP1( vec , int );
-	void	_dP2( vec , int );
-	void	_dCx( vec , int );
-	void	_dPhi( vec , int );
-
-public:
-	psdlag( lcurve , lcurve , vec );
-	virtual ~psdlag();
-
-	void	step_pars( int , vec& , vec& );
-	void	print_pars( vec& , vec& );
-};
-
-class psdlag10 : public mod {
-
-	int		n1;
-	ivec	ip1,ip2,icx,iphi;
-
-	void	_C( vec );
-	void	_dC( vec , int );
-	void	_dP1( vec , int );
-	void	_dP2( vec , int );
-	void	_dCx( vec , int );
-	void	_dPhi( vec , int );
-
-public:
-	psdlag10( lcurve , lcurve , vec );
-	virtual ~psdlag10();
-
-	void	step_pars( int , vec& , vec& );
-	void	print_pars( vec& , vec& );
-	void	what_pars( int&, int& );
-};
-
-#endif /* PSDLAG_HPP_ */
diff --git a/src/inc/psdlag_rms.hpp b/src/inc/psdlag_rms.hpp
deleted file mode 100644
index 81efdf8..0000000
--- a/src/inc/psdlag_rms.hpp
+++ /dev/null
@@ -1,56 +0,0 @@
-/*
- * psdlag.hpp
- *
- *  Created on: Jun 1, 2013
- *      Author: azoghbi
- */
-
-#include "mod.hpp"
-
-#ifndef PSDLAG_RMS_HPP_
-#define PSDLAG_RMS_HPP_
-
-class psdlag_rms : public mod {
-
-	int		n1;
-	ivec	ip1,ip2,icx,iphi;
-	double	mean1,mean2,mean1sq,mean2sq,mean12;
-
-	void	_C( vec );
-	void	_dC( vec , int );
-	void	_dP1( vec , int );
-	void	_dP2( vec , int );
-	void	_dCx( vec , int );
-	void	_dPhi( vec , int );
-
-public:
-	psdlag_rms( lcurve , lcurve , vec );
-	virtual ~psdlag_rms();
-
-	void	step_pars( int , vec& , vec& );
-	void	print_pars( vec& , vec& );
-};
-
-class psdlag10_rms : public mod {
-
-	int		n1;
-	ivec	ip1,ip2,icx,iphi;
-	double	mean1,mean2,mean1sq,mean2sq,mean12;
-
-	void	_C( vec );
-	void	_dC( vec , int );
-	void	_dP1( vec , int );
-	void	_dP2( vec , int );
-	void	_dCx( vec , int );
-	void	_dPhi( vec , int );
-
-public:
-	psdlag10_rms( lcurve , lcurve , vec );
-	virtual ~psdlag10_rms();
-
-	void	step_pars( int , vec& , vec& );
-	void	print_pars( vec& , vec& );
-	void	what_pars( int&, int& );
-};
-
-#endif /* PSDLAG_RMS_HPP_ */
diff --git a/src/lag.cpp b/src/lag.cpp
deleted file mode 100644
index 3c1cde4..0000000
--- a/src/lag.cpp
+++ /dev/null
@@ -1,227 +0,0 @@
-/*
- * lag.cpp
- *
- *  Created on: Jun 1, 2013
- *      Author: azoghbi
- */
-
-#include "inc/lag.hpp"
-
-lag::lag( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
-
-	// ----------- initial parameters ------------ //
-	n1		=	lc1.len;
-	n		=	n1 + lc2.len;
-	dt		=	lc1.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	lc1.demean(); lc2.demean();
-	int		i;
-	for( i=0 ; i<n1 ; i++ ){
-		lc[i]	=	lc1.lc[i];
-		lce[i]	=	lc1.lce[i]*lc1.lce[i];
-		t[i]	=	lc1.t[i];
-	}
-	for( i=n1 ; i<n ; i++ ){
-		lc[i]	=	lc2.lc[i-n1];
-		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
-		t[i]	=	lc2.t[i-n1];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 2 );
-	// ------------------------------------------ //
-
-	// --------- constants and indices ---------- //
-	p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
-	for( i=0 ; i<nfq ; i++ ){
-		p1[i]	=	pars[i];
-		p2[i]	=	pars[i+nfq];
-		icx[i]	=	i;
-		iphi[i]	=	i+nfq;
-	}
-	// ------------------------------------------ //
-
-}
-
-lag::~lag() {}
-
-void lag::_C( vec p ){
-	int	i,j,k; double d;
-	for( i=0 ; i<n1 ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); } C(i,j) = d;}
-	}
-	for( i=n1 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); } C(i,j) = d;}
-		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
-			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
-	}
-}
-
-void lag::_dC( vec p , int k ){
-	if( k<nfq ){	_dCx( p , k );
-	}else{			_dPhi( p , k-nfq );}
-}
-
-void lag::_dCx( vec p , int k ){
-	int	i,j; double phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void lag::_dPhi( vec p , int k ){
-	int	i,j; double cx=p[k],phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-
-void lag::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
-		pars[i] += dpar[i]/((n<10)?10:1);
-	}
-}
-
-
-void lag::print_pars( vec& pars , vec& errs ){
-	for( int i=0 ; i<nfq ; i++ ){
-		if( pars[i]<0 ){ pars[i]*=-1; pars[i+nfq] += M_PI;}
-		while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
-		while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
-	}
-	mod::print_pars( pars , errs );
-}
-
-
-// ****************************************** //
-
-
-lag10::lag10( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
-
-	// ----------- initial parameters ------------ //
-	n1		=	lc1.len;
-	n		=	n1 + lc2.len;
-	dt		=	lc1.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	lc1.demean(); lc2.demean();
-	int		i;
-	for( i=0 ; i<n1 ; i++ ){
-		lc[i]	=	lc1.lc[i];
-		lce[i]	=	lc1.lce[i]*lc1.lce[i];
-		t[i]	=	lc1.t[i];
-	}
-	for( i=n1 ; i<n ; i++ ){
-		lc[i]	=	lc2.lc[i-n1];
-		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
-		t[i]	=	lc2.t[i-n1];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 2 );
-	// ------------------------------------------ //
-
-	// --------- constants and indices ---------- //
-	p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
-	for( i=0 ; i<nfq ; i++ ){
-		p1[i]	=	pow(10,pars[i]);
-		p2[i]	=	pow(10,pars[i+nfq]);
-		icx[i]	=	i;
-		iphi[i]	=	i+nfq;
-	}
-	// ------------------------------------------ //
-
-}
-
-lag10::~lag10() {}
-
-void lag10::_C( vec p ){
-	int	i,j,k; double d; for(i=0;i<nfq;i++){p[i] = pow(10,p[i]);}
-	for( i=0 ; i<n1 ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); }
-		d += f1*p1[0]*Cfq2[0](i,i) + f2*p1[nfq-1]*Cfq2[2](i,i);
-		d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); }
-		d += f1*p1[0]*Cfq2[0](i,j) + f2*p1[nfq-1]*Cfq2[2](i,j);
-		C(i,j) = d;}
-	}
-	for( i=n1 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); }
-		d += f1*p2[0]*Cfq2[0](i,i) + f2*p2[nfq-1]*Cfq2[2](i,i);
-		d += lce[i]; C(i,i) = d;
-		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); }
-		d += f1*p2[0]*Cfq2[0](i,j) + f2*p2[nfq-1]*Cfq2[2](i,j);
-		C(i,j) = d;}
-		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
-			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); }
-		d += f1*p[icx[0]] * ( Cfq2[0](i,j)*cos(p[iphi[0]]) - Sfq2[0](i,j)*sin(p[iphi[0]]) );
-		d += f2*p[icx[nfq-1]] * ( Cfq2[2](i,j)*cos(p[iphi[nfq-1]]) - Sfq2[2](i,j)*sin(p[iphi[nfq-1]]) );
-		C(i,j) = d;}
-	}
-}
-
-void lag10::_dC( vec p , int k ){
-	if( k<nfq ){	_dCx( p , k );
-	}else{			_dPhi( p , k-nfq );}
-}
-
-void lag10::_dCx( vec p , int k ){
-	int	i,j; double phi=p[iphi[k]],cx = log(10)*pow(10,p[k]);
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void lag10::_dPhi( vec p , int k ){
-	int	i,j; double cx=pow(10,p[k]),phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-
-void lag10::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
-		//pars[i] += dpar[i];
-		pars[i] += dpar[i]/((n<5)?10:1);
-	}
-}
-
-
-void lag10::print_pars( vec& pars , vec& errs ){
-	for( int i=0 ; i<nfq ; i++ ){
-		while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
-		while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
-	}
-	mod::print_pars( pars , errs );
-}
-
-void lag10::what_pars( int& ip1 , int& ip2 ){
-	ip1 = nfq; ip2 = npar;
-}
diff --git a/src/lag_rms.cpp b/src/lag_rms.cpp
deleted file mode 100644
index 9ac355a..0000000
--- a/src/lag_rms.cpp
+++ /dev/null
@@ -1,230 +0,0 @@
-/*
- * lag.cpp
- *
- *  Created on: Jun 1, 2013
- *      Author: azoghbi
- */
-
-#include "inc/lag_rms.hpp"
-
-lag_rms::lag_rms( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
-
-	// ----------- initial parameters ------------ //
-	n1		=	lc1.len;
-	n		=	n1 + lc2.len;
-	dt		=	lc1.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	mean1	=	lc1.demean(); mean2	=	lc2.demean();
-	mean1sq	=	mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
-	int		i;
-	for( i=0 ; i<n1 ; i++ ){
-		lc[i]	=	lc1.lc[i];
-		lce[i]	=	lc1.lce[i]*lc1.lce[i];
-		t[i]	=	lc1.t[i];
-	}
-	for( i=n1 ; i<n ; i++ ){
-		lc[i]	=	lc2.lc[i-n1];
-		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
-		t[i]	=	lc2.t[i-n1];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 2 );
-	// ------------------------------------------ //
-
-	// --------- constants and indices ---------- //
-	p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
-	for( i=0 ; i<nfq ; i++ ){
-		p1[i]	=	pars[i] * mean1sq;
-		p2[i]	=	pars[i+nfq] * mean2sq;
-		icx[i]	=	i;
-		iphi[i]	=	i+nfq;
-	}
-	// ------------------------------------------ //
-
-}
-
-lag_rms::~lag_rms() {}
-
-void lag_rms::_C( vec p ){
-	int	i,j,k; double d;
-	for( i=0 ; i<nfq ; i++ ){ p[icx[i]] *= mean12;}
-	for( i=0 ; i<n1 ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); } C(i,j) = d;}
-	}
-	for( i=n1 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); } C(i,j) = d;}
-		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
-			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
-	}
-}
-
-void lag_rms::_dC( vec p , int k ){
-	if( k<nfq ){	_dCx( p , k );
-	}else{			_dPhi( p , k-nfq );}
-}
-
-void lag_rms::_dCx( vec p , int k ){
-	int	i,j; double phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void lag_rms::_dPhi( vec p , int k ){
-	int	i,j; double cx=p[k]*mean12,phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-
-void lag_rms::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
-		pars[i] += dpar[i]/((n<10)?10:1);
-	}
-}
-
-
-void lag_rms::print_pars( vec& pars , vec& errs ){
-	for( int i=0 ; i<nfq ; i++ ){
-		if( pars[i]<0 ){ pars[i]*=-1; pars[i+nfq] += M_PI;}
-		while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
-		while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
-	}
-	mod::print_pars( pars , errs );
-}
-
-
-// ****************************************** //
-
-
-lag10_rms::lag10_rms( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
-
-	// ----------- initial parameters ------------ //
-	n1		=	lc1.len;
-	n		=	n1 + lc2.len;
-	dt		=	lc1.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	mean1	=	lc1.demean(); mean2	=	lc2.demean();
-	mean1sq	=	mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
-	int		i;
-	for( i=0 ; i<n1 ; i++ ){
-		lc[i]	=	lc1.lc[i];
-		lce[i]	=	lc1.lce[i]*lc1.lce[i];
-		t[i]	=	lc1.t[i];
-	}
-	for( i=n1 ; i<n ; i++ ){
-		lc[i]	=	lc2.lc[i-n1];
-		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
-		t[i]	=	lc2.t[i-n1];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 2 );
-	// ------------------------------------------ //
-
-	// --------- constants and indices ---------- //
-	p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
-	for( i=0 ; i<nfq ; i++ ){
-		p1[i]	=	pow(10,pars[i]) * mean1sq;
-		p2[i]	=	pow(10,pars[i+nfq]) * mean2sq;
-		icx[i]	=	i;
-		iphi[i]	=	i+nfq;
-	}
-	// ------------------------------------------ //
-
-}
-
-lag10_rms::~lag10_rms() {}
-
-void lag10_rms::_C( vec p ){
-	int	i,j,k; double d; for(i=0;i<nfq;i++){p[icx[i]] = pow(10,p[i]) * mean12;}
-	for( i=0 ; i<n1 ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); }
-		d += f1*p1[0]*Cfq2[0](i,i) + f2*p1[nfq-1]*Cfq2[2](i,i);
-		d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); }
-		d += f1*p1[0]*Cfq2[0](i,j) + f2*p1[nfq-1]*Cfq2[2](i,j);
-		C(i,j) = d;}
-	}
-	for( i=n1 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); }
-		d += f1*p2[0]*Cfq2[0](i,i) + f2*p2[nfq-1]*Cfq2[2](i,i);
-		d += lce[i]; C(i,i) = d;
-		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); }
-		d += f1*p2[0]*Cfq2[0](i,j) + f2*p2[nfq-1]*Cfq2[2](i,j);
-		C(i,j) = d;}
-		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
-			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); }
-		d += f1*p[icx[0]] * ( Cfq2[0](i,j)*cos(p[iphi[0]]) - Sfq2[0](i,j)*sin(p[iphi[0]]) );
-		d += f2*p[icx[nfq-1]] * ( Cfq2[2](i,j)*cos(p[iphi[nfq-1]]) - Sfq2[2](i,j)*sin(p[iphi[nfq-1]]) );
-		C(i,j) = d;}
-	}
-}
-
-void lag10_rms::_dC( vec p , int k ){
-	if( k<nfq ){	_dCx( p , k );
-	}else{			_dPhi( p , k-nfq );}
-}
-
-void lag10_rms::_dCx( vec p , int k ){
-	int	i,j; double phi=p[iphi[k]],cx = log(10)*pow(10,p[k])*mean12;
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void lag10_rms::_dPhi( vec p , int k ){
-	int	i,j; double cx=pow(10,p[k])*mean12,phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-
-void lag10_rms::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
-		//pars[i] += dpar[i];
-		pars[i] += dpar[i]/((n<5)?10:1);
-	}
-}
-
-
-void lag10_rms::print_pars( vec& pars , vec& errs ){
-	for( int i=0 ; i<nfq ; i++ ){
-		while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
-		while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
-	}
-	mod::print_pars( pars , errs );
-}
-
-void lag10_rms::what_pars( int& ip1 , int& ip2 ){
-	ip1 = nfq; ip2 = npar;
-}
diff --git a/src/main.cpp b/src/main.cpp
deleted file mode 100644
index a17e05c..0000000
--- a/src/main.cpp
+++ /dev/null
@@ -1,221 +0,0 @@
-/*
- * main.cpp
- *
- *  Created on: May 31, 2013
- *      Author: azoghbi
- */
-
-#include "inc/main.hpp"
-
-/**************************************
- * Main function
- * Checks if there is an input and
- * raises an error.
- * otherwise, call do_work(fname)
- **************************************/
-int main( int argc , char* argv[] ){
-	if ( argc < 2 ){
-		cerr << "** Error ** No input file given." << endl;
-		usage();
-	}else {
-		do_work( argv[1] );
-	}
-	return 0;
-}
-
-
-
-
-/**************************************
- * Main function that does the work
- * INPUT:	char*	fname: input file name
- * RETURN:	void
- **************************************/
-
-void do_work( char* fname ){
-
-	/************   Reading the input file  ***********/
-	ifstream fp(fname);
-	string line,mcmcfile,colon; stringstream ss;
-	int i,nfiles,nfq,mode,npar,bin1,bin2,fit_type,nrun,nburn,nwk,ref; bool strict;
-	getline(fp,line);ss.str(line); ss >> nfiles;								ss.clear();
-	vector<string> files(nfiles);ivec secL;secL.setlength(nfiles);for(i=0;i<nfiles;i++){
-	getline(fp,line); ss.str(line); ss >> files[i] >> secL[i];					ss.clear();}
-	getline(fp,line); ss.str(line); ss >> nfq;vec fqL;fqL.setlength(nfq+1);
-	for(i=0;i<=nfq;i++){ss >> fqL[i];}											ss.clear();
-	getline(fp,line); ss.str(line); ss >> mode;									ss.clear();
-	npar = (mode==0)?4*nfq:nfq; vec pars; pars.setlength(npar);
-	getline(fp,line); ss.str(line); for(i=0;i<npar;i++){ss >> pars[i];}			ss.clear();
-	getline(fp,line); ss.str(line); ss >> ref >> colon >> bin2;			ss.clear();
-	if(colon.compare(0,1,":") == 0 ){
-		bin1	=	atoi(colon.substr(1).c_str());
-	}else{
-		bin1	=	ref; ref = -10; bin2 = atoi(colon.c_str());
-	}
-	if(mode!=0){bin1 = (mode==-1)?-1:(mode-1);bin2=-1;}
-	getline(fp,line); ss.str(line); ss >> fit_type;								ss.clear();
-	getline(fp,line); ss.str(line); ss >> i; strict	= i!=0;						ss.clear();
-	getline(fp,line); ss.str(line); ss >> nrun >> nburn >> nwk >> mcmcfile;		ss.clear();
-	fp.close();
-	/*********   END Reading the input file  **********/
-
-	/**********  Read the light curves  **************/
-	vector<vector<lcurve> > LC,LC_ref;
-	if ( ref == -10 ){
-		for(i=0;i<nfiles;i++)  readLC( LC , files[i] , secL[i] , bin1 , bin2 , strict );
-	}else {
-		for(i=0;i<nfiles;i++) {
-			if( i!=ref ) { readLC( LC , files[i] , secL[i] , bin2 , -1 , strict ); }
-		}
-		readLC( LC_ref , files[ref] , secL[ref] , bin1 , -1 , strict );
-	}
-	/********  END Read the light curves  ************/
-
-	if( mode > 0 or mode==-1 ){
-		vec errs; errs.setlength(nfq);
-		vector<lcurve> lc1; for( i=0 ; i<int(LC.size()) ; i++ ){ lc1.push_back( LC[i][0]); }
-		Mod<psd10_rms>	p1( lc1 , fqL );
-
-		p1.optimize( pars , errs );
-		if( fit_type==1 ) p1.errors( pars , errs );
-	}else if( mode == 0 ){
-		vec ec,e1,e2,errs; ec.setlength(2*nfq);e1.setlength(nfq);e2.setlength(nfq);errs.setlength(4*nfq);
-		vector<lcurve> lc1,lc2;
-		if( ref==-10 ){ // standard way of getting reference
-			for( i=0 ; i<int(LC.size()) ; i++ ){ lc1.push_back( LC[i][0]); lc2.push_back( LC[i][1]); }
-		}else{
-			for( i=0 ; i<int(LC_ref.size()) ; i++ ){ lc2.push_back( LC[i][0]);}
-			for( i=0 ; i<int(LC.size()) ; i++ ){ lc1.push_back( LC_ref[i][0]);}
-		}
-		Mod<psd10_rms>	p1( lc1 , fqL ), p2( lc2 , fqL );
-		vec		pars1,pars2;pars1.setlength(nfq);pars2.setlength(nfq);for(i=0;i<nfq;i++){pars1[i]=pars[i];pars2[i]=pars[i+nfq];}
-
-		p1.optimize( pars1 , e1 );
-		p2.optimize( pars2 , e2 );
-		if( fit_type==1 ) {
-			p1.errors( pars1 , e1 );
-			p2.errors( pars2 , e2 );
-		}
-
-		vec		pc; pc.setlength(2*nfq);
-		for(i=0;i<nfq;i++){pc[i]=pars[i+2*nfq];pc[i+nfq]=pars[i+3*nfq];pars[i]=pars1[i];pars[i+nfq]=pars2[i];pc[i]=(pars1[i]+pars2[i])/2;}
-		Mod<lag10_rms>		l( lc1 , lc2 , fqL , pars );
-		l.optimize( pc , ec );
-		if( fit_type==1 ) l.errors( pc , ec );
-		if( fit_type==2 ) {
-			mcmc mc( 2*nfq , mcmc_lag10 , (void*)&l );
-			mc.nrun=nrun; mc.nburn=nburn; mc.nwk=nwk;
-			mc.run( pc , ec , mcmcfile.c_str() );
-		}
-
-		for(i=0;i<nfq;i++){ pars[i+2*nfq]=pc[i]; pars[i+3*nfq]=pc[i+nfq];}
-		for(i=0;i<nfq;i++){ errs[i] = e1[i]; errs[i+nfq] = e2[i]; errs[i+2*nfq]=ec[i]; errs[i+3*nfq]=ec[i+nfq];}
-
-		Mod<psdlag10_rms>	pl( lc1 , lc2 , fqL );
-		if( fit_type==3 ) pl.errors( pars , errs );
-		if( fit_type==4 ) {
-			mcmc mc( 4*nfq , mcmc_psdlag10 , (void*)&pl );
-			mc.nrun=nrun; mc.nburn=nburn; mc.nwk=nwk;
-			mc.run( pars , errs , mcmcfile.c_str() );
-		}
-
-		pl.print_pars( pars , errs );
-	}
-}
-
-
-
-
-
-
-
-/**************************************
- * Usage function
- * prints a sample input text
- * INPUT:	none
- * RETURN:	void
- **************************************/
-void usage(){
-	cerr << setw(40) << setfill('*') << "\n";
-	cout << setw(30) << left << "1"							<< "# of lc"			<< endl;
-	cout << setw(30) << left << "lc.dat 0"					<< "name1 secL1"		<< endl;
-	cout << setw(30) << left << "3 1e-5 1e-4 5e-4 1e-3"		<< "nfq fqL1 fqL2.." 	<< endl;
-	cout << setw(30) << left << "1"							<< "mode? 0:lag, n:psd" << endl;
-	cout << setw(30) << left << "1 1 1 1 1 1 1 1 1 1 1 1"	<< "init. pars" 		<< endl;
-	cout << setw(30) << left << "0 1"						<< "bin1, bin2,-1: all" << endl;
-	cout << setw(30) << left << "0"							<< "fit_type" 			<< endl;
-	cout << setw(30) << left << "0"							<< "strict" 			<< endl;
-	cout << setw(30) << left << "100 50 50 mcmc.dat"		<< "nrun,nburn,nwk,file"<< endl;
-	cerr << setw(40) << setfill('*') << "\n";
-}
-
-
-
-
-
-/***************************************
- * READ LIGHTCURVES
- * Read a standard light curve file that has a description at the top
- * the first line is something like: # dT nlc b0 b1 b2 ... b_lc (b0 .. are bins boundaries)
- * Then the content with: time lc1 err1 lc2 err2 ... etc
- *
- * INPUT:
- * 		LC		: a vector of vectors of light curves, where the first dimension is the number of
- * 				light curves (i.e. segments), this is incremented in this function when reading a new file.
- * 				The other dimension is either 1 for psd and 2 for lags, which contains either one light curve
- * 				or two for the lag. Passed by reference
- * 		fname	: name of the light curve file
- * 		secL	: the section length to split the light curves to.
- * 		b1		: the bin1 to read for the first light curve. 0-based, -1 for all
- * 		b2		: the bin2 to read, 0-based, -1 for none (.i.e just psd)
- * 		strict	: a boolean on whether to be strict or not in segmenting the light curve
- * 				e.g. for secL=220: 0 gives 100,100 and 1 gives 100,120
- *
- **************************************/
-void readLC( vector<vector<lcurve> >&LC , string fname , int secL , int b1 , int b2 , bool strict ){
-
-	/*		Initialize file stream, define some variables		*/
-	ifstream fp(fname.c_str()); string line,sdum; stringstream ss;
-	int i,j,nlc,n,nsec,ns,sl; double dt;
-
-	/*		read dt from the description line 					*/
-	getline(fp,line); ss.str(line); ss >> sdum >> dt >> nlc; ss.clear();
-
-	/*		get the number of lines 							*/
-	n=0;while(getline(fp,line)){if(line!="")n++;} fp.clear(); fp.seekg(0,ios::beg);getline(fp,line);
-
-	/*		calculate section lengths							*/
-	if(secL==0){secL=n;}
-	nsec = n/secL; vector<int> secl(nsec,secL); if(not strict) secl[nsec-1] += n-nsec*secL;
-
-	/*		Start looping through the number of sections		*/
-	for(ns=0;ns<nsec;ns++){
-		sl	=	secl[ns]; vector<lcurve> Lc;
-		vec t; vector<vec> lc(nlc),lce(nlc); t.setlength( sl ); for(i=0;i<nlc;i++){lc[i].setlength(sl);lce[i].setlength(sl);}
-
-		// Read the first chunk
-		for( i=0 ; i<sl ; i++ ){
-			getline(fp,line);ss.str(line);ss>>t[i];for(j=0;j<nlc;j++){ss>>lc[j][i]>>lce[j][i];} ss.clear();
-		}
-
-		// if we want the whole light curves in bin1
-		if( b1==-1 ){
-			vec lc_,lce_; lc_.setlength( sl );lce_.setlength( sl );for(i=0;i<sl;i++){lc_[i]=0;lce_[i]=0;}
-			for( i=0 ; i<sl ; i++ ){
-				for( j=0 ; j<nlc ; j++ ){ if(j==b2){continue;} lc_[i] += lc[j][i]; lce_[i] += lce[j][i]*lce[j][i]; }
-				lce_[i] = sqrt(lce_[i]);
-			}
-			Lc.push_back(lcurve(t,lc_,lce_,dt));
-
-		// else: i.e. light curve 1 is bin1
-		}else {
-			Lc.push_back(lcurve(t,lc[b1],lce[b1],dt));
-		}
-		// Do we need a second light curve?
-		if(b2!=-1){Lc.push_back(lcurve(t,lc[b2],lce[b2],dt));}
-
-		// push the resulting vector<lcurve> to LC
-		LC.push_back(Lc);
-	}
-	fp.close();
-}
diff --git a/src/makefile b/src/makefile
deleted file mode 100644
index 08bd484..0000000
--- a/src/makefile
+++ /dev/null
@@ -1,13 +0,0 @@
-
-# original
-#libdir='/eos/azoghbi/soft/usr/lib'
-#incdir='/eos/azoghbi/soft/usr/include'
-# most systems
-#libdir='/home/caes/science/psdlag-agn/src/inc'
-#incdir='/home/caes/science/psdlag-agn/src/inc'
-# thor
-libdir='/home/research/xaw5719/psdlag-agn/src/inc'
-alglibdir='/home/research/xaw5719/psdlag-agn/src/inc/alglib'
-incdir='/home/research/xaw5719/psdlag-agn/src/inc'
-psdlag:
-	g++ *cpp inc/alglib/*cpp -o psdlag -O3 -Wall -I${incdir} -L${libdir} -L${alglibdir} -I${alglibdir}
diff --git a/src/mcmc.cpp b/src/mcmc.cpp
deleted file mode 100644
index 565061b..0000000
--- a/src/mcmc.cpp
+++ /dev/null
@@ -1,73 +0,0 @@
-/*
- * mcmc.cpp
- *
- *  Created on: May 14, 2013
- *      Author: abduz
- */
-
-#include "inc/mcmc.hpp"
-
-mcmc::mcmc(int np,double (*f)(vec&,void*),void *p) {
-	npar	=	np;
-	loglikelihood	=	f;
-	ptr		=	p;
-	hqrndrandomize(rnd);
-	nrun = 100; nburn = 100; nwk = 100; ncheck = 10;avalue = 2.0;
-}
-
-mcmc::~mcmc() {}
-
-void mcmc::run( vec& p, vec& pe, const char*fname){
-	int		i,j,k,nacc,nr,n,nsum;
-	vec2	pars; pars.setlength(nwk,npar);
-	vec		prob,currp,bestp,sum,sum2;
-	ivec2	compj;
-	double	bestlike,tmplike=-1e20,z,q,a1=1./sqrt(avalue),a2=sqrt(avalue);;
-	prob.setlength(nwk);currp.setlength(npar);bestp.setlength(npar);sum.setlength(npar);sum2.setlength(npar);
-	compj.setlength(nwk,nwk-1);
-	nsum=0;for(i=0;i<npar;i++){sum[i]=0;sum2[i]=0;}
-	ofstream fout(fname);
-
-	for(k=0;k<nwk;k++){j=0;for(i=0;i<k;i++){compj[k][j]=i;j++;}for(i=k+1;i<nwk;i++){compj[k][j]=i;j++;}}
-
-	bestlike	=	-1e20;
-	for(k=0;k<nwk;k++){
-		//pars[k][0] = p[0];currp[0] = p[0];
-		for(i=0;i<npar;i++){pars[k][i] = p[i]*(1.0+hqrndnormal(rnd)*0.05); currp[i]=pars[k][i];}
-		prob[k]	=	loglikelihood( currp , ptr );
-		if(bestlike<prob[k]){bestlike=prob[k];for(i=0;i<npar;i++){bestp[i]=pars[k][i];}}
-	}
-
-	nr = 0;nacc = 0;
-	for(n=1;n<=(nrun+nburn);n++){
-		for(k=0;k<nwk;k++){
-			nr++;j = compj[k][hqrnduniformi(rnd,nwk-1)];
-			z	=	hqrnduniformr(rnd)*(a2-a1) + a1; z = z*z;
-			for(i=0;i<npar;i++){currp[i]=pars[j][i] + z*(pars[k][i]-pars[j][i]);}
-			tmplike	=	loglikelihood(currp,ptr);
-			q = pow(z,npar-1)*exp(tmplike-prob[k]);
-			if( hqrnduniformr(rnd)<=q and tmplike>-1e10 ){
-				nacc++;for(i=0;i<npar;i++){pars[k][i]=currp[i];} prob[k] = tmplike;
-				if(bestlike<tmplike){bestlike=tmplike;for(i=0;i<npar;i++){bestp[i]=currp[i];}}
-				if(n>nburn) {
-					for(i=0;i<npar;i++){
-						fout << currp[i] << "  ";
-						sum[i] += currp[i];sum2[i] += currp[i]*currp[i]; nsum++;
-					} fout << tmplike << endl;
-				}
-			}
-		}
-		if(n%ncheck==0){
-			q = 1.0*nacc/nr; nr=0;nacc=0;
-			cout << setw(8) << n << setw(12) << q << setw(12) << tmplike << setw(12) << bestlike << endl;
-		}
-	}
-	for(i=0;i<npar;i++){
-		p[i]=bestp[i];
-		pe[i] = sqrt((sum2[i] - sum[i]*sum[i]/nsum)/(nsum+1));
-	}cout << endl;
-
-
-	fout.close();
-
-}
diff --git a/src/mod.cpp b/src/mod.cpp
deleted file mode 100644
index c2dc7b0..0000000
--- a/src/mod.cpp
+++ /dev/null
@@ -1,236 +0,0 @@
-/*
- * mod.cpp
- *
- *  Created on: May 31, 2013
- *      Author: azoghbi
- */
-
-#include "inc/mod.hpp"
-
-mod::mod() {
-	n = 0 ; nfq = 0 ; npar = 0; dt = 0; info = 0;f1=0;f2=0;
-
-}
-
-mod::~mod() {}
-
-void mod::init( vec fqL , int fac ){
-
-	// --------------- Freq stuff----------------- //
-	int		i,j;
-	nfq		=	fqL.length() - 1;
-	npar	=	nfq * fac ;
-	FqL.setlength(nfq+1); for( i=0 ; i<=nfq ; i++ ){ FqL[i] = fqL[i]; }
-	// ------------------------------------------ //
-
-
-	// ---------------- Integrals ---------------- //
-	Cfq.resize( nfq );Sfq.resize( nfq );
-	for( i=0 ; i<nfq ; i++ ){ Cfq[i].setlength( n , n ); Sfq[i].setlength( n , n );}
-	_CSfq( fqL , Cfq , Sfq );
-
-	Cfq2.resize( 3 );Sfq2.resize( 3 );
-	for( i=0 ; i<3 ; i++ ){ Cfq2[i].setlength( n , n ); Sfq2[i].setlength( n , n );}
-	vec		fqL2;fqL2.setlength(4); fqL2[0] = 1e-1*fqL[0]; fqL2[1] = fqL[0]; fqL2[2] = fqL[nfq]; fqL2[3] = 10*fqL[nfq];
-	//f1 = 1.5; f2 = 0.6;
-	f1 = 0.5; f2 = .0;
-	_CSfq( fqL2 , Cfq2 , Sfq2 );
-	// ------------------------------------------ //
-
-
-	// -------------- Containers ---------------- //
-	C.setlength( n , n );
-	Ci.setlength( n , n );
-	I.setlength( n , n );
-	for( i=0 ; i<n ; i++ ){ I(i,i) = 1.0; for( j=0 ; j<i ; j++ ){ I(i,j)=0; I(j,i)=0;}}
-	Cilc.setlength( n );
-	yyT.setlength( n , n );
-	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ yyT(i,j) = lc[i]*lc[j];}}
-	yyTmC.setlength( n , n );
-	dC.setlength( n , n );
-	Cv.resize( npar ); for( i=0 ; i<npar ; i++ ){ Cv[i].setlength( n , n ); }
-	C2.setlength( n , n );
-	C3.setlength( n , n );
-
-	// ------------------------------------------ //
-
-}
-
-void mod::_CSfq( vec fqL , vector<vec2>&cfq , vector<vec2>& sfq ){
-	//fqL[0] *= 1e-6; fqL[nfq]*=2;
-	int		i,j,k, nfq=fqL.length()-1;
-	double	tt,norm,t1,t2,pi=M_PI,pidt=pi*dt,sdt2,sdt1,cdt2,cdt1,s2dt2,s2dt1,c2dt2,c2dt1,dtmt,dtpt;
-	double	st2,st1,ct2,ct1, sm2,sm1,cm2,cm1,  sp2,sp1,cp2,cp1;
-	vec		w; w.setlength(nfq+1); for( k=0 ; k<=nfq ; k++ ){ w[k] = 2*pi*fqL[k]; }
-	norm	=	1./( pidt*pidt );
-
-	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<n ; j++ ){
-		tt	=	t[j] - t[i];
-		dtmt = dt - tt; dtpt = dt+tt;
-		for( k=0 ; k<nfq ; k++ ){
-			if( tt==0 ){
-				sinecosineintegrals( w[k+1]*dt , sdt2 , cdt2 );
-				sinecosineintegrals( w[k]  *dt , sdt1 , cdt1 );
-				t1	=	( cos( w[k+1]*dt ) - 1 ) / ( 2 * fqL[k+1] ) + pidt * sdt2;
-				t1	-=	( cos( w[k]  *dt ) - 1 ) / ( 2 * fqL[k]   ) + pidt * sdt1;
-				cfq[k](i,j)		=	norm * t1;
-				sfq[k](i,j)		=	0;
-			}else if( fabs(tt)==dt ){
-				sinecosineintegrals( w[k+1]*dt , sdt2 , cdt2 );
-				sinecosineintegrals( w[k]  *dt , sdt1 , cdt1 );
-				sinecosineintegrals( 2*w[k+1]*dt , s2dt2 , c2dt2 );
-				sinecosineintegrals( 2*w[k]  *dt , s2dt1 , c2dt1 );
-				t1	=	-( cos( w[k+1]*dt )*pow( sin( 0.5*w[k+1]*dt ) , 2 ) ) / fqL[k+1];
-				t1	-=	-( cos( w[k]  *dt )*pow( sin( 0.5*w[k]  *dt ) , 2 ) ) / fqL[k]  ;
-				t2	=	pidt * ( s2dt2 - sdt2 );
-				t2	-=	pidt * ( s2dt1 - sdt1 );
-				cfq[k](i,j)		=	norm * ( t1 + t2 );
-
-				t1	=	-( 2*cos( 0.5*w[k+1]*dt )*pow( sin( 0.5*w[k+1]*dt ) , 3 ) ) / fqL[k+1];
-				t1	-=	-( 2*cos( 0.5*w[k]  *dt )*pow( sin( 0.5*w[k]  *dt ) , 3 ) ) / fqL[k]  ;
-				t2	=	pidt * ( cdt2 - c2dt2 );
-				t2	-=	pidt * ( cdt1 - c2dt1 );
-				sfq[k](i,j)		=	norm * fabs(tt)/tt * ( t1 + t2 );
-			}else{
-				sinecosineintegrals( w[k+1]*tt , st2 , ct2 );
-				sinecosineintegrals( w[k]  *tt , st1 , ct1 );
-				sinecosineintegrals( w[k+1]*dtmt , sm2 , cm2 );
-				sinecosineintegrals( w[k]  *dtmt , sm1 , cm1 );
-				sinecosineintegrals( w[k+1]*dtpt , sp2 , cp2 );
-				sinecosineintegrals( w[k]  *dtpt , sp1 , cp1 );
-				t1	=	( cos( w[k+1]*dtmt ) - 2*cos( w[k+1]*tt ) + cos( w[k+1]*dtpt ) )/(4*fqL[k+1]);
-				t1	-=	( cos( w[k]  *dtmt ) - 2*cos( w[k]  *tt ) + cos( w[k]  *dtpt ) )/(4*fqL[k]  );
-				t2	=	( dtmt*sm2 - 2*tt*st2 + dtpt*sp2 ) * (pi/2);
-				t2	-=	( dtmt*sm1 - 2*tt*st1 + dtpt*sp1 ) * (pi/2);
-				cfq[k](i,j)		=	norm * ( t1 + t2 );
-
-				t1	=	( -sin( w[k+1]*dtmt ) - 2*sin( w[k+1]*tt ) + sin( w[k+1]*dtpt ) )/(4*fqL[k+1]);
-				t1	-=	( -sin( w[k]  *dtmt ) - 2*sin( w[k]  *tt ) + sin( w[k]  *dtpt ) )/(4*fqL[k]  );
-				t2	=	( dtmt*cm2 + 2*tt*ct2 - dtpt*cp2 ) * (pi/2);
-				t2	-=	( dtmt*cm1 + 2*tt*ct1 - dtpt*cp1 ) * (pi/2);
-				sfq[k](i,j)		=	norm * ( t1 + t2 );
-			}
-		}
-
-	}}
-}
-
-
-void mod::_C( vec ){
-	cerr << endl << "mod::_C is an empty function. It should not be called directly" << endl << endl;;
-	exit(1);
-}
-
-void mod::_dC( vec , int ){
-	cerr << endl << "mod::_C is an empty function. It should not be called directly" << endl << endl;;
-	exit(1);
-}
-
-double mod::loglikelihood( vec p ){
-	double		logdet , chi2 , loglike ; int i;
-	_C( p );
-	if (not spdmatrixcholesky( C , n , false )) {return -1e20;};
-	logdet = 0 ; for( i=0 ; i<n ; i++ ){ logdet += log( C(i,i) ); } logdet *= 2;
-	spdmatrixcholeskysolvem( C , n , false , I , n , info , rep , Ci );
-	rmatrixmv( n , n , Ci , 0 , 0 , 0 , lc , 0 , Cilc , 0 );
-	chi2 = 0 ; for( i=0 ; i<n ; i++ ){ chi2 += Cilc[i]*lc[i] ; }
-	loglike	=	-0.5 * ( chi2 + logdet + n*log(2*M_PI) );
-	return loglike;
-}
-
-void mod::dlikelihood( vec p , double& loglike , vec& grad , vec2& hess ){
-	double		logdet , chi2 , res; int i,j,k;
-	_C( p );
-	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ yyTmC(i,j) = yyT(i,j) - C(i,j); yyTmC(j,i) = yyTmC(i,j);}}
-	spdmatrixcholesky( C , n , false );
-	logdet = 0 ; for( i=0 ; i<n ; i++ ){ logdet += log( C(i,i) ); } logdet *= 2;
-	spdmatrixcholeskysolvem( C , n , false , I , n , info , rep , Ci );
-	rmatrixmv( n , n , Ci , 0 , 0 , 0 , lc , 0 , Cilc , 0 );
-	chi2 = 0 ; for( i=0 ; i<n ; i++ ){ chi2 += Cilc[i]*lc[i] ; }
-	loglike	=	-0.5 * ( chi2 + logdet + n*log(2*M_PI) );
-
-
-	for( i=0 ; i<npar ; i++ ){
-		_dC( p , i );
-		rmatrixgemm( n , n , n , 1.0 , Ci , 0 , 0 , 0 , dC , 0 , 0 , 0 , 0.0 , Cv[i] , 0 , 0 );
-	}
-
-	for( i=0 ; i<npar ; i++ ){
-		rmatrixgemm( n , n , n , 1.0 , Cv[i] , 0 , 0 , 0 , Ci , 0 , 0 , 0 , 0.0 , C2 , 0 , 0 );
-		rmatrixgemm( n , n , n , 1.0 , yyTmC , 0 , 0 , 0 , C2 , 0 , 0 , 0 , 0.0 , C3 , 0 , 0 );
-		res = 0; for( k=0 ; k<n ; k++ ){ res += C3(k,k); }
-		grad[i]		=	0.5 * res;
-		for( j=0 ; j<=i ; j++ ){
-			rmatrixgemm( n , n , n , 1.0 , Cv[i] , 0 , 0 , 0 , Cv[j] , 0 , 0 , 0 , 0.0 , C3 , 0 , 0 );
-			res = 0; for( k=0 ; k<n ; k++ ){ res += C3(k,k); }
-			hess(i,j)	=	0.5 * res; hess(j,i) = 0.5*res;
-		}
-	}
-}
-
-
-void mod::optimize( vec& pars , vec& errs ){
-	int		nmax = 200;
-	double	tol = 1e-3;
-
-	int		i,j,n;
-	double	loglike,dpmax;
-	vec		tmppars,dpar,grad; vec2 hess,hessi,ii;
-	tmppars.setlength(npar); dpar.setlength(npar); grad.setlength(npar);
-	hess.setlength( npar , npar );ii.setlength(npar,npar);hessi.setlength(npar,npar);
-	for( i=0 ; i<npar ; i++ ){
-		tmppars[i] = pars[i]; dpar[i] = 0.01*pars[i]; ii(i,i)=1;for(j=0;j<i;j++){ii(i,j)=0;ii(j,i)=0;}
-	}
-
-	for( n=1 ; n<= nmax ; n++ ){
-
-		dlikelihood( tmppars , loglike , grad , hess );
-		spdmatrixcholesky( hess , npar , false );
-		spdmatrixcholeskysolvem( hess , npar , false , ii , npar , info , rep , hessi );
-
-		dpmax = -1e20;
-		for( i=0 ; i<npar ; i++ ){
-			dpar[i] = 0; for( j=0 ; j<npar ; j++ ){ dpar[i] += grad[j] * hessi(i,j);}
-			dpmax	=	max( dpmax , fabs(dpar[i]) );
-		}
-		step_pars( n , dpar , tmppars );
-
-		cout << n << " " << dpmax << " " << loglike << " "; for(j=0;j<npar;j++){cout << tmppars[j] << " ";} cout << endl;
-		if ( dpmax<tol ) break;
-	}
-	for( i=0 ; i<npar ; i++ ){
-		pars[i] = tmppars[i];
-		errs[i]	= sqrt(hessi(i,i));
-	}
-	cout << setfill('-') << setw(40) << "\n";
-	for(i=0;i<npar;i++){cout << pars[i] << " ";} cout << endl;
-	for(i=0;i<npar;i++){cout << errs[i] << " ";} cout << endl;
-	for(i=0;i<npar;i++){cout << grad[i] << " ";} cout << endl;
-	cout << setfill('-') << setw(40) << "\n" << setfill(' ');
-	print_pars( pars, errs );
-}
-
-
-void mod::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		pars[i] += dpar[i];
-	}
-}
-
-void mod::print_pars( vec& pars, vec& errs ){
-	cout << setfill('*') << setw(40) << "\n" << setfill(' ');
-	cout << endl << "# "; for(int i=0;i<=nfq;i++){cout << FqL[i] << " ";} cout << endl;
-	for(int i=0;i<npar;i++){ printf( "%3.3e %3.3e\n" , pars[i] , errs[i] );}
-	cout << setfill('*') << setw(40) << "\n" << setfill(' ');
-}
-
-void mod::print_pars( vec& pars, vec& errs1 , vec& errs2 ){
-	cout << setfill('*') << setw(40) << "\n" << setfill(' ');
-	cout << endl << "# "; for(int i=0;i<=nfq;i++){cout << FqL[i] << " ";} cout << endl;
-	for(int i=0;i<npar;i++){ printf( "%3.3e %3.3e %3.3e\n" , pars[i] , errs1[i] , errs2[i] );}
-	cout << setfill('*') << setw(40) << "\n" << setfill(' ');
-}
-
-void mod::what_pars( int& ip1 , int& ip2 ){
-	ip1 = 0; ip2 = npar;
-}
diff --git a/src/psd.cpp b/src/psd.cpp
deleted file mode 100644
index 25ec0df..0000000
--- a/src/psd.cpp
+++ /dev/null
@@ -1,105 +0,0 @@
-/*
- * psd.cpp
- *
- *  Created on: May 31, 2013
- *      Author: azoghbi
- */
-
-#include "inc/psd.hpp"
-
-psd::psd( lcurve inlc , vec fqL ) {
-
-	// ----------- initial parameters ------------ //
-	n		=	inlc.len;
-	dt		=	inlc.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	inlc.demean();
-	int		i;
-	for( i=0 ; i<n ; i++ ){
-		lc[i]	=	inlc.lc[i];
-		lce[i]	=	inlc.lce[i]*inlc.lce[i];
-		t[i]	=	inlc.t[i];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 1 );
-	// ------------------------------------------ //
-}
-
-psd::~psd() {}
-
-void psd::_C( vec p ){
-	int	i,j,k; double d;
-	for( i=0 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,i); }
-		d += f1*p[0]*Cfq2[0](i,i) + f2*p[nfq-1]*Cfq2[2](i,i);
-		d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,j); }
-		d += f1*p[0]*Cfq2[0](i,j) + f2*p[nfq-1]*Cfq2[2](i,j);
-		C(i,j) = d;}
-	}
-}
-
-void psd::_dC( vec p , int k ){
-	int	i,j;
-	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);} }
-}
-
-
-// ********************************************** //
-// ********************************************** //
-
-psd10::psd10( lcurve inlc , vec fqL ) {
-
-	// ----------- initial parameters ------------ //
-	n		=	inlc.len;
-	dt		=	inlc.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	inlc.demean();
-	int		i;
-	for( i=0 ; i<n ; i++ ){
-		lc[i]	=	inlc.lc[i];
-		lce[i]	=	inlc.lce[i]*inlc.lce[i];
-		t[i]	=	inlc.t[i];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 1 );
-	// ------------------------------------------ //
-}
-
-psd10::~psd10() {}
-
-void psd10::_C( vec p ){
-	int	i,j,k; double d; for(i=0;i<npar;i++){p[i]=pow(10,p[i]);}
-	for( i=0 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,i); }
-		d += f1*p[0]*Cfq2[0](i,i) + f2*p[nfq-1]*Cfq2[2](i,i);
-		d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,j); }
-		d += f1*p[0]*Cfq2[0](i,j) + f2*p[nfq-1]*Cfq2[2](i,j);
-		C(i,j) = d;}
-	}
-}
-
-void psd10::_dC( vec p , int k ){
-	int	i,j; double d = log(10)*pow(10,p[k]);
-	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = d*Cfq[k](i,j); dC(j,i) = dC(i,j);} }
-}
-
-void psd10::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
-		pars[i] += dpar[i];
-	}
-}
diff --git a/src/psd_rms.cpp b/src/psd_rms.cpp
deleted file mode 100644
index 0c90a2e..0000000
--- a/src/psd_rms.cpp
+++ /dev/null
@@ -1,108 +0,0 @@
-/*
- * psd.cpp
- *
- *  Created on: May 31, 2013
- *      Author: azoghbi
- */
-
-#include "inc/psd_rms.hpp"
-
-psd_rms::psd_rms( lcurve inlc , vec fqL ) {
-
-	// ----------- initial parameters ------------ //
-	n		=	inlc.len;
-	dt		=	inlc.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	mean	=	inlc.demean();
-	mean2	=	mean*mean;
-	int		i;
-	for( i=0 ; i<n ; i++ ){
-		lc[i]	=	inlc.lc[i];
-		lce[i]	=	inlc.lce[i]*inlc.lce[i];
-		t[i]	=	inlc.t[i];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 1 );
-	// ------------------------------------------ //
-}
-
-psd_rms::~psd_rms() {}
-
-void psd_rms::_C( vec p ){
-	int	i,j,k; double d;
-	for( i=0 ; i<npar ; i++ ){ p[i] *= mean2; }
-	for( i=0 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,i); }
-		d += f1*p[0]*Cfq2[0](i,i) + f2*p[nfq-1]*Cfq2[2](i,i);
-		d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,j); }
-		d += f1*p[0]*Cfq2[0](i,j) + f2*p[nfq-1]*Cfq2[2](i,j);
-		C(i,j) = d;}
-	}
-}
-
-void psd_rms::_dC( vec p , int k ){
-	int	i,j;
-	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);} }
-}
-
-
-// ********************************************** //
-// ********************************************** //
-
-psd10_rms::psd10_rms( lcurve inlc , vec fqL ) {
-
-	// ----------- initial parameters ------------ //
-	n		=	inlc.len;
-	dt		=	inlc.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	mean		=	inlc.demean();
-	mean2	=	mean*mean;
-	int		i;
-	for( i=0 ; i<n ; i++ ){
-		lc[i]	=	inlc.lc[i];
-		lce[i]	=	inlc.lce[i]*inlc.lce[i];
-		t[i]	=	inlc.t[i];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 1 );
-	// ------------------------------------------ //
-}
-
-psd10_rms::~psd10_rms() {}
-
-void psd10_rms::_C( vec p ){
-	int	i,j,k; double d; for(i=0;i<npar;i++){p[i]=pow(10,p[i])*mean2;}
-	for( i=0 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,i); }
-		d += f1*p[0]*Cfq2[0](i,i) + f2*p[nfq-1]*Cfq2[2](i,i);
-		d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[k]*Cfq[k](i,j); }
-		d += f1*p[0]*Cfq2[0](i,j) + f2*p[nfq-1]*Cfq2[2](i,j);
-		C(i,j) = d;}
-	}
-}
-
-void psd10_rms::_dC( vec p , int k ){
-	int	i,j; double d = log(10)*pow(10,p[k])*mean2;
-	for( i=0 ; i<n ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = d*Cfq[k](i,j); dC(j,i) = dC(i,j);} }
-}
-
-void psd10_rms::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
-		pars[i] += dpar[i];
-	}
-}
diff --git a/src/psdlag.cpp b/src/psdlag.cpp
deleted file mode 100644
index 79f68fd..0000000
--- a/src/psdlag.cpp
+++ /dev/null
@@ -1,252 +0,0 @@
-/*
- * psdlag.cpp
- *
- *  Created on: Jun 1, 2013
- *      Author: azoghbi
- */
-
-#include "inc/psdlag.hpp"
-
-psdlag::psdlag( lcurve lc1, lcurve lc2 , vec fqL ) {
-
-	// ----------- initial parameters ------------ //
-	n1		=	lc1.len;
-	n		=	n1 + lc2.len;
-	dt		=	lc1.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	lc1.demean(); lc2.demean();
-	int		i;
-	for( i=0 ; i<n1 ; i++ ){
-		lc[i]	=	lc1.lc[i];
-		lce[i]	=	lc1.lce[i]*lc1.lce[i];
-		t[i]	=	lc1.t[i];
-	}
-	for( i=n1 ; i<n ; i++ ){
-		lc[i]	=	lc2.lc[i-n1];
-		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
-		t[i]	=	lc2.t[i-n1];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 4 );
-	// ------------------------------------------ //
-
-	// --------- constants and indices ---------- //
-	ip1.setlength(nfq); ip2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
-	for( i=0 ; i<nfq ; i++ ){
-		ip1[i]	=	i;
-		ip2[i]	=	i+nfq;
-		icx[i]	=	i+2*nfq;
-		iphi[i]	=	i+3*nfq;
-	}
-	// ------------------------------------------ //
-
-}
-
-psdlag::~psdlag() {}
-
-void psdlag::_C( vec p ){
-	int	i,j,k; double d;
-	for( i=0 ; i<n1 ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,j); } C(i,j) = d;}
-	}
-	for( i=n1 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,j); } C(i,j) = d;}
-		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
-			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
-	}
-}
-
-void psdlag::_dC( vec p , int k ){
-	if( k<nfq ){			_dP1( p , k );
-	}else if( k<2*nfq ){	_dP2( p , k%nfq );
-	}else if( k<3*nfq ){	_dCx( p , k%nfq );
-	}else{					_dPhi( p , k%nfq );}
-}
-
-void psdlag::_dP1( vec p , int k ){
-	int	i,j;
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag::_dP2( vec p , int k ){
-	int	i,j;
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);}
-		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag::_dCx( vec p , int k ){
-	int	i,j; double phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag::_dPhi( vec p , int k ){
-	int	i,j; double cx=p[icx[k]],phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
-		pars[i] += dpar[i]/((n<10)?10:1);
-	}
-}
-
-
-void psdlag::print_pars( vec& pars , vec& errs ){
-	for( int i=0 ; i<nfq ; i++ ){
-		if( pars[i+2*nfq]<0 ){ pars[i+2*nfq]*=-1; pars[i+3*nfq] += M_PI;}
-		while( pars[i+3*nfq] > M_PI ){ pars[i+3*nfq] -= 2*M_PI; }
-		while( pars[i+3*nfq] <-M_PI ){ pars[i+3*nfq] += 2*M_PI; }
-		mod::print_pars( pars , errs );
-	}
-}
-
-// ++++++++++++++++++++++++++++++++++++++++++++++++++++++ //
-
-psdlag10::psdlag10( lcurve lc1, lcurve lc2 , vec fqL ) {
-
-	// ----------- initial parameters ------------ //
-	n1		=	lc1.len;
-	n		=	n1 + lc2.len;
-	dt		=	lc1.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	lc1.demean(); lc2.demean();
-	int		i;
-	for( i=0 ; i<n1 ; i++ ){
-		lc[i]	=	lc1.lc[i];
-		lce[i]	=	lc1.lce[i]*lc1.lce[i];
-		t[i]	=	lc1.t[i];
-	}
-	for( i=n1 ; i<n ; i++ ){
-		lc[i]	=	lc2.lc[i-n1];
-		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
-		t[i]	=	lc2.t[i-n1];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 4 );
-	// ------------------------------------------ //
-
-	// --------- constants and indices ---------- //
-	ip1.setlength(nfq); ip2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
-	for( i=0 ; i<nfq ; i++ ){
-		ip1[i]	=	i;
-		ip2[i]	=	i+nfq;
-		icx[i]	=	i+2*nfq;
-		iphi[i]	=	i+3*nfq;
-	}
-	// ------------------------------------------ //
-
-}
-
-psdlag10::~psdlag10() {}
-
-void psdlag10::_C( vec p ){
-	int	i,j,k; double d; for(i=0;i<3*nfq;i++){ p[i] = pow(10,p[i]);}
-	for( i=0 ; i<n1 ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,j); } C(i,j) = d;}
-	}
-	for( i=n1 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,j); } C(i,j) = d;}
-		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
-			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
-	}
-}
-
-void psdlag10::_dC( vec p , int k ){
-	if( k<nfq ){			_dP1( p , k );
-	}else if( k<2*nfq ){	_dP2( p , k%nfq );
-	}else if( k<3*nfq ){	_dCx( p , k%nfq );
-	}else{					_dPhi( p , k%nfq );}
-}
-
-void psdlag10::_dP1( vec p , int k ){
-	int	i,j; double p1 = log(10)*pow(10,p[k]);
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = p1 * Cfq[k](i,j); dC(j,i) = dC(i,j);} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag10::_dP2( vec p , int k ){
-	int	i,j; double p2 = log(10) * pow(10,p[ip2[k]]);
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = p2 * Cfq[k](i,j); dC(j,i) = dC(i,j);}
-		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag10::_dCx( vec p , int k ){
-	int	i,j; double phi=p[iphi[k]], cx=log(10)*pow(10,p[icx[k]]);
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag10::_dPhi( vec p , int k ){
-	int	i,j; double cx=pow(10,p[icx[k]]),phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag10::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
-		//pars[i] += dpar[i];
-		pars[i] += dpar[i]/((n<5)?10:1);
-	}
-}
-
-
-void psdlag10::print_pars( vec& pars , vec& errs ){
-	for( int i=0 ; i<nfq ; i++ ){
-		while( pars[i+3*nfq] > M_PI ){ pars[i+3*nfq] -= 2*M_PI; }
-		while( pars[i+3*nfq] <-M_PI ){ pars[i+3*nfq] += 2*M_PI; }
-	}
-	mod::print_pars( pars , errs );
-}
-
-void psdlag10::what_pars( int& ip1 , int& ip2 ){
-	ip1 = 3*nfq; ip2 = 4*nfq;
-}
diff --git a/src/psdlag_rms.cpp b/src/psdlag_rms.cpp
deleted file mode 100644
index 1f78912..0000000
--- a/src/psdlag_rms.cpp
+++ /dev/null
@@ -1,255 +0,0 @@
-/*
- * psdlag.cpp
- *
- *  Created on: Jun 1, 2013
- *      Author: azoghbi
- */
-
-#include "inc/psdlag_rms.hpp"
-
-psdlag_rms::psdlag_rms( lcurve lc1, lcurve lc2 , vec fqL ) {
-
-	// ----------- initial parameters ------------ //
-	n1		=	lc1.len;
-	n		=	n1 + lc2.len;
-	dt		=	lc1.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	mean1	=	lc1.demean(); mean2	=	lc2.demean();
-	mean1sq	=	mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
-	int		i;
-	for( i=0 ; i<n1 ; i++ ){
-		lc[i]	=	lc1.lc[i];
-		lce[i]	=	lc1.lce[i]*lc1.lce[i];
-		t[i]	=	lc1.t[i];
-	}
-	for( i=n1 ; i<n ; i++ ){
-		lc[i]	=	lc2.lc[i-n1];
-		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
-		t[i]	=	lc2.t[i-n1];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 4 );
-	// ------------------------------------------ //
-
-	// --------- constants and indices ---------- //
-	ip1.setlength(nfq); ip2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
-	for( i=0 ; i<nfq ; i++ ){
-		ip1[i]	=	i;
-		ip2[i]	=	i+nfq;
-		icx[i]	=	i+2*nfq;
-		iphi[i]	=	i+3*nfq;
-	}
-	// ------------------------------------------ //
-
-}
-
-psdlag_rms::~psdlag_rms() {}
-
-void psdlag_rms::_C( vec p ){
-	int	i,j,k; double d;
-	for( i=0 ; i<nfq ; i++ ){ p[ip1[i]] *= mean1sq; p[ip2[i]] *= mean2sq; p[icx[i]] *= mean12;}
-	for( i=0 ; i<n1 ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,j); } C(i,j) = d;}
-	}
-	for( i=n1 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,j); } C(i,j) = d;}
-		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
-			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
-	}
-}
-
-void psdlag_rms::_dC( vec p , int k ){
-	if( k<nfq ){			_dP1( p , k );
-	}else if( k<2*nfq ){	_dP2( p , k%nfq );
-	}else if( k<3*nfq ){	_dCx( p , k%nfq );
-	}else{					_dPhi( p , k%nfq );}
-}
-
-void psdlag_rms::_dP1( vec p , int k ){
-	int	i,j;
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag_rms::_dP2( vec p , int k ){
-	int	i,j;
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = Cfq[k](i,j); dC(j,i) = dC(i,j);}
-		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag_rms::_dCx( vec p , int k ){
-	int	i,j; double phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag_rms::_dPhi( vec p , int k ){
-	int	i,j; double cx=p[icx[k]]*mean12,phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag_rms::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
-		pars[i] += dpar[i]/((n<10)?10:1);
-	}
-}
-
-
-void psdlag_rms::print_pars( vec& pars , vec& errs ){
-	for( int i=0 ; i<nfq ; i++ ){
-		if( pars[i+2*nfq]<0 ){ pars[i+2*nfq]*=-1; pars[i+3*nfq] += M_PI;}
-		while( pars[i+3*nfq] > M_PI ){ pars[i+3*nfq] -= 2*M_PI; }
-		while( pars[i+3*nfq] <-M_PI ){ pars[i+3*nfq] += 2*M_PI; }
-		mod::print_pars( pars , errs );
-	}
-}
-
-// ++++++++++++++++++++++++++++++++++++++++++++++++++++++ //
-
-psdlag10_rms::psdlag10_rms( lcurve lc1, lcurve lc2 , vec fqL ) {
-
-	// ----------- initial parameters ------------ //
-	n1		=	lc1.len;
-	n		=	n1 + lc2.len;
-	dt		=	lc1.dt;
-	// ------------------------------------------ //
-
-
-	// ----------- light curve setup ------------ //
-	setlc();
-	mean1	=	lc1.demean(); mean2	=	lc2.demean();
-	mean1sq	=	mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
-	int		i;
-	for( i=0 ; i<n1 ; i++ ){
-		lc[i]	=	lc1.lc[i];
-		lce[i]	=	lc1.lce[i]*lc1.lce[i];
-		t[i]	=	lc1.t[i];
-	}
-	for( i=n1 ; i<n ; i++ ){
-		lc[i]	=	lc2.lc[i-n1];
-		lce[i]	=	lc2.lce[i-n1]*lc2.lce[i-n1];
-		t[i]	=	lc2.t[i-n1];
-	}
-	// ------------------------------------------ //
-
-	// ----------- parent initilizer ------------ //
-	init( fqL , 4 );
-	// ------------------------------------------ //
-
-	// --------- constants and indices ---------- //
-	ip1.setlength(nfq); ip2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
-	for( i=0 ; i<nfq ; i++ ){
-		ip1[i]	=	i;
-		ip2[i]	=	i+nfq;
-		icx[i]	=	i+2*nfq;
-		iphi[i]	=	i+3*nfq;
-	}
-	// ------------------------------------------ //
-
-}
-
-psdlag10_rms::~psdlag10_rms() {}
-
-void psdlag10_rms::_C( vec p ){
-	int	i,j,k; double d; for(i=0;i<nfq;i++){ p[ip1[i]] = pow(10,p[ip1[i]])*mean1sq; p[ip2[i]] = pow(10,p[ip2[i]])*mean2sq; p[icx[i]] = pow(10,p[icx[i]])*mean12;}
-	for( i=0 ; i<n1 ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip1[k]]*Cfq[k](i,j); } C(i,j) = d;}
-	}
-	for( i=n1 ; i<n ; i++ ){
-		d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
-		for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p[ip2[k]]*Cfq[k](i,j); } C(i,j) = d;}
-		for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
-			d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
-	}
-}
-
-void psdlag10_rms::_dC( vec p , int k ){
-	if( k<nfq ){			_dP1( p , k );
-	}else if( k<2*nfq ){	_dP2( p , k%nfq );
-	}else if( k<3*nfq ){	_dCx( p , k%nfq );
-	}else{					_dPhi( p , k%nfq );}
-}
-
-void psdlag10_rms::_dP1( vec p , int k ){
-	int	i,j; double p1 = log(10)*pow(10,p[k])*mean1sq;
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = p1 * Cfq[k](i,j); dC(j,i) = dC(i,j);} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag10_rms::_dP2( vec p , int k ){
-	int	i,j; double p2 = log(10) * pow(10,p[ip2[k]])*mean2sq;
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = p2 * Cfq[k](i,j); dC(j,i) = dC(i,j);}
-		for( j=0 ; j<n1 ; j++ ){ dC(i,j) = 0; dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag10_rms::_dCx( vec p , int k ){
-	int	i,j; double phi=p[iphi[k]], cx=log(10)*pow(10,p[icx[k]])*mean12;
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag10_rms::_dPhi( vec p , int k ){
-	int	i,j; double cx=pow(10,p[icx[k]])*mean12,phi=p[iphi[k]];
-	for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
-	for( i=n1 ; i<n ; i++ ){
-		for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
-		for( j=0 ; j<n1 ; j++ ){
-			dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
-	}
-}
-
-void psdlag10_rms::step_pars( int n , vec& dpar , vec& pars ){
-	for( int i=0 ; i<npar ; i++ ){
-		if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
-		//pars[i] += dpar[i];
-		pars[i] += dpar[i]/((n<5)?10:1);
-	}
-}
-
-
-void psdlag10_rms::print_pars( vec& pars , vec& errs ){
-	for( int i=0 ; i<nfq ; i++ ){
-		while( pars[i+3*nfq] > M_PI ){ pars[i+3*nfq] -= 2*M_PI; }
-		while( pars[i+3*nfq] <-M_PI ){ pars[i+3*nfq] += 2*M_PI; }
-	}
-	mod::print_pars( pars , errs );
-}
-
-void psdlag10_rms::what_pars( int& ip1 , int& ip2 ){
-	ip1 = 3*nfq; ip2 = 4*nfq;
-}