/* simulate.c
   ----------
   Calculates meas array info for all cell sites, using the partition model calculated by main.c.
   Generates 8 .MAT files containing 'double' arrays which need to be loaded and casted by hande.
   
   Uses:
   Matlab External Interface Library - to load and save .MAT data
*/

#include "include\mat.h"
#include "stdlib.h"
#include "conio.h"
#include "math.h"
#include "string.h"

#define PI 3.1415926535897932384626433832795

/* file locations */
#pragma comment(lib, "C:\\Appl\\MatlabR14\\extern\\lib\\win32\\microsoft\\msvc60\\libmat.lib")
#pragma comment(lib, "C:\\Appl\\MatlabR14\\extern\\lib\\win32\\microsoft\\msvc60\\libmx.lib")
#pragma comment(lib, "C:\\Appl\\MatlabR14\\extern\\lib\\win32\\microsoft\\msvc60\\libeng.lib")
#pragma comment(lib, "libgsl.lib")
#define CELL_INFO_FILE "c:\\temp\\3065\\cell_info.mat"

/* partition coefficients */
#define DISTANCE_COEFF (-0.00090122)
#define BLOCKS_COEFF   (-4.0490)
#define EXCESS_COEFF   (-100.3736)
#define ANGLE_COEFF    (-4.1160)
#define DIFF_COEFF     (-1.0625)
#define AZIMUTH_COEFF  (9.0302)
#define ANGULAR_SPREAD (120)
#define PR_FLOOR       (-126)

void getCellInfo(MATFile *pmat, char *cell_name, short **terrain, double *azimuth, double *eirp, double *BS_x, double *BS_y, double *freq, double *cellSize, double *height, int *w, int *h);
int numBlocks(short *terrain, int tx, int ty, int tz, int rx, int ry, int rz, int h, int cellSize, double wavelength, double *diff_score);
double *simulateCell(short *terrain, double azimuth, double eirp, double BS_x, double BS_y, double freq, double cellSize, double height, int w, int h);
void save(MATFile *pmat, char *cell_name, double *sim, int w, int h, int meas_exists);
void report(char *message, int num, int max_count);

int main ( int argc, char **argv ) {
	/* cell info data */
	MATFile *pmat;
	int wa=0, wb=0, wc=0, wd=0, we, wf, wg, wh;
	int ha=0, hb=0, hc=0, hd=0, he, hf, hg, hh;

	/* cell data */
	short *terrain;
	double azimuth;
	double eirp;
	double BS_x;
	double BS_y;
	double freq;
	double cellSize;
	double height;

	/* simulation data */
	double *sim_a;
 	double *sim_b;
	double *sim_c;
	double *sim_d;
	double *sim_e;
	double *sim_f;
	double *sim_g;
	double *sim_h;

	/* load MATLAB file */
	printf("Loading Cell_Info.mat...");
	pmat = matOpen(CELL_INFO_FILE, "r");
	if (!pmat) {
		printf("Couldn't load file\n");
		getch();
		return 0;
	}
	printf("done\n");

	/* Simulate Cell A */
	printf("Simulating cells\n");
	getCellInfo(pmat, "Cell_Info_A", &terrain, &azimuth, &eirp, &BS_x, &BS_y, &freq, &cellSize, &height, &wa, &ha);
	printf("Simulating Cell A");
	sim_a = simulateCell(terrain, azimuth, eirp, BS_x, BS_y, freq, cellSize, height, wa, ha);
	printf("done\n");

	/* Simulate Cell B */
	getCellInfo(pmat, "Cell_Info_B", &terrain, &azimuth, &eirp, &BS_x, &BS_y, &freq, &cellSize, &height, &wb, &hb);
	printf("Simulating Cell B...");
	sim_b = simulateCell(terrain, azimuth, eirp, BS_x, BS_y, freq, cellSize, height, wb, hb);
	printf("done\n");

	/* Simulate Cell C */
	getCellInfo(pmat, "Cell_Info_C", &terrain, &azimuth, &eirp, &BS_x, &BS_y, &freq, &cellSize, &height, &wc, &hc);
	printf("Simulating Cell C...");
	sim_c = simulateCell(terrain, azimuth, eirp, BS_x, BS_y, freq, cellSize, height, wc, hc);
	printf("done\n");

	/* Simulate Cell D */
	printf("Simulating Cell D...");
	getCellInfo(pmat, "Cell_Info_D", &terrain, &azimuth, &eirp, &BS_x, &BS_y, &freq, &cellSize, &height, &wd, &hd);
	sim_d = simulateCell(terrain, azimuth, eirp, BS_x, BS_y, freq, cellSize, height, wd, hd);
	printf("done\n");

	/* Simulate Cell E */
	getCellInfo(pmat, "Cell_Info_E", &terrain, &azimuth, &eirp, &BS_x, &BS_y, &freq, &cellSize, &height, &we, &he);
	printf("Simulating Cell E...");
	sim_e = simulateCell(terrain, azimuth, eirp, BS_x, BS_y, freq, cellSize, height, we, he);
	printf("done\n");
	
	/* Simulate Cell F */
	getCellInfo(pmat, "Cell_Info_F", &terrain, &azimuth, &eirp, &BS_x, &BS_y, &freq, &cellSize, &height, &wf, &hf);
	printf("Simulating Cell F...");
	sim_f = simulateCell(terrain, azimuth, eirp, BS_x, BS_y, freq, cellSize, height, wf, hf);
	printf("done\n");
	
	/* Simulate Cell G */
	getCellInfo(pmat, "Cell_Info_G", &terrain, &azimuth, &eirp, &BS_x, &BS_y, &freq, &cellSize, &height, &wg, &hg);
	printf("Simulating Cell G...");
	sim_g = simulateCell(terrain, azimuth, eirp, BS_x, BS_y, freq, cellSize, height, wg, hg);
	printf("done\n");
	
	/* Simulate Cell H */
	getCellInfo(pmat, "Cell_Info_H", &terrain, &azimuth, &eirp, &BS_x, &BS_y, &freq, &cellSize, &height, &wh, &hh);
	printf("Simulating Cell H...");
	sim_h = simulateCell(terrain, azimuth, eirp, BS_x, BS_y, freq, cellSize, height, wh, hh);

	printf("done\n");

	/* Write data back to MAT files */
	printf("Saving MAT file...");
	save(pmat, "Cell_Info_A", sim_a, wa, ha, 0);
	save(pmat, "Cell_Info_B", sim_b, wb, hb, 0);
	save(pmat, "Cell_Info_C", sim_c, wc, hc, 0);
	save(pmat, "Cell_Info_D", sim_d, wd, hd, 0);
	save(pmat, "Cell_Info_E", sim_e, we, he, 0);
	save(pmat, "Cell_Info_F", sim_f, wf, hf, 0);
	save(pmat, "Cell_Info_G", sim_g, wg, hg, 0);
	save(pmat, "Cell_Info_H", sim_h, wh, hh, 0);
	printf("done\n");

	/* close MAT file */
	matClose(pmat);

	printf("\nAll finished\n");
	getch();

	return 0;
}

