/*! \file Commons.cpp
 *
 *	DEVELOPMENT NOTE:
 *	The construction of common blocks requires some information
 *	that is stored in configuration objects and is needed to compute
 *	the allocation size of vectors and matrices. Currently, this
 *	information is passed as arguments to the constructors of the
 *	common blocks. A simpler and more logical way to operate is
 *	to design the constructors to take as arguments only pointers
 *	to the configuration objects. These, on their turn, need to
 *	expose methods to access the relevant data in read-only mode.
 */
#include <complex>

#ifndef INCLUDE_COMMONS_H
#include "../include/Commons.h"
#endif

using namespace std;

C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
	nlmmt = 2 * (l_max * (l_max + 2));
	nsph = ns;
	lm = l_max;

	rmi = new complex<double>*[lm];
	rei = new complex<double>*[lm];
	for (int ri = 0; ri < lm; ri++) {
		rmi[ri] = new complex<double>[nsph]();
		rei[ri] = new complex<double>[nsph]();
	}
	w = new complex<double>*[nlmmt];
	for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4]();
	vints = new complex<double>*[nsph];
	rc = new double*[nsph];
	nshl = new int[nsph]();
	iog = new int[nsph]();
	for (int vi = 0; vi < nsph; vi++) {
		rc[vi] = new double[_nshl[vi]]();
		vints[vi] = new complex<double>[16]();
		nshl[vi] = _nshl[vi];
		iog[vi] = _iog[vi];
	}
	fsas = new complex<double>[nsph]();
	sscs = new double[nsph]();
	sexs = new double[nsph]();
	sabs = new double[nsph]();
	sqscs = new double[nsph]();
	sqexs = new double[nsph]();
	sqabs = new double[nsph]();
	gcsv = new double[nsph]();
	rxx = new double[nsph]();
	ryy = new double[nsph]();
	rzz = new double[nsph]();
	ros = new double[nsph]();

	sas = new complex<double>**[nsph];
	for (int si = 0; si < nsph; si++) {
		sas[si] = new complex<double>*[2];
		sas[si][0] = new complex<double>[2]();
		sas[si][1] = new complex<double>[2]();
	}
}

C1::~C1() {
	delete[] rmi;
	delete[] rei;
	for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
	for (int vi = nsph - 1; vi > - 1; vi--) {
		delete[] rc[vi];
		delete[] vints[vi];
	}
	delete[] rc;
	delete[] vints;
	for (int si = nsph - 1; si > -1; si--) {
		delete[] sas[si][1];
		delete[] sas[si][0];
		delete[] sas[si];
	}
	delete[] sas;
	delete[] fsas;
	delete[] sscs;
	delete[] sexs;
	delete[] sabs;
	delete[] sqscs;
	delete[] sqexs;
	delete[] sqabs;
	delete[] gcsv;
	delete[] rxx;
	delete[] ryy;
	delete[] rzz;
	delete[] ros;
	delete[] iog;
	delete[] nshl;
}

C1_AddOns::C1_AddOns(C4 *c4) {
	nsph = c4->nsph;
	lmpo = c4->lmpo;
	nlemt = 2 * c4->nlem;
	vh = new complex<double>[(nsph * nsph - 1) * c4->litpo]();
	vj0 = new complex<double>[nsph * c4->lmtpo]();
	vj = new complex<double>[1](); // QUESTION: is 1 really enough for a general case?
	vyhj = new complex<double>[(nsph * nsph - 1) * c4->litpos]();
	vyj0 = new complex<double>[nsph * c4->lmtpos]();
	am0m = new complex<double>*[nlemt];
	for (int ai = 0; ai < nlemt; ai++) {
		am0m[ai] = new complex<double>[nlemt]();
	}
	vint = new complex<double>[16](); // QUESTION: is dimension 16 generally fixed?
	vintm = new complex<double>[16]();
	vintt = new complex<double>[16]();
	vints = new complex<double>*[nsph];
	for (int vi = 0; vi < nsph; vi++) vints[vi] = new complex<double>[16]();
	fsac = new complex<double>*[2];
	sac = new complex<double>*[2];
	fsacm = new complex<double>*[2];
	for (int fi = 0; fi < 2; fi++) {
		fsac[fi] = new complex<double>[2]();
		sac[fi] = new complex<double>[2]();
		fsacm[fi] = new complex<double>[2]();
	}
	scscp = new complex<double>[2]();
	ecscp = new complex<double>[2]();
	scscpm = new complex<double>[2]();
	ecscpm = new complex<double>[2]();
	allocate_vectors(c4);
	sscs = new double[nsph]();
	ecscm = new double[2]();
	scscm = new double[2]();
	scsc = new double[2]();
	ecsc = new double[2]();
}

