/*-------------------------------------------------------------------------*/
/**
   @file	spectral_lines.c
   @author	N. Devillard
   @date	November 1999
   @version	$Revision: 1.14 $
   @brief	spectrum line handling routines
*/
/*--------------------------------------------------------------------------*/

/*
	$Id: spectral_lines.c,v 1.14 2001/11/12 13:41:29 yjung Exp $
	$Author: yjung $
	$Date: 2001/11/12 13:41:29 $
	$Revision: 1.14 $
*/

/*---------------------------------------------------------------------------
   								Includes
 ---------------------------------------------------------------------------*/

#include <math.h>

#include "spectral_lines.h"
#include "comm.h"
#include "emission_lines.h"

/*---------------------------------------------------------------------------
   								Defines
 ---------------------------------------------------------------------------*/

#ifndef ASCIILINESZ
#define ASCIILINESZ		1024
#endif

/*---------------------------------------------------------------------------
  							Function prototypes
 ---------------------------------------------------------------------------*/

static int emission_lines_compare(const void * e1, const void * e2)  ;

/*---------------------------------------------------------------------------
  							Function codes
 ---------------------------------------------------------------------------*/


/*-------------------------------------------------------------------------*/
/**
  @brief	Initialize a spectral line table.
  @param	path	Name of the table to initialize.
  @return	1 newly allocated spectral_table object.

  The received path string indicates where the spectral table should be
  loaded from.

  \begin{tabular}{ll}
  "oh"			&	OH table \\
  "Xe"   		&	Xenon table \\
  "Ar"   		&	Argon table \\
  "Xe+Ar"	    &	Xenon+Argon table \\
  /path/file	&	User-specified external table
  \end{tabular}

  If the given table name corresponds to a user-specified table, it is 
  loaded from the disk and dynamically allocated. 
  See spectral_table_parse_list() for the acceptable
  The returned table must be deallocated using spectral_table_destroy()
 */
