#include "header.h"
#include "proto.h"
#include "output.h"


/*
 * This file is highly data dependent.  Any change in the structure or
 * organization of the database will cause this stuff to stop working.
*/

typedef struct _vdata {
	int count;
	double scale;
	PTR ptr;
} VDATA;

VDATA * GetVarDataInit(PTR vraw, VARDATA *vardata, VDATA *v);
DATA GetVarDataValue(PTR vraw, VARDATA *vardata, VDATA *v, int offset);
double compute_bbr(int chan, double radiance, double temp);
void ff_emissivity(OSTRUCT **o);
void make_tesx();

double *tesw[3] = { 0, 0, 0 };
double bbr(double wn, double temp);

FIELD *
FindFakeField(char *name, LIST *tables)
{
	FIELD *f, *f0, *f1, *f2;
	FAKEFIELD *ff;
	LIST *list;
    char field_name[256], *p;

    strcpy(field_name, name);

    /**
    ** If this field name includes a dimension, get rid of it
    **/
    if ((p = strchr(field_name, '[')) != NULL) {
        *p = '\0';
    }

	if (!strcasecmp(field_name, "emissivity")) {
		if ((f0 = FindField("obs.scan_len", tables)) == NULL) {
			fprintf(stderr, "Unable to find field: obs.scan_len, needed "
							"by derived field: emissivity\n");
			exit(1);
		}
		if ((f1 = FindField("rad.target_temp", tables)) == NULL) {
			fprintf(stderr, "Unable to find field: rad.target_temp, needed "
							"by derived field: emissivity\n");
			exit(1);
		}
		if ((f2 = FindField("rad.cal_rad", tables)) == NULL) {
			fprintf(stderr, "Unable to find field: rad.cal_rad, needed "
							"by derived field: emissivity\n");
			exit(1);
		}

		f = (FIELD *)calloc(1, sizeof(FIELD));
		ff = (FAKEFIELD *)calloc(1, sizeof(FAKEFIELD));

		list = new_list();
		list_add(list, f0);
		list_add(list, f1);
		list_add(list, f2);
		ff->fields = (FIELD **)list->ptr;
		ff->nfields = list->number;
		free(list);
			
		ff->fptr = ff_emissivity;

		/*
		** This field looks and acts like var_data
		*/
		f->name = strdup("emissivity");
		f->dimension = ff->fields[2]->dimension;
		f->vardata = ff->fields[2]->vardata;
		f->iformat = REAL;
		f->size = 4;
		f->fakefield = ff;

		make_tesx();

		return(f);
	}
	return(NULL);
}


void
ff_emissivity(OSTRUCT **o)
{
	int i, n, end, scan_len;
	double  target_temp, result, rad;
	DATA data, dvar;
	PTR raw, var;
	VDATA vdata;
	/*
	** Next N ostructs' are dependent fields.  Calculate the value 
	** and output it.
	*/

	/*
	** Get scan_len
	*/
	raw = o[1]->slice->start_rec + o[1]->field->start;
	scan_len = (*ConvertData(raw, o[1]->field).str == '1' ? 1 : 2);

	/*
	** Get target_temp
	*/
	raw = o[2]->slice->start_rec + o[2]->field->start;
	target_temp = ConvertAndScaleData(raw, o[2]->field);

	/*
	** Get byte pointer to varfield data
	*/
	raw = o[3]->slice->start_rec + o[3]->field->start;
	data = ConvertData(raw, o[3]->field);

	if (data.i == -1) {
		printf("N/A");
		return;
	}

	var = GiveMeVarPtr(o[3]->slice->start_rec, o[3]->table, data.i);
	GetVarDataInit(var, o[3]->field->vardata, &vdata);

	n = vdata.count;
	end = o[0]->cooked.range.end;
	if (end < 0) end = n;

	for (i = o[0]->cooked.range.start ; i <= end ; i++) {
		if (i <= n)  {
			dvar = GetVarDataValue(var, o[3]->field->vardata, &vdata, i-1);
			rad = bbr(tesw[scan_len][i-1], target_temp);
			if (rad > 0) dvar.r = dvar.r / rad;
			else dvar.r = 0;
			PRINT_REAL(dvar);
		} else {
			printf("N/A");
		}
		if (i <= n) printf("%c", O_DELIM);
	}
}


/**
 **
 ** These routines aren't written very well.
 ** GetVarDataValue assumes all its input values are safe
 **
 **/

DATA
GetVarDataValue(PTR vraw, VARDATA *vardata, VDATA *v, int offset)
{
	DATA data;

	data = ConvertVarData(v->ptr + offset * vardata->size, vardata);
	if (vardata->type == Q15) {
		data.r = data.i * v->scale;
	}
	return(data);
}

VDATA *
GetVarDataInit(PTR vraw, VARDATA *vardata, VDATA *v)
{
	v->count = ConvertVaxVarByteCount(vraw, vardata) / vardata->size;

	switch (vardata->type) {
		case Q15:
			vraw += 2;	/* Skip the byte size */

			/* Assumption: Q15 exponent is always a signed integer */
			/* Read the Q15 (signed int) exponent from the first 2-bytes */

			v->scale = pow(2, (double)ConvertVarData(vraw, vardata).i - 15);

			/* Subtract 1 to exclude the quotient from the count of elements */
			v->count--;
			break;

		case VAX_VAR:
			v->scale = 1;
	}
	vraw += 2;	/* Skip the exponent */
	v->ptr = vraw;

	return(v);
}

/**
 ** black-body radiance support function
 ** 
 **/

#define C1 1.1909e-12
#define C2 1.43879

double 
bbr(double wn, double temp)
{
	/**
	 ** wn = wavenumber
	 ** temp = temperature in Kelvin
	 **/
    if (temp <= 0)
        return (0.0);
    return ((C1 * (wn * wn * wn)) / (exp(C2 * wn / temp) - 1.0));
}


#define START_CHAN_148  28
#define SCAN_STEP_148   5.29044097730104556043
#define wavenum_148(i,nw)       (floor(((START_CHAN_148 + (296/nw * i)) \
								 * SCAN_STEP_148)*10000.0)/ 10000.0)
#define channum_148(j,nw)       ((int)(((j/SCAN_STEP_148 - START_CHAN_148) / \
							  (296/nw)) + 0.5))


void
make_tesx()
{
	int i;

	if (tesw[1] == NULL) {
		tesw[1] = calloc(149, sizeof(double));
		tesw[2] = calloc(149*2, sizeof(double));

		for (i = 0 ; i < 148 ; i++) {
			tesw[1][i] = wavenum_148(i, 148);
			tesw[2][i] = wavenum_148(i, 296);
		}
	}
}