
//${SPLAT_DIR}/splatsh $0 ${1+"$@"}; exit;
// Note first line is blank, do not remove it as this starts "sh", which
// runs the next line, which runs splatsh on this file. This header
// section is ignored as it is a beanshell comment, that "sh" never sees.

/**
 * Name:
 *    fitgauss2
 *
 * Purpose:
 *    Measure a series of spectral-lines by fitting gaussians.
 *
 * Usage:
 *    fitgauss2 positions_file width results_file spectrum
 *
 * Description:
 *    Fit a background subtracted spectrum with a series of gaussians
 *    using the specified expected positions and line widths to
 *    constrain the fits.
 *
 * Arguments
 *    positions_file
 *       A simple text file that contains the expected wavelengths of
 *       the lines. Each wavelength should be on a separate line.
 *    width
 *       A width over which each line should fitted. This is
 *       centered on the wavelength in the positions table.
 *    results_file
 *       A file for the results of fitting each line.
 *    spectrum
 *       The spectrum containing the data to be fitted. This should be
 *       background subtracted such that the lines have positive data
 *       values.
 *
 * Language:
 *    Beanshell (Java-based scripting language).
 *
 * @since $Date: 2003-06-04 10:20:35 +0100 (Wed, 04 Jun 2003) $
 * @since 26-NOV-2001
 * @author Peter W. Draper
 * @version $Id: fitgauss2 1030 2003-06-04 09:20:35Z pwd $
 * @copyright Copyright (C) 2001 Central Laboratory of the Research Councils */

// Import any classes we're using from SPLAT.
import uk.ac.starlink.splat.data.SpecData;
import uk.ac.starlink.splat.data.SpecDataFactory;
import uk.ac.starlink.splat.util.GaussianFitter;
import uk.ac.starlink.splat.util.QuickLineFitter;

// Method to print the usage message and exit.
usage()
{
    print( "Usage: fitgauss2 positions_file width results_file spectrum" );
    exit();
}

// Basic check of input args.
if ( bsh.args.length != 4 ) {
    usage();
}

// Gather args, doing some simple type checking.

// Positions file. Must exist already.
positionsFile = pathToFile( bsh.args[0] );
if ( ! positionsFile.exists() || ! positionsFile.canRead() ) {
    print( "!! Cannot read positions file: " + bsh.args[0] );
    usage();
}

// Check width is a floating point value.
try {
    fitWidth = Double.parseDouble( bsh.args[1] );
}
catch (Exception badNumber) {
    print( badNumber.printStackTrace() );
    usage();
}

// Results file.
resultsFile = pathToFile( bsh.args[2] );

// Access the spectrum. This is opened by the SpecDataFactory
// class, which understands the various specifications and data
// types (i.e. you can use any specification that SPLAT understands).
try {
    spectrum = SpecDataFactory.getReference().get( bsh.args[3] );
}
catch (Exception badSpecification) {
    print( badSpecification.printStackTrace() );
    print( "!! Failed to access spectrum: " + bsh.args[3] );
    usage();
}

// Open the input and results files. These are wrapped so we can write
// and read them a line at a time.
try {
    inputStream = new BufferedReader( new FileReader( positionsFile ) );
}
catch (Exception fileNotOpened) {
    print( fileNotOpened.printStackTrace() );
    print( "!! Failed to open input file" );
    exit();
}

try {
    outputStream = new PrintStream( new FileOutputStream( resultsFile ) );
}
catch (Exception fileNotOpened) {
    print( fileNotOpened.printStackTrace() );
    print( "!! Failed to open output file" );
    exit();
}

// Method to print a line to standard output and file.
record ( line ) {
    print( line );
    outputStream.println( line );
}

//  Add header section to output file.
record( "# Gaussian fit results:" );
record( "# Line ID, Centre, Peak, Sigma, Flux, Rms" );

//  Write one lines worth of results to the terminal and the results
//  file.
//  ** Modify this and the header section if you want a different format. **
recordLine( line, gaussFitter, gaussRms ) {
    record( line + ", " + 
            gaussFitter.getCentre() + ", " +
            gaussFitter.getScale() + ", " +
            gaussFitter.getSigma() + ", " +
            gaussFitter.getFlux() + ", " +
            gaussRms );
    record( "" );
}

// Now read the input file a line at a time.
nLines = 0;
nFit = 0;
while ( ( line = inputStream.readLine() ) != null ) {
    nLines++;

    // Should be a double precision number. If fails complain and pass
    // on to the next.
    try {
        guess = Double.parseDouble( line );
    }
    catch (Exception numberConversion) {
        print( numberConversion.printStackTrace() );
        continue;
    }

    //  Get a sub-spectrum that just encompasses the line data within
    //  the width.
    range = new double[2];
    range[0] = guess - fitWidth;
    range[1] = guess + fitWidth;
    subSpectrum = spectrum.getSect( "Line: " + line, range );

    //  No data, then complain and pass on.
    if ( subSpectrum == null ) {
        print( "" );
        print( "!! Failed to find any data for line: " + line );
        print( "" );
        continue;
    }

    //  Access the spectral raw data values.
    specCoords = subSpectrum.getXData();
    specData = subSpectrum.getYData();
    specDataErrors = subSpectrum.getYDataErrors();

    //  Fit the line! First by the QuickFitter to get an estimate of
    //  the peak value.
    quickFitter = new QuickLineFitter( specCoords, specData, null );
    if ( quickFitter.isAbsorption() ) {
        quickPeak = -quickFitter.getPeak();
    }
    else {
        quickPeak = quickFitter.getPeak();
    }

    //  Finally by the Gaussian.
    if ( specDataErrors == null ) {
        gaussFitter = new GaussianFitter( specCoords, specData,
                                          quickPeak,
                                          quickFitter.getCentre(),
                                          0.5 * quickFitter.getWidth() );

    }
    else {
        gaussFitter = new GaussianFitter( specCoords, specData,
                                          specDataErrors,
                                          quickFitter.getPeak(),
                                          quickFitter.getCentre(),
                                          0.5 * quickFitter.getWidth() );
    }

    //  Record the results.
    recordLine( line, gaussFitter, 
                gaussFitter.calcRms( specCoords, specData ) );
    nFit++;
}

print( "Number of lines fitted: " + nFit + " (out of " + nLines + ")" );

// Exit, closing all files.
inputStream.close();
outputStream.close();
exit();