/*--------------------------------------------------------------------------*/
spectral_table * spectral_table_init(char * path)
{
	spectral_table *	spt1 ;
	spectral_table *	spt2 ;
	spectral_table *	spt ;

	if (!strcmp(path, "oh")) {
		spt = spectral_table_select((spectral_table*)&emission_lines_table, 
									"oh") ;
	} else if (!strcmp(path, "Xe")) {
		spt = spectral_table_select((spectral_table*)&emission_lines_table, 
									"Xe") ;
	} else if (!strcmp(path, "Ar")) {
		spt = spectral_table_select((spectral_table*)&emission_lines_table, 
									"Ar") ;
	} else if (!strcmp(path, "Xe+Ar")) {
		spt1 = spectral_table_select((spectral_table*)&emission_lines_table, 
									"Xe") ;
		spt2 = spectral_table_select((spectral_table*)&emission_lines_table, 
									"Ar") ;
		spt = spectral_table_merge(spt1, spt2) ;
		spectral_table_destroy(spt1) ;
		spectral_table_destroy(spt2) ;
	} else {
		/* The specified list exists: load it */
		spt = spectral_table_parse_list(path);
		if (spt == NULL) {
			e_error("parsing file [%s]", path);
		}
	}
	spectral_table_sort(spt) ;
	return spt ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Create a spectral table
  @param	size	Size of the created table
  @return	the allocated spectral table
 */
/*--------------------------------------------------------------------------*/
spectral_table * spectral_table_create(int size)
{
	spectral_table	*	spt ;

	spt = malloc(sizeof(spectral_table)) ;
	spt->nlines = size ;
	spt->lines = malloc(size * sizeof(emission_line)) ;
	return spt ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Deallocate a spectral table.
  @param	spt	Spectral table to deallocate.
  @return	void
 */
/*--------------------------------------------------------------------------*/
void spectral_table_destroy(spectral_table * spt)
{
	if (spt == NULL) return ;
	free(spt->lines) ;
	free(spt) ;
	return ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Sort a spectral table 
  @param	table	spectral table to sort (modified) 
  @return	0 if Ok, -1 otherwise
 */
/*--------------------------------------------------------------------------*/
int spectral_table_sort(spectral_table * table)
{
	if (table == NULL) return -1 ;

	qsort((void*)table->lines, (size_t)table->nlines, sizeof(emission_line),
			emission_lines_compare) ;
	return 0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Compare two emission lines 
  @param	e1	first emission line
  @param	e2	second emission line
  @return	int 1 if e1->wavel < e2->wavel, -1 otherwise.
 */
/*--------------------------------------------------------------------------*/
static int emission_lines_compare(const void * e1, const void * e2) 
{
	if (((emission_line*)e1)->wavel < ((emission_line*)e2)->wavel)
		return -1 ;
    else 
		return 1 ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Dump a spectral table
  @param	table	spectral table to dump
  @param	out     Opened file pointer
  @return	0 if Ok, -1 otherwise
 */
/*--------------------------------------------------------------------------*/
int spectral_table_dump(
		spectral_table	*	table,
		FILE			*	out)
{
	int		i ;
		
	/* Check input parameters */
    if (table==NULL) return -1 ;
    if (out==NULL) out=stdout ;

	for (i=0 ; i<table->nlines ; i++) {
		fprintf(out, "%g\t%g\t%s\n", (table->lines[i]).wavel,
				(table->lines[i]).intens, (table->lines[i]).type) ;
	}
    return 0 ;  
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Merge two spectral tables
  @param	spt1	First spectral table
  @param	spt2	Second spectral table
  @return	a spectral table composed with the two input ones
 */
/*--------------------------------------------------------------------------*/
spectral_table * spectral_table_merge(
		spectral_table	*	spt1,
		spectral_table	*	spt2)
{
	spectral_table	*	spt ;
	int					i ;

	spt = spectral_table_create(spt1->nlines + spt2->nlines) ;
	
	/* Merge the tw input tables */
	if (spt != NULL) {
		for (i=0 ; i<spt1->nlines ; i++) {
			(spt->lines[i]).wavel = (spt1->lines[i]).wavel ;
			(spt->lines[i]).intens = (spt1->lines[i]).intens ;
			strcpy((spt->lines[i]).type, (spt1->lines[i]).type) ;
		}
		for (i=0 ; i<spt2->nlines ; i++) {
			(spt->lines[i+spt1->nlines]).wavel = (spt2->lines[i]).wavel ;
			(spt->lines[i+spt1->nlines]).intens = (spt2->lines[i]).intens ;
			strcpy((spt->lines[i+spt1->nlines]).type, (spt2->lines[i]).type) ;
		}
	}
	return spt ;
}

/*-------------------------------------------------------------------------*/
/**
  @brief	Select lines in a spectral table
  @param	ref		reference spectral table
  @param	type	type of lines selected
  @return	a spectral table created with selected lines
 */
/*--------------------------------------------------------------------------*/
spectral_table * spectral_table_select(
		spectral_table	*	ref,
		char			*	type)
{
	spectral_table	*	selected ;
	int					nb_selected ;
	int					i, j ;

	/* Initialize */
	nb_selected = 0 ;

	/* Count the number of selected lines */
	for (i=0 ; i<ref->nlines ; i++) {
		if (!strcmp(type, (ref->lines[i]).type)) nb_selected ++ ;
	}
	
	if (nb_selected == 0) return NULL ;
	
	/* Allocate the output table */
	selected = spectral_table_create(nb_selected) ;

	/* Copy the selected lines in the output spectral table */
	j = 0 ;
	for (i=0 ; i<ref->nlines ; i++) {
		if (!strcmp(type, (ref->lines[i]).type)) {
			(selected->lines[j]).wavel = (ref->lines[i]).wavel ;
			(selected->lines[j]).intens = (ref->lines[i]).intens ;
			strcpy((selected->lines[j]).type, (ref->lines[i]).type) ;
			j++ ;
		}
	}
	return selected ;
}

		
/*-------------------------------------------------------------------------*/
/**
  @brief	Read a spectral table from an external file.
  @param	path	Full path name of the file to load.
  @return	1 pointer to newly allocated spectral table.

  This function loads an external spectral table into a spectral_table
  object. The file must be an ASCII file, following these properties:

  \begin{itemize}
  \item Lines starting with a '#' are comments and ignored.
  \item Blank lines are ignored.
  \item Spectral lines are given as two values per line, separated by any
  number of blanks or tabs. The first value gives the wavelength in
  Angstroems, the second gives the relative intensity, which must be
  consistent with all lines in the same table.
  \end{itemize}

  The returned object must be deallocated using spectral_table_destroy().
 */
/*--------------------------------------------------------------------------*/
spectral_table * spectral_table_parse_list(char * path)
{
	spectral_table 	*	spt ;
	FILE	 		*	listfile ;
	int					i ;
	int					values_on_line ;
	int					nlines ;
	double				wave,
						irel ;
	int					lineno ;
	char				line[ASCIILINESZ];

	if ((listfile=fopen(path, "r"))==NULL) {
		e_error("cannot open file [%s]", path);
		return NULL ;
	}

	/* Count how many lines in the given file */
	nlines = 0 ;
	lineno = 0 ;
	while (fgets(line, ASCIILINESZ, listfile)!=NULL) {
		lineno++ ;
		if (line[0]!='#') {
			values_on_line = sscanf(line, "%lg %lg", &wave, &irel);
			if (values_on_line != 2) {
				e_error("in file %s (%d): expected two values",
						path, lineno);
				return NULL ;
			}
			nlines++ ;
		}
	}

	/* Allocate output table */
	spt = malloc(sizeof(spectral_table));
	spt->nlines = nlines ;
	spt->lines = malloc(nlines * sizeof(emission_line)) ;

	/* Now load line information */
	rewind(listfile);
	i=0 ;
	while (fgets(line, ASCIILINESZ, listfile)!=NULL) {
		if (line[0]!='#') {
			values_on_line = sscanf(line, "%lg %lg", &wave, &irel);
			if (values_on_line==2) {
				(spt->lines[i]).wavel = wave ;
				(spt->lines[i]).intens = irel ;
				strcpy((spt->lines[i]).type, "EF") ;
				i++ ;
			}
		}
	}
	fclose(listfile);
	return spt ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Build a 1d signal from a spectral table 
  @param	spt			Spectral table to use.
  @param	wave_min	Lower bound of the spectral range.
  @param	wave_max	Upper bound of the spectral range.
  @param	size		Size of the signal to generate.
  @return	1 pointer to a newly allocated array of 'size' doubles.

  Provide an allocated spectral table and a wavelength range in Angstroems
  as [wavemin, wavemax], and the size in pixels of the signal to generate.
  The returned array is a 1d signal (of doubles) containing the spectral
  lines as pixels.

  Returns NULL in case of error (e.g. invalid wavelength range).
 */
/*--------------------------------------------------------------------------*/
double * spectral_table_build_signal(
		spectral_table	*	spt,
		double				wave_min,
		double				wave_max,
		int					size)
{
	double	*	arr ;
	double	*	smooth ;
	int			i, j ;
	double		wdist ;
	double		wapp ;
	int			i_min, i_max ;

    /* Locate indexes of min and max wavelengths */
    i=0 ;
    while (((spt->lines[i]).wavel<wave_min) && (i<spt->nlines)) i++ ;
    i_min=i ;
    while (((spt->lines[i]).wavel<wave_max) && (i<spt->nlines)) i++ ;
    i_max=i ;

	if (i_min==spt->nlines) {
		e_error("wavelength range [%g %g] outside known bounds",
				wave_min, wave_max);
		return NULL ;
	}

    /* Build up a signal from the list of lines */
    arr = malloc(size * sizeof(double));
    wdist = (wave_max - wave_min) / (double)(size-1) ;
    for (i=0 ; i<size ; i++) {
        arr[i] = (double)0;
        wapp = wave_min + (double)i * wdist ;
        for (j=i_min ; j<i_max ; j++) {
            if (fabs((spt->lines[j]).wavel - wapp) < wdist) {
                arr[i] += (double)(spt->lines[j]).intens ;
            }
        }
    }
	/* Convolve with a gaussian to smooth out the signal */
	smooth = malloc(size*sizeof(double));
	smooth[0] = (2.0 * arr[0] + arr[1]) / 3.0 ;
	smooth[size-1] = (arr[size-2] + 2.0 * arr[size-1]) / 3.0 ;
	for (i=1 ; i<(size-1) ; i++) {
		smooth[i] = (arr[i-1] + 2.0 * arr[i] + arr[i+1]) * 0.25 ;
	}
	free(arr);
    return smooth ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Output a list of table lines as a signal to a file.
  @param	table_name		Name of the table to output.
  @param	outfilename		Name of the output file.
  @param	wave_min		Lower bound of the wavelength range.
  @param	wave_max		Upper bound of the wavelength range.
  @param	nsamples		Number of samples to produce.
  @return	void

  This function builds a 1d signal based on the requested table name and
  wavelength range, and outputs it to an ASCII file. This is useful to dump
  out lists of lines for e.g. debugging purposes, to compare a calibrated
  signal with a list of reference lines.

  Provide NULL or the character string "STDOUT" to output the list to
  stdout. Any other name specifies a file. If another file by the same name
  already exists, it will be overwritten.
 */
/*--------------------------------------------------------------------------*/
void spectral_table_build_spectrum(
		char	*	table_name,
		char	*	outfilename,
		double		wave_min,
		double		wave_max,
		int			nsamples)
{
	spectral_table	*	spt;
	FILE			*	out ;
	double			*	spectral_signal ;
	int					i ;
	double				wave_step ;

	spt = spectral_table_init(table_name);
	spectral_signal = spectral_table_build_signal(spt,
												  wave_min,
												  wave_max,
												  nsamples);
	spectral_table_destroy(spt);
	wave_step = (wave_max - wave_min) / (double)(nsamples-1);

	if (outfilename == NULL) {
		out = stdout ;
	} else if (!strcmp(outfilename, "STDOUT")) {
		out = stdout ;
	} else {
		out = fopen(outfilename, "w");
	}
	if (out==NULL) {
		e_error("cannot create file [%s]", outfilename);
		free(spectral_signal);
		return ;
	}
	for (i=0 ; i<nsamples ; i++) {
		fprintf(out, "%g\t%g\n",
				wave_min + (double)i * wave_step,
				spectral_signal[i]);
	}
	if (out!=stdout)
		fclose(out);
	free(spectral_signal);
	return ;
}