/* getCelInfo - extracts useful information from Cell_Info file */
void getCellInfo(MATFile *pmat, char *cell_name, short **terrain, double *azimuth, double *eirp, double *BS_x, double *BS_y, double *freq, double *cellSize, double *height, int *w, int *h) {
	/* matlab objects */
    mxArray *cell_info;
	mxArray *terrain_array;
	mxArray *cell_site;
	mxArray *azimuth_array;
	mxArray *eirp_array;
	mxArray *BS_x_array;
	mxArray *BS_y_array;
	mxArray *freq_array;
	mxArray *cellSize_array;
	mxArray *height_array;
	const int *terrain_dim;

	/* get matlab arrays */
	cell_info = matGetVariable(pmat, cell_name);
	terrain_array = mxGetField(cell_info, 0, "terrain");
	cell_site = mxGetField(cell_info, 0, "cellSite");
	azimuth_array = mxGetField(cell_site, 0, "azimuth");
	eirp_array = mxGetField(cell_site, 0, "eirp");
	BS_x_array = mxGetField(cell_site, 0, "BSx");
	BS_y_array = mxGetField(cell_site, 0, "BSy");
	freq_array = mxGetField(cell_site, 0, "freq");
	cellSize_array = mxGetField(cell_site, 0, "cellSize");
	height_array = mxGetField(cell_site, 0, "height");
	terrain_dim = mxGetDimensions(terrain_array);

	/* extract and export data */
	*w = terrain_dim[1];
	*h = terrain_dim[0];
	*terrain =(short *)mxGetPr(terrain_array);
	*azimuth = mxGetScalar(azimuth_array);
	*eirp = mxGetScalar(eirp_array);
	*BS_x = mxGetScalar(BS_x_array);
	*BS_y = mxGetScalar(BS_y_array);
	*freq = mxGetScalar(freq_array);
	*cellSize = mxGetScalar(cellSize_array);
	*height = mxGetScalar(height_array);
}

