#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include "fake.h"
#include "ff_irtm.h"

/* Exported */
void ff_irtm(OSTRUCT **o);
void cook_print_ff_irtm(OSTRUCT * o, int scaled, int frame_size);
void print_header_ff_irtm(OSTRUCT *o);
FIELD *make_ff_irtm(LIST *tables);
char *nameof_ff_irtm = "irtm";

/* Local */
static double ff_irtm_i(
	int     irtm_ch,     /* IRTM channel number */
	int     scan_len,    /* TES scan length */
	int     st_offset,   /* offset into TES response array */
	int     n,           /* Number of calrad[] spectra */
	double  calrad[],    /* Actual calibrated radiance values */
	int    *status       /* Success/failure flag */
);
static void ff_irtm_fill_tes_wavelengths();
static double interpolate_tes_response(double irtm_resp[],
                                       double irtm_wavelength[],
                                       int n, double tes_wavelength,
                                       int *status);

#define PRINT_ARRAY(array, st, ed) {\
	int l;\
\
	for(l = st; l <= ed; l++){\
		fprintf(stderr, "%g ", (array)[l]);\
		if ((l-st+1)%4 == 0){ fprintf(stderr, "\n"); }\
	}\
}

#define TRACE(x) 

#define SCAN_LENGTH_MULTIPLIER 143
#define MAX_SCAN_LEN 2

/* Number of TES channels for the given scan-length, where
	scan-length equals 1 or 2 */
#define N_TES_CHAN(scan_length) ((scan_length)*SCAN_LENGTH_MULTIPLIER)

/* Brightness temperatures calculated from TES spectra */
static double btemps[N_TES_CHAN(MAX_SCAN_LEN)];

/* Calculate wave number given a wave length */
#define WAVE_NO(wave_length) (1e4/(wave_length))

/* A low non-zero value to be used instead of zero calrads */
#define LOW_NON_ZERO_VALUE 1e-30



FIELD *
make_ff_irtm(LIST *tables)
{
	int        i;
	FIELD     *f;
	FAKEFIELD *ff;
	char      *dep[2] = {
		"obs.scan_len",
		"rad.cal_rad"
	};

	assert((f = (FIELD *)calloc(sizeof(FIELD), 1)) != NULL);
	f->name = strdup(nameof_ff_irtm);
	f->iformat = REAL;
	f->size = 4;
	f->dimension = N_IRTM_CHAN;

	assert((ff = (FAKEFIELD *)calloc(sizeof(FAKEFIELD), 1)) != NULL);
	ff->print_header_function = print_header_ff_irtm;
	ff->cook_function = cook_print_ff_irtm;
	ff->fptr = ff_irtm;

	/* Add in dependencies */
	ff->nfields = 2;
	assert((ff->fields = (FIELD **)calloc(sizeof(FIELD *), ff->nfields)) != NULL);

	for (i = 0; i < ff->nfields; i++){
		ff->fields[i] = FindField(dep[i], tables);
		if (ff->fields[i] == NULL){
			fprintf(stderr,
				"Dependent field %s not found while processing fake field %s.\n",
				dep[i], f->name);
			return NULL;
		}
	}

	f->fakefield = ff;

	return f;
}

void
cook_print_ff_irtm(OSTRUCT *o, int scaled, int frame_size)
{
	o->cooked.frame_size = frame_size;
	o->cooked.print_func = ff_irtm;
	o->cooked.flags = OF_FAKEFIELD;

	o->cooked.range.start = MAX(1, o->range.start);
	o->cooked.range.end = MIN(o->field->dimension, o->range.end);
	if(o->cooked.range.end < 0){
		o->cooked.range.end = o->field->dimension;
	}

	if (o->cooked.range.start > o->cooked.range.end){
		fprintf(stderr, "Fake field %s has switched ranges [%d,%d].\n",
				o->text, o->cooked.range.start, o->cooked.range.end);
		abort();
	}
}

void
print_header_ff_irtm(OSTRUCT *o)
{
	int i;

	for(i = o->cooked.range.start; i <= o->cooked.range.end; i++){
		if (o->cooked.frame_size == TEXT_OUTPUT_FRAME_SIZE){
			fprintf(ofptr, "%s[%d]", o->text, i);
			if (i < o->cooked.range.end){
				print_field_delim();
			}
		}
		else {
			print_output_mask_real();
		}
	}
}


