/*
*C-code for generation of single-particle basis and many-body basis.
*Contain 2 functions, one called genSpBasis that generates the 
*single-particle basis and another called genMpBasis that 
*generates the many-body basis.
*/

#include <stdio.h> 
#include <stdlib.h>
#include <string.h>
#include <math.h>

typedef struct basis {
	int** basis;
	int noStates;
	int maxNoStates;
}basis;

typedef struct mpState {
	int* states;
	int noParticles;
} mpState;

basis* genSpBasis(int maxE);
void genMpBasis(int K, basis* spBasis, int maxE, int M, int P, int e, int m, int l,
		mpState mpState, basis* mpBasis, int position, int lastSpState);
int initBasis(basis*, int, int);
int addState(basis*, int *);

int main() {
	char e[256];
	char noParticles[256];
	char M[256];
	char P[256];
	int maxESpStates[] = {0, 0, 1, 2, 3, 4, 5, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 33, 36};

	printf("Enter maximum energy over the energy groundstate for the system: ");
	gets(e);

	printf("Enter number of particles for the system: ");
	gets(noParticles);

	printf("Enter parity (0 for positive parity and 1 for negative parity) for the system: ");
	gets(P);

	printf("Enter magnetic projection (0 or 1) for the system: ");
	gets(M);

	basis* spBasis;
	basis* mpBasis;

	int maxE = atoi(e) + maxESpStates[atoi(noParticles)-1];
	spBasis = genSpBasis(maxE - maxESpStates[atoi(noParticles)-2]);
	
	mpBasis = malloc(sizeof(basis));
	initBasis(mpBasis, 100, atoi(noParticles));

	mpState mpState;
	mpState.noParticles = atoi(noParticles);
	mpState.states = malloc(atoi(noParticles) * sizeof(int));

	genMpBasis(atoi(noParticles), spBasis, maxE, atoi(M), atoi(P), 0, 0, 0,
	 		     mpState, mpBasis, 0, 0);
	
	// FILE *fp;	

	// Print out mpBasis to file "mpBasis.txt"
	/*
	fp = fopen("mpBasis.txt", "w");
	fprintf(fp, "Number of states in mp-basis: %d\n", mpBasis->noStates);
	for (int i = 0; i < mpBasis->noStates; i++) {
		for (int j = 0; j < atoi(noParticles); j++) {
			fprintf(fp, "%d ", mpBasis->basis[i][j]);
		}
		fprintf(fp, "\n");
	}
	*/

	// Print out spBasis to file "spBasis.txt"	
	/*
	fp = fopen("spBasis.txt", "w");
	fprintf(fp, "Number of states in sp-basis: %d\n", spBasis->noStates);
	for (int i = 0; i < spBasis->noStates; i++) {
		for (int j = 0; j < atoi(noParticles); j++) {
			fprintf(fp, "%d ", spBasis->basis[i][j]);
		}
		fprintf(fp, "\n");
	}
	*/

	// Print out the single particle states
	/*
	for (int i = 0; i < spBasis->noStates; i++) {
		for (int j = 0; j < 4; j++) {
			printf("%d ", spBasis->basis[i][j]);
		}
		printf("\n");
	} 
	*/
	// Print out the many-particle basis
	/*for (int i = 0; i < mpBasis->noStates; i++) {
		for (int j = 0; j < atoi(noParticles); j++) {
			printf(" %d ", mpBasis->basis[i][j]);
		}
		printf("\n");
	} */

	printf("Number of single-particle states in single-particle basis: %d\n", spBasis->noStates);
	printf("Number of many-body states in many-body basis%d\n", mpBasis->noStates);
	
	//fclose(fp);
	return 0;
}

basis* genSpBasis(int maxE) {
	basis* spBasis = malloc(sizeof(basis));
	initBasis(spBasis, 500, 4);

	for (int n = 0; n <= maxE; n++) {
		for (int l = n % 2; l <= n; l += 2) {
			for (int m = -(2 * l + 1); m <= (2 * l + 1); m += 2) {
				int* spState = malloc(4 * sizeof(int));
				spState[0] = n;       
        			spState[1] = l;       
        			spState[2] = 2*l+1;       
				spState[3] = m;

				addState(spBasis, spState);
			}

			for (int m = -(2 * l - 1); m <= (2 * l - 1); m += 2) {
				int* spState = malloc(4 * sizeof(int));
				spState[0] = n;       
        			spState[1] = l;       
        			spState[2] = 2*l-1;       
				spState[3] = m;
 
				addState(spBasis, spState);
			}
		}
	}

	return spBasis;
}

void genMpBasis(int K, basis* spBasis, int maxE, int M, int P, int e, int m, int l, mpState mpState, 
		basis* mpBasis, int position, int lastSpState) {
	int min;

	if (position == 0) {
		min = 0;
	} else {
		min = lastSpState + 1;
	}

	for (int i = min; i < spBasis->noStates; i++) {
		mpState.states[position] = i;
		
		int newE = e + spBasis->basis[i][0];
		int newM = m + spBasis->basis[i][3];
		int newL = l + spBasis->basis[i][1];

		if (newE <= maxE) {
			if (position < K-1) {
				genMpBasis(K, spBasis, maxE, M, P, newE, newM, newL, mpState, mpBasis, 
					position+1, i);
			} else if (newM / 2 == M && newL % 2 == P) {
				int* mpStateP = malloc(mpState.noParticles * sizeof(int));			
				for (int i = 0; i < mpState.noParticles; i++) {
					mpStateP[i] = mpState.states[i];	
				}
				addState(mpBasis, mpStateP);
			} 	
		} else {
			break;
		}
	}
}


int initBasis(basis* basis, int initialSize, int sizeX) {
	basis->basis = malloc(initialSize * sizeof(int*));
	for (int i = 0; i < initialSize; i++) {
		basis->basis[i] = malloc(sizeX * sizeof(int));
	}
	basis->noStates = 0;
	basis->maxNoStates = initialSize;

	return 1;
}

int addState(basis* basis, int *state) {
	if (basis->noStates == basis->maxNoStates) {
		basis->maxNoStates *= 2;
		basis->basis = realloc(basis->basis, basis->maxNoStates * sizeof(int*));
	}
	basis->basis[basis->noStates++] = state;	
}