#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "catalog.c"
#define PI 3.1415926535897932384626433832795028842

double sphericalDistance(double ra1, double dec1, double ra2, double dec2) {
	return acos(sin(dec1) * sin(dec2) + cos(dec1) * cos(dec2) * cos(ra2 - ra1));
}

double starScore(double rao, double deco, double ras, double decs, double lums) {
	#define METHOD 2
	#define LUMEXP 2
	
	#if     LUMEXP == 0
	#define LUMW 1
	#elif   LUMEXP == 1
	#define LUMW lums
	#elif LUMEXP == 2
	#define LUMW (lums * lums)
	#elif LUMEXP == 3
	#define LUMW (lums * lums * lums)
	#endif
	
	#if METHOD == 2
	#define SCO(d, l) (cos(dist) * l)
	#elif METHOD == 1
	#define SCO(d, l) ((dist < PI / 2) ? (cos(dist) * l) : 0)
	#endif
	
	double dist = sphericalDistance(rao, deco, ras, decs);
	return SCO(dist, LUMW); // (dist < PI / 2) ? (cos(dist) * lums * lums) : 0;
}

double totalScore(double ra, double dec, const double catalog[][3], int nstars) {
	double score = 0;
	for(int i = 0; i < nstars; i++)
		score += starScore(ra, dec, catalog[i][0], catalog[i][1], catalog[i][2]);
	return score;
}

void findMaximumScore(double *ra, double *dec, double *score, const double catalog[][3], int nstars) {
	double cra = PI;
	double cdec = 0;
	
	double prevscore, cscore, cgradra, dra, ddec, cgraddec, normalgrad;
	for(int i = 0; ; i++) {
		prevscore = cscore;
		cscore = totalScore(cra, cdec, catalog, nstars);
		if(prevscore > cscore) { // half steps back if the score is now lower than it was before
			cra -= dra / 2;
			cdec -= ddec / 2;
			cscore = totalScore(cra, cdec, catalog, nstars);
		}
		cgradra = totalScore(cra + 0.001, cdec, catalog, nstars) - cscore;
		cgraddec = totalScore(cra, cdec + 0.001, catalog, nstars) - cscore;
		
		if(i == 0) normalgrad = sqrt(cgradra * cgradra + cgraddec * cgraddec);
		
		dra = cgradra / normalgrad;
		ddec = cgraddec / normalgrad;
		/*dra = fabs(dra) > 0.5 ? 0.5 * (dra > 0 ? 1 : -1) : dra;
		ddec = fabs(ddec) > 0.5 ? 0.5 * (ddec > 0 ? 1 : -1) : ddec;*/
		
		cra += dra;
		cdec += ddec;
		
		cra = fmod(cra, 2 * PI);
		cdec = fmod(cdec - PI / 2, PI) + PI / 2;
		
		if(fabs(dra) < 0.001 && fabs(ddec) < 0.001) break;
	}
	
	*ra = cra;
	*dec = cdec;
	*score = cscore;
}

double makeScoreGrid(double rai, double deci, double raf, double decf, int nra, int ndec, const double catalog[][3], int nstars) {
	double maxsc = -1e+308;
	double minsc = 1e+308;
	double *scoregrid = (double *)malloc(nra * ndec * sizeof(double));
	
	for(int i = 0; i < ndec; i++) {
		for(int j = 0; j < nra; j++) {
			double rao = rai + j * (raf - rai) / (nra - 1);
			double deco = deci + i * (decf - deci) / (ndec - 1);
			double cscore = 0;
			for(int k = 0; k < nstars; k++) cscore += starScore(rao, deco, catalog[k][0], catalog[k][1], catalog[k][2]);
			if(cscore > maxsc) maxsc = cscore;
			if(cscore < minsc) minsc = cscore;
			scoregrid[i * nra + j] = cscore;
		}
	}
	
	for(int i = 0; i < ndec; i++) {
		for(int j = 0; j < nra; j++) {
			double sc = (scoregrid[i * nra + j] - minsc) / (maxsc - minsc);
			putchar(sc > 0.99 ? '!' : (sc > 0.01 ? '.' : ' '));
		}
		putchar('\n');
	}
}

int main() {
	int catalogsize = sizeof(catalog) / (3 * sizeof(double));
	// makeScoreGrid(0.0, -PI / 2, PI, PI / 2, 40, 20, catalog, catalogsize);
	
	// #define DOLUMINOSITYDISTRIBUTION
	#ifdef DOLUMINOSITYDISTRIBUTION
	double lumdist[10] = { 0 };
	int lumdistn[10] = { 0 };
	for(int i = 0; i < catalogsize; i++) {
		double lum = catalog[i][2];
		int lumn = (int)floor(-2.5 * log10(lum)) + 2;
		lumdist[lumn] += lum * lum; // related to the luminosity weighting in the starScore function
		lumdistn[lumn]++;
	}
	printf("Luminosity distribution:\n"
		"\tmag-2~-1-> %f : %i\n"
		"\tmag-1~0 -> %f : %i\n"
		"\tmag 0~1 -> %f : %i\n"
		"\tmag 1~2 -> %f : %i\n"
		"\tmag 2~3 -> %f : %i\n"
		"\tmag 3~4 -> %f : %i\n"
		"\tmag 4~5 -> %f : %i\n"
		"\tmag 5~6 -> %f : %i\n",
		lumdist[0], lumdistn[0], lumdist[1], lumdistn[1], lumdist[2], lumdistn[2], lumdist[3], lumdistn[3],
		lumdist[4], lumdistn[4], lumdist[5], lumdistn[5], lumdist[6], lumdistn[6], lumdist[7], lumdistn[7]
	);
	#endif
	
	double ra, dec, score;
	findMaximumScore(&ra, &dec, &score, catalog, catalogsize);
	makeScoreGrid(
		ra - 0.1, dec - 0.1,
		ra + 0.1, dec + 0.1,
		11, 11,
		catalog, catalogsize
	);
	printf("\t\t<td>\n\t\t\tRA %.1f&deg; Dec %.1f&deg;<br/>\n\t\t\t%.2f ~ %.2f\n\t\t</td>\n", ra * 180 / PI, dec * 180 / PI, score, totalScore(fmod(ra + PI, 2 * PI), -dec, catalog, catalogsize));
}