/*
** o[0] points to fake field named "irtm";
** o[1] points to rad.scan_len;
** o[2] points to rad.calrad vardata field;
*/
void
ff_irtm(OSTRUCT **o)
{
	double   calrad[N_TES_CHAN(MAX_SCAN_LEN)];
	int      i, j, k, n;
	int      st_idx, ed_idx; /* TES spectrum index range for current IRTM ch */
	int      st_offset; 
	int      scan_len;
	int      scan_len_0; /* scan_len -1 */
	DATA     data;
	VDATA    vdata;
	char    *var;
	char    *raw;
	int      status;
	double   irtm;

	/* variables for caching only */
	int      step;
	VARDATA *v = o[2]->field->vardata;
	int      t = v->type;

	/*
	** Get scan_len
	*/
	raw = o[1]->slice->start_rec + o[1]->field->start;
	data = ConvertData(raw, o[1]->field);
	scan_len = *data.str == '1' ? 1 : 2;
	scan_len_0 = scan_len - 1;
	TRACE(fprintf(stderr, "\nff_irtm.c: scan_len=%d\n", scan_len))

	/*
	** Get byte pointer to calrad vardata
	*/
	raw = o[2]->slice->start_rec + o[2]->field->start;
	data = ConvertData(raw, o[2]->field);
	TRACE(fprintf(stderr, "ff_irtm.c: cal_rad[] byte offset=%#x\n", data.i))

	if (data.i > -1) {
		/* if we have a good byte offset, use it to get the calibrated
			radiance */
		var = GiveMeVarPtr(o[2]->slice->start_rec, o[2]->table, data.i);
		GetVarDataInit(var, o[2]->field->vardata, &vdata);
		n = vdata.count;
	}
	else {
		n = 0; /* Assume #(cal_rad[]) == 0 */
	}
	TRACE(fprintf(stderr, "ff_irtm.c: #(cal_rad[])==%d\n", n))

	/* Process each of the IRTM channels one by one. */
	for (i = o[0]->cooked.range.start; i <= o[0]->cooked.range.end; i++){

		/* Note: i & both the start & end indices are 1-based */
		st_idx = tes_resp[i-1].tes_ch[scan_len_0].start;
		ed_idx = tes_resp[i-1].tes_ch[scan_len_0].end;
		step = tes_resp_wt_step[scan_len_0];

		/* If we are working with scan-len=1 then the start index into
			the TES-response array may be one greater than that of
			scan-len=2 */
		st_offset = st_idx - tes_resp[i-1].tes_ch[MAX_SCAN_LEN-1].start;

		if ((ed_idx/step) > n){
			/* Available calrad values are insufficient */
			irtm = 0.0;
		}
		else {
			/* Extract the required portion of calibrated radiance */
			for (j=st_idx, k=0; j<=ed_idx; j+=step, k++){
				data = GetVarDataValue(var, v, &vdata, (j-1)/step);
				/* Note k+st_idx-1 = (j-1)/step */

				/* Store calibrated radiance compactly. */
				calrad[k] = ((t == Q15) ? data.r : (double)data.i);
				/* Note: k = (j-st_idx)/step */
			}

			/* k = 1+(ed_idx-st_idx+1)/step */
			irtm = ff_irtm_i(i, scan_len, st_offset, 
									(ed_idx-st_idx+1)/step, calrad, &status);
		}
		print_real(irtm, o[0]->cooked.frame_size);

		/* Print the field delimiter if we are processing IRTM value
			other than last IRTM requested */
		if ((i < o[0]->cooked.range.end) &&
			(o[0]->cooked.frame_size == TEXT_OUTPUT_FRAME_SIZE)){
			print_field_delim();
		}
	}
}

#if 0
	static double H = 6.626e-34; /* Plank's constant J-s */
	static double C = 2.998E8;   /* Speed of light m/s */
	static double K = 1.381e-23; /* Boltzmann's contstant J/K */
#endif


static double
ff_irtm_i(
	int     irtm_ch,     /* IRTM channel number */
	int     scan_len,    /* TES scan length */
	int     st_offset,   /* offset into TES response array */
	int     n,           /* Number of calrad[] spectra */
	double  calrad[],    /* Actual calibrated radiance values */
	int    *status       /* Success/failure flag */
)
{
	int     irtm_ch_0 = irtm_ch - 1;   /* irtm_ch -1 */
	int     scan_len_0 = scan_len - 1; /* scan_len -1 */
	double  irtm = 0.0;
	int     i, j;
	int     step;

	/* Skip every other tes_resp_wt[] index while processing calrad 
		data with scan-len=1 */
	step = tes_resp_wt_step[scan_len_0];
	/* n = (ed_idx - st_idx + 1)/step; */

	/* Take care of potentially bad data in calibrated radiance.
		Note that calrad[] is already compacted, i.e. all required
		elements of calrad[] are consecutive even when scan-len=2 */
	for(i = 0; i < n; i++){
		if (calrad[i] == 0){
			/* ln(0) = (kaboom!) so replace zeroes with low non-zero values */
			calrad[i] = LOW_NON_ZERO_VALUE;
		}
		else if (calrad[i] < 0){
			/* ln(<0) = (kaboom!) assume bad data */
			*status = 0; /* IRTM value returned is no good */
			return 0.0;
		}
	}

	/* Compute spectral brightness temperatures from the TES spectra */
	for(i=0, j=st_offset; i<n; i++, j+=step){
		btemps[i] = 
		 	btemp(WAVE_NO(tes_resp[irtm_ch_0].wl_wt_pairs[j].wl), calrad[i]);

			#if 0
			/* H*C*100.0*WAVE_NO(tes_resp[irtm_ch_0].wl_wt_pairs[j].wl)/
				(K*log(2*H*pow(C*100.0, 2.0)*
				pow(WAVE_NO(tes_resp[irtm_ch_0].wl_wt_pairs[j].wl), 3.0)/
				calrad[i] + 1.0)); */
			#endif
	}

	irtm = 0.0;
	for(i=0, j=st_offset; i<n; i++, j+=step){
		/* Normalize TES (response) weights and integrate over the TES spectra */
		irtm += btemps[i] *
				tes_resp[irtm_ch_0].wl_wt_pairs[j].wt / 
				tes_resp[irtm_ch_0].sum_tes_resp_wt[scan_len_0];
	}

	*status = 1; /* IRTM value is good */
	return irtm;
}