/***************************************************************************
 *		  Copyright (C) 1994  Charles P. Peterson                  *
 *	     4007 Enchanted Sun, San Antonio, Texas 78244-1254             *
 *              Email: Charles_P_Peterson@fcircus.sat.tx.us                *
 *                                                                         *
 *		  This is free software with NO WARRANTY.                  *
 *	      See gfft.c, or run program itself, for details.              *
 *		      Support is available for a fee.                      *
 ***************************************************************************
 *
 * Program:     gfft--General FFT analysis
 * File:        calibrat.c
 * Purpose:     calibration(s) 
 * Author:      Charles Peterson (CPP)
 * History:     15-March-1994 CPP; Created.
 *              5-Aug-94 CPP (1.05); add new calibrations to end of list
 * Comment:     Simple 'relative' calibration with linear interpolation.
 *              Any number of calibration(s) may be in effect.
 *              NOT complex deconvolution (maybe in a few years).
 *              Linear interpolation, and flat extension are used
 *              (Flat extension helps with fuzzy endpoints).
 */

#include <stdio.h>
#include <stdlib.h>  /* strtod */
#include <string.h>  /* strtok */
#include <math.h>    /* log10, pow */

#include "gfft.h"
#include "settings.h"
#include "wbench.h"

#define CAL_BUFSIZ 133
#define PROGRESS_INCREMENT 128

int calibration_list__count (struct cal_st **clistp)
{
    int count = 0;
    struct cal_st *clist = *clistp;

    while (clist)
    {
	count++;
	clist = clist->next;
    }
    return count;
}


void calibration_list__cancel (struct cal_st **clistp)
{
    struct cal_st *clist = *clistp;
/*
 * Free tail of list recursively, then free this node
 */
    if (clist)
    {
	if (clist->next)
	{
	    calibration_list__cancel (&clist->next);
	}
	gfree (clist->frequencies);
	gfree (clist->amplitudes);
/*	gfree (clist->filename); potentially dangerous */
	gfree (clist);
    }
    *clistp = NULL;
}


void calibration_list__add (struct cal_st **clistp,
			    char *cname, 
			    FILE *cfile,
			    BOOLEAN db_scale)
{
    int progress_count = 0;
    struct cal_st *clist = *clistp;
    char cal_buf[CAL_BUFSIZ];
    char *token;
    char *beyondpointer;
    double *freq_array = NULL;
    float *ampl_array = NULL;
    double frequency;
    float amplitude;
    long size = 0;
    struct cal_st *new_cal;
/*
 * First, read in file data
 */
    while (fgets (cal_buf, CAL_BUFSIZ, cfile))
    {
	if (progress_count++ % PROGRESS_INCREMENT == 0)
	{
	    progress_requester__update (-1);
	}
    /*
     * Skip over comment lines
     */
	if (cal_buf[0] == '#' || cal_buf[0] == ';' || cal_buf[0] == '!')
	{
	    continue;
	}
    /*
     * get frequency value
     */
	token = strtok (cal_buf, " ");
	if (!token)
	{
	    error_message (INVALID_CALIBRATION_FILE);
	    RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
	}
	frequency = strtod (token, &beyondpointer);
	if (*beyondpointer != '\0')
	{
	    error_message (INVALID_CALIBRATION_FILE);
	    RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
	}
	token = strtok (NULL, " ");
	if (!token)
	{
	    error_message (INVALID_CALIBRATION_FILE);
	    RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
	}
	amplitude = (float) strtod (token, &beyondpointer);
	if (*beyondpointer != '\0' && *beyondpointer != '\n')
	{
	    error_message (INVALID_CALIBRATION_FILE);
	    RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
	}
    /*
     * Add new values to arrays
     */
	size++;
	freq_array = grealloc (freq_array, (size * sizeof frequency), 
			       NOTHING_SPECIAL);
	ampl_array = grealloc (ampl_array, (size * sizeof amplitude), 
			       NOTHING_SPECIAL);
	freq_array[size-1] = frequency;
	ampl_array[size-1] = amplitude;
    }
    if (size <= 0)
    {
	error_message (INVALID_CALIBRATION_FILE);
	RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
    }

/*
 * create new node
 */
    new_cal = gmalloc (sizeof (struct cal_st), NOTHING_SPECIAL);
    new_cal->next = NULL;
    new_cal->filename = cname;
    new_cal->size = size;
    new_cal->db = db_scale;
    new_cal->frequencies = freq_array;
    new_cal->amplitudes = ampl_array;
    new_cal->index = 0;
/*
 * insert new node in list, or at start
 */
    if (clist)
    {
	while (clist->next)
	{
	    clist = clist->next;
	}
	clist->next = new_cal;
    }
    else
    {
	*clistp = new_cal;
    }
}

void calibration_list__reset (struct cal_st **clistp)
{
    if (*clistp)
    {
	(*clistp)->index = 0;
	calibration_list__reset (&((*clistp)->next));
    }
}

void calibration_list__write (FILE *fp, char *cstr, struct cal_st **clistp)
{
    struct cal_st *clist = *clistp;

    while (clist)
    {
	if (clist->db)
	{
	    fprintf (fp, "%sDBCalibrate %s\n", cstr, clist->filename);
	}
	else
	{
	    fprintf (fp, "%sCalibrate %s\n", cstr, clist->filename);
	}
	clist = clist->next;
    }
}


float calibration_list__apply (struct cal_st **clistp,
			       float value, double frequency)
{
    struct cal_st *clist = *clistp;
    long index;
    BOOLEAN outrange = TRUE;
    double factor;

/*
 * Assumption is made that frequencies are applied in an increasing
 * order.  When starting over, call calibration_list__reset
 *
 * Make local copy of 'index'
 */
    index = clist->index;
/*
 * Find applicable index
 */
    while (frequency > clist->frequencies[index])
    {
	outrange = FALSE;
	if (++index >= clist->size)
	{
	    index--;
	    outrange = TRUE;
	    break;
	}
    }
    clist->index = index;  /* update */
    
/*
 * Compute the factor to use
 */

    if (frequency == clist->frequencies[index] || outrange)
    {
	factor = clist->amplitudes[index];
    }
    else
    {
    /*
     * linear interpolation is fun
     * (Sorry I don't bother with anything more sophisticated yet)
     */
	double x1 = clist->frequencies[index-1];
	double x2 = clist->frequencies[index];
	double y1 = clist->amplitudes[index-1];
	double y2 = clist->amplitudes[index];
	double slope = (y2 - y1) / (x2 - x1);
	double offset = y1 - (slope * x1);

	factor = (slope * frequency) + offset;
    }

/*
 * Apply
 */
    if (Db)
    {
	if (clist->db)
	{
	    value -= factor;
	}
	else
	{
	    double dbfactor;
	    if (factor <= 0.0) RAISE_ERROR (NOTHING_SPECIAL);
	    dbfactor = 20.0 * log10 (factor);
	    value -= dbfactor;
	}
    }
    else
    {
	if (!clist->db)
	{
	    if (factor == 0.0) RAISE_ERROR (NOTHING_SPECIAL);
	    value /= factor;
	}
	else
	{
	    value /= pow (10.0, (factor * 0.05));
	}
    }
/*
 * Do next calibration in list
 */
    if (clist->next)
    {
	value = calibration_list__apply (&clist->next, value, frequency);
    }

    return value;
}