/* numBlocks - calculates how much space between (tx, ty, tz) & (rx, ry, rz) is blocked */
int numBlocks(short *terrain, int tx, int ty, int tz, int rx, int ry, int rz, int h, int cellSize, double wavelength, double *diff_score) {
	double distance2d;
	double angle2d;
	double zslope;
	double expected_z;
	double actual_z;
	
	int thisx, thisy, r;
	int blocks=0;

	double ob_h, d1, d2, v;

	int dx = rx - tx;
	int dy = ry - ty;
	int dz = rz - tz;
	
	distance2d = sqrt(dx*dx+dy*dy);
	angle2d = atan2(dy, dx);
	zslope = (double)dz / distance2d;
	*diff_score = 0;

	/* iterate from transmitter to receiver, checking terrain */
	for (r=0; r<(int)distance2d; r+=cellSize) {
		thisx = (int)(cos(angle2d)*r)+tx;
		thisy = (int)(sin(angle2d)*r)+ty;

		/* calculate expected z as a line of slope "zslope" */
		expected_z = zslope*r + tz;
		actual_z = (double)terrain[thisx*h+thisy];

		/* if the terrain is heigher then the expected z, obstruction */
		if (actual_z > expected_z) {

			/* compute diffraction score if h << distance */
			if ( abs((int)actual_z - tz) < (int)(distance2d / 100)) {
				ob_h = actual_z - rz;
				d1 = r;
				d2 = distance2d - r;
				v = ob_h*(d1+d2)/(wavelength*d1*d2);
				*diff_score = v;
			}
			blocks = 1;
		}
	}

	return blocks;
}

/* simulateCell - calculates the "meas" array for a given set of cell data */
double *simulateCell(short *terrain, double azimuth, double eirp, double BS_x, double BS_y, double freq, double cellSize, double height, int w, int h) {
	double *sim = malloc(w*h*sizeof(double));
	int tx = (int)(BS_x-1);
	int ty = (int)(BS_y-1);
	int tz = terrain[tx*h+ty] + (int)height;
	int x, y, z;
	int dx, dy, dz;
	int blocks;
	int in_azimuth_angle;
	double d, r;
	double diff_score;
	double angle;
	double wavelength = 300000000/(freq*1000000);

	/* azimuth angle data */
	double azimuth_rad;
	double azimuth_x;
	double azimuth_y;
	double angle2d_rad;
	double angle2d_x;
	double angle2d_y;
	double angular_separation;

	/* prelim azimuth calculations */
	azimuth_rad = (-azimuth + 90)*PI/180;
	azimuth_x = cos(azimuth_rad);
	azimuth_y = sin(azimuth_rad);

	/* calculate "meas" for each element in terrain array */
	for (y=0; y<h; y++) {
		for(x=0; x<w; x++) {
			z = terrain[x*h+y];
			dx = (int)((x - tx)*cellSize);
			dy = (int)((y - ty)*cellSize);
			dz = z - tz;
			d = sqrt(dx*dx + dy*dy + dz*dz);
			r = sqrt(dx*dx + dy*dy);
			angle = atan2(dz, r)*180/PI;
			blocks = numBlocks(terrain, tx, ty, tz, x, y, z, h, (int)cellSize, wavelength, &diff_score);

			/* azimuth calcuations */
			angle2d_rad = atan2(dy, dx);
			angle2d_x = cos(angle2d_rad);
			angle2d_y = sin(angle2d_rad);
			angular_separation = acos((azimuth_x * angle2d_x) + (azimuth_y * angle2d_y)) * 180 / PI;
			if (angular_separation < (ANGULAR_SPREAD / 2)) in_azimuth_angle = 1;
			else in_azimuth_angle = 0;

			/* the linear partition model */
			sim[x*h+y] = DISTANCE_COEFF*d;// + BLOCKS_COEFF*(double)blocks + ANGLE_COEFF*angle + DIFF_COEFF*diff_score + EXCESS_COEFF;
			sim[x*h+y] += BLOCKS_COEFF*(double)blocks;
			sim[x*h+y] += ANGLE_COEFF*angle;
			sim[x*h+y] += DIFF_COEFF*diff_score;
			sim[x*h+y] += AZIMUTH_COEFF*(double)in_azimuth_angle;
			sim[x*h+y] += EXCESS_COEFF;

			/* cap off at noise floor */
			if (sim[x*h+y] < PR_FLOOR) sim[x*h+y] = PR_FLOOR;
		}
		if (!(y % 100)) printf(".");
	}


	return sim;
}

void save(MATFile *pmat, char *cell_name, double *sim, int w, int h, int meas_exists) {
	MATFile *out;
	mxArray *meas_info;
	char cell_file[256];

	/* transfer data to a MATLAB array */
	meas_info = mxCreateDoubleMatrix(h, w, mxREAL);
	mxSetPr(meas_info, sim);
	
	/* generate file name */
	memcpy(cell_file, cell_name, strlen(cell_name));
	memcpy(cell_file+strlen(cell_name), "_meas.mat", 9);
	cell_file[strlen(cell_name)+9] = 0;

	/* create and save file */
	out = matOpen(cell_file, "w");
	matPutVariableAsGlobal(out, "meas", meas_info);
	matClose(out);
}

/* report - little "helper" assertion routine */
void report(char *message, int num, int max_count) {
	static count=0;

	if (count == 0) printf("\n");

	if (count < max_count) {
		printf("%s(%d): %d\n", message, count, num);
		count++;
	}
	
}