C1_AddOns::~C1_AddOns() {
	delete[] sscs;
	delete[] vyj0;
	delete[] vyhj;
	for (int ii = lmpo - 1; ii > -1; ii--) delete[] ind3j[ii];
	delete[] ind3j;
	delete[] v3j0;
	delete[] vh;
	delete[] vj0;
	for (int ai = nlemt - 1; ai > -1; ai--) {
		delete[] am0m[ai];
	}
	delete am0m;
	delete[] vint;
	delete[] vintm;
	delete[] vintt;
	for (int vi = nsph - 1; vi > -1; vi--) delete[] vints[vi];
	delete[] vints;
	for (int fi = 1; fi > -1; fi--) {
		delete[] fsac[fi];
		delete[] sac[fi];
		delete[] fsacm[fi];
	}
	delete[] fsac;
	delete[] sac;
	delete[] fsacm;
	delete[] scscp;
	delete[] ecscp;
	delete[] scscpm;
	delete[] ecscpm;
	delete[] ecscm;
	delete[] scscm;
	delete[] scsc;
	delete[] ecsc;
}

void C1_AddOns::allocate_vectors(C4 *c4) {
	// This calculates the size of v3j0
	int lm = (c4->li > c4->le) ? c4->li : c4->le;
	const int nv3j = c4->nv3j;
	v3j0 = new double[nv3j]();
	ind3j = new int*[lmpo];
	for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
	// This calculates the size of vyhj
	int ivy = (nsph * nsph - 1) * c4->litpos;
	vyhj = new complex<double>[ivy]();
	// This calculates the size of vyj0
	ivy = c4->lmtpos * c4->nsph;
	vyj0 = new complex<double>[ivy]();
}

C2::C2(int ns, int nl, int npnt, int npntts) {
	nsph = ns;
	int max_n = (npnt > npntts) ? npnt : npntts;
	nhspo = 2 * max_n - 1;
	ris = new complex<double>[nhspo]();
	dlri = new complex<double>[nhspo]();
	vkt = new complex<double>[nsph]();
	dc0 = new complex<double>[nl]();
	vsz = new double[nsph]();
}

C2::~C2() {
	delete[] ris;
	delete[] dlri;
	delete[] vkt;
	delete[] dc0;
	delete[] vsz;
}

C3::C3() {
	tsas = new complex<double>*[2];
	tsas[0] = new complex<double>[2];
	tsas[1] = new complex<double>[2];
	tfsas = complex<double>(0.0, 0.0);
	gcs = 0.0;
	scs = 0.0;
	ecs = 0.0;
	acs = 0.0;
}

C3::~C3() {
	delete[] tsas[1];
	delete[] tsas[0];
	delete[] tsas;
}

C6::C6(int lmtpo) {
	rac3j = new double[lmtpo]();
}

C6::~C6() {
	delete[] rac3j;
}

C9::C9(int ndi, int nlem, int ndit, int nlemt) {
	gis = new complex<double>*[ndi];
	gls = new complex<double>*[ndi];
	for (int gi = 0; gi < ndi; gi++) {
		gis[gi] = new complex<double>[nlem]();
		gls[gi] = new complex<double>[nlem]();
	}
	sam = new complex<double>*[ndit];
	for (int si = 0; si < ndit; si++) sam[si] = new complex<double>[nlemt]();
	gis_size_0 = ndi;
	sam_size_0 = ndit;
}

C9::~C9() {
	for (int gi = gis_size_0 - 1; gi > -1; gi--) {
		delete[] gis[gi];
		delete[] gls[gi];
	}
	delete[] gis;
	delete[] gls;
	for (int si = sam_size_0 - 1; si > -1; si--) delete[] sam[si];
	delete[] sam;
}
