/* kernel_readin.c */

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#include <string.h>
#include "mex.h"
#include "defs.h"

#define MATLAB_KERNELS  (prhs[0])
#define GEOMETRY        (prhs[1])
#define DOSE            (plhs[0])
/* #define TERMA           (plhs[1])
#define KERMAC			(plhs[2])
#define DEFF            (plhs[3])
#define TERMA_MASK      (plhs[4])
#define DOSE_MASK       (plhs[5]) */

//function prototypes
void check_MATLAB_KERNELS(const mxArray *);
void check_GEOMETRY(const mxArray *);
void load_kernels(const mxArray *, MONO_KERNELS *);
void load_Geometry(const mxArray *, DOUBLE_GRID *, BEAM *);
int calc_deff(DOUBLE_GRID *,DOUBLE_GRID *, DOUBLE_GRID *, BEAM *);
int terma_kerma(DOUBLE_GRID *,DOUBLE_GRID *,DOUBLE_GRID *,MONO_KERNELS *,DOUBLE_GRID *);
int make_poly_kernel(MONO_KERNELS *, KERNEL *);
int calc_dose(DOUBLE_GRID *,DOUBLE_GRID *,DOUBLE_GRID *,KERNEL *,BEAM *, DOUBLE_GRID *);
int terma_dose_masks(DOUBLE_GRID *, DOUBLE_GRID *, BEAM *);

char errstr[200];  // error string that all routines have access to

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	int dims[3];  // dimensions of input arrays
	// matlab arrays for easy output if need be (just comment these lines out and uncomment
	// the #define lines above and the mxDestroy lines below).
	mxArray *TERMA, *KERMAC, *DEFF, *TERMA_MASK, *DOSE_MASK;

	char tmpstr[200];

	DOUBLE_GRID density;
	DOUBLE_GRID terma_mask;
	DOUBLE_GRID dose_mask;
	DOUBLE_GRID deff;
	DOUBLE_GRID terma;
	DOUBLE_GRID kermac;
	DOUBLE_GRID dose;

	MONO_KERNELS mono_kernels;
	KERNEL poly_kernel;

	BEAM beam;

	// check and load-in input data
	check_MATLAB_KERNELS(MATLAB_KERNELS); // see function comments for proper input format
	check_GEOMETRY(GEOMETRY);             // see function comments for proper input format

	load_kernels(MATLAB_KERNELS, &mono_kernels);
	load_Geometry(GEOMETRY,&density,&beam);

	// copy density grid geometry to all other grids on interest

	copy_grid_geometry(&density,&terma_mask);
	copy_grid_geometry(&density,&dose_mask);
	copy_grid_geometry(&density,&deff);
	copy_grid_geometry(&density,&terma);
	copy_grid_geometry(&density,&kermac);
	copy_grid_geometry(&density,&dose);

	// Allocate memory for all of the grids
	
	// dimensions of all the grids
	dims[0] = density.x_count;
	dims[1] = density.y_count;
	dims[2] = density.z_count;

	// Allocate memory for all of the data grids using Matlab mxArray functions. 
	// The idea behind this is for the outputs to send the outputs directly to Matlab.
	if ((TERMA_MASK = mxCreateNumericArray(3,dims,mxDOUBLE_CLASS,mxREAL))==NULL)
	  mexErrMsgTxt("Could not allocate memory for terma_mask grid"); 

	terma_mask.matrix = mxGetPr(TERMA_MASK);

	if ((DOSE_MASK = mxCreateNumericArray(3,dims,mxDOUBLE_CLASS,mxREAL))==NULL)
	  mexErrMsgTxt("Could not allocate memory for dose_mask grid"); 

	dose_mask.matrix = mxGetPr(DOSE_MASK);

	if ((DEFF = mxCreateNumericArray(3,dims,mxDOUBLE_CLASS,mxREAL))==NULL)
	  mexErrMsgTxt("Could not allocate memory for deff grid"); 

	deff.matrix = mxGetPr(DEFF);

	if ((TERMA = mxCreateNumericArray(3,dims,mxDOUBLE_CLASS,mxREAL))==NULL)
	  mexErrMsgTxt("Could not allocate memory for terma grid"); 

	terma.matrix = mxGetPr(TERMA);

	if ((KERMAC = mxCreateNumericArray(3,dims,mxDOUBLE_CLASS,mxREAL))==NULL)
	  mexErrMsgTxt("Could not allocate memory for kermac grid"); 

	kermac.matrix = mxGetPr(KERMAC);

	if ((DOSE = mxCreateNumericArray(3,dims,mxDOUBLE_CLASS,mxREAL))==NULL)
	  mexErrMsgTxt("Could not allocate memory for dose matrix"); 

	// Point dose matrix pointer at DOSE Matlab matrix data
	dose.matrix = mxGetPr(DOSE);

	/* start calculations */

	// create terma and dose masks
	if (SUCCESS != terma_dose_masks(&terma_mask,&dose_mask,&beam))
	{
		sprintf(tmpstr,"Failed in terma_dose_masks!\n");
		strcat(tmpstr,errstr);
		strcpy(errstr,tmpstr);
		mexErrMsgTxt(errstr);
	}

	//create polyenergetic kernel from mono kernels and fluence,mu data
	if (SUCCESS != make_poly_kernel(&mono_kernels,&poly_kernel) )
	{
		sprintf(tmpstr,"Failed making polyenergetic kernel!\n");
		strcat(tmpstr,errstr);
		strcpy(errstr,tmpstr);
		mexErrMsgTxt(errstr);
	}
	
	//create effective depth array from density array
	if (SUCCESS != calc_deff(&density,&deff,&terma_mask,&beam))
	{
		sprintf(tmpstr,"Failed in calc_deff!\n");
		strcat(tmpstr,errstr);
		strcpy(errstr,tmpstr);
		mexErrMsgTxt(errstr);
	}

	//create kerma and terma arrays
	//note kerma is collision kerma and is used for a kernel hardening correction
	if (SUCCESS != terma_kerma(&deff,&terma,&kermac,&mono_kernels,&terma_mask))
	{
		sprintf(tmpstr,"Failed in terma_kerma calculation!\n");
		strcat(tmpstr,errstr);
		strcpy(errstr,tmpstr);
		mexErrMsgTxt(errstr);
	}

	//use all this stuff to calculate dose
	if ( (SUCCESS != calc_dose(&density,&terma,&dose,&poly_kernel,&beam,&dose_mask)) ) 
	{
		sprintf(tmpstr,"Failed calculating dose!\n");
		strcat(tmpstr,errstr);
		strcpy(errstr,tmpstr);
		mexErrMsgTxt(errstr);
	}

	// Do not free any of the matrices because they are all MATLAB mxArrays that 
	// are used in the outputs.

	// destroy the arrays created with mxCreate
	mxDestroyArray(DEFF);
	mxDestroyArray(TERMA_MASK);
	mxDestroyArray(DOSE_MASK);
	mxDestroyArray(TERMA);
	mxDestroyArray(KERMAC);
}