/* massXpert - the true massist's program.
   --------------------------------------
   Copyright(C) 2006,2007 Filippo Rusconi

   http://www.filomace.org/massXpert

   This file is part of the massXpert project.

   The massxpert project is the successor to the "GNU polyxmass"
   project that is an official GNU project package(see
   www.gnu.org). The massXpert project is not endorsed by the GNU
   project, although it is released ---in its entirety--- under the
   GNU General Public License. A huge part of the code in massXpert
   is actually a C++ rewrite of code in GNU polyxmass. As such
   massXpert was started at the Centre National de la Recherche
   Scientifique(FRANCE), that granted me the formal authorization to
   publish it under this Free Software License.

   This software is free software; you can redistribute it and/or
   modify it under the terms of the GNU  General Public
   License version 3, as published by the Free Software Foundation.
   

   This software is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this software; if not, write to the

   Free Software Foundation, Inc.,

   51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.

  
   
   NOTE : This work is based on a previous work by Dirk Nolting,
   licensed under the GNU General Public License:
   
   ***************************************************************************
   *  file                 :  ipc.c                                          *
   *  copyright            :(C) 2001-2005 by Dirk Nolting                   *
   *  email                : nolting@uni-duesseldorf.de                      *
   ***************************************************************************

   ***************************************************************************
   *                                                                         *
   *   This program is free software; you can redistribute it and/or modify  *
   *   it under the terms of the GNU General Public License as published by  *
   *   the Free Software Foundation; either version~2 of the License, or     *
   *  (at your option) any later version.                                   *
   *                                                                         *
   ***************************************************************************
   *
   */


/////////////////////// Qt includes
#include <QtGui>


/////////////////////// Std includes
#include <math.h>
#include <algorithm>
#include <limits> // for std::numeric_limits

using namespace std;

/////////////////////// Local includes
#include "isotopicPatternCalculationDlg.hpp"
#include "isotopicPeak.hpp"
#include "application.hpp"


namespace massXpert
{

  IsotopicPatternCalculationDlg::IsotopicPatternCalculationDlg 
 (QWidget *parent, const QList<Atom *> &atomList)
    : QDialog(parent), m_atomList(atomList)
  {
    Q_ASSERT(parent);

    m_aborted = false;
    m_filePath = "";
  
    m_ui.setupUi(this);

    // Set the name of the polymer chemistry definition.

    const PolChemDef & polChemDef = 
      static_cast<CalculatorWnd *>(parent)->polChemDef();
    m_ui.polChemDefLabel->setText(polChemDef.name());
  


    m_ui.minimumProbabilityDoubleSpinBox->setRange(0, 1000000);
    m_ui.minimumProbabilityDoubleSpinBox->setValue(0.0001);
    m_minimumProbability = 0.0001;
  
    m_ui.maximumPeaksSpinBox->setRange(0, 1000000);
    m_ui.maximumPeaksSpinBox->setValue(1000);
    m_maximumPeaks = 1000;
  
    m_ui.chargeSpinBox->setRange(0, 1000000);
    m_ui.chargeSpinBox->setValue(1);
    m_charge = 1;
  
    m_ui.resolutionSpinBox->setRange(0, 1000000);
    m_ui.resolutionSpinBox->setValue(15000);
    m_resolution = 15000;
  
    m_ui.progressBar->setRange(0, 100);
    m_ui.progressBar->setValue(0);

    QSettings settings 
     (static_cast<Application *>(qApp)->configSettingsFilePath(), 
       QSettings::IniFormat);
  
    settings.beginGroup("isotopic_pattern_calculation_dlg");

    restoreGeometry(settings.value("geometry").toByteArray());

    settings.endGroup();

    connect(m_ui.executePushButton,
	     SIGNAL(clicked()),
	     this,
	     SLOT(execute()));
  
    connect(m_ui.abortPushButton,
	     SIGNAL(clicked()),
	     this,
	     SLOT(abort()));

    connect(m_ui.outputFilePushButton,
	     SIGNAL(clicked()),
	     this,
	     SLOT(outputFile()));

  }


  void 
  IsotopicPatternCalculationDlg::closeEvent(QCloseEvent *event)
  {
    if (event)
      printf("%s", "");
  
    QSettings settings 
     (static_cast<Application *>(qApp)->configSettingsFilePath(), 
       QSettings::IniFormat);
  
    settings.beginGroup("isotopic_pattern_calculation_dlg");
  
    settings.setValue("geometry", saveGeometry());
  
    settings.endGroup();
  }


  IsotopicPatternCalculationDlg::~IsotopicPatternCalculationDlg()
  {
  
  
  }


  bool
  IsotopicPatternCalculationDlg::fetchValidateInputData()
  {
    Application *application = static_cast<Application *>(qApp);

    QString formula = m_ui.formulaLineEdit->text();
  
    if (!formula.isEmpty())
      m_formula.setFormula(formula);
    else
      {
	QMessageBox::warning(0,
			      tr("massXpert: Isotopic Pattern Calculation"),
			      tr("Enter a valid formula."),
			      QMessageBox::Ok);
	return false;
      }
  
    // We want to validate the formula and in the mean time construct
    // the list of all the AtomCount objects(first true), and since
    // the formula is reused we also ensure that that list is reset
    //(second true).
  
    if (!m_formula.validate(m_atomList, true, true))
      {
	QMessageBox::warning(0,
			      tr("massXpert: Isotopic Pattern Calculation"),
			      tr("Formula(%1) is not valid.")
			      .arg(formula),
			      QMessageBox::Ok);
	return false;
      }
  
    // Shall we use localization for all the numerical output?
    bool locale = m_ui.localeCheckBox->checkState() == Qt::Checked ?
      true : false;
  

    int totAtoms = m_formula.totalAtoms();
    m_ui.progressBar->setRange(0, totAtoms);
  
    int totIsotopes = m_formula.totalIsotopes(m_atomList);

    m_mono = 0;
    m_avg = 0;
  
    if (!m_formula.accountMasses(m_atomList, &m_mono, &m_avg, 1))
      return false;
  
    QString textEditText;

    if (locale)
      textEditText += tr("INPUT\n=====\n\nformula: %1\n"
			  "Mono Mass: %2 \t Avg mass: %3\n"
			  "Total number of atoms: %4\n"
			  "Total number of isotopes: %5")
	.arg(formula)
	.arg(application->locale().
	      toString(m_mono, 'f', MXP_OLIGOMER_DEC_PLACES))
	.arg(application->locale().
	      toString(m_avg, 'f', MXP_OLIGOMER_DEC_PLACES))
	.arg(totAtoms)
	.arg(totIsotopes);
    else
      textEditText += tr("INPUT\n=====\n\nformula: %1\n"
			  "Mono mass: %2 \t Avg mass: %3\n"
			  "Total number of 1) atoms: %4 ; 2) isotopes: %5\n\n")
	.arg(formula)
	.arg(QString().setNum(m_mono, 'f', MXP_OLIGOMER_DEC_PLACES))
	.arg(QString().setNum(m_avg, 'f', MXP_OLIGOMER_DEC_PLACES))
	.arg(totAtoms)
	.arg(totIsotopes);

  
    // The list of atomCount objects has to be made with fully detailed
    // isotopes...
    if (!m_formula.deepAtomCopy(m_atomList))
      {
	QMessageBox::warning(0,
			      tr("massXpert: Isotopic Pattern Calculation"),
			      tr("Failed to deep-copy atom list."),
			      QMessageBox::Ok);

	textEditText += tr("Failed to deep-copy atom list.\n");
      
	// Out put some text in the textEdit:
	m_ui.resultTextEdit->append(textEditText);
      
	return false;
      }

    m_charge = m_ui.chargeSpinBox->value();
    m_resolution = m_ui.resolutionSpinBox->value();
    m_maximumPeaks = m_ui.maximumPeaksSpinBox->value();
 
    textEditText += tr("Charge: %1 ; Resolution: %2 ; Maximum peaks: %3 ;")
      .arg(m_charge)
      .arg(m_resolution)
      .arg(m_maximumPeaks);
    
    m_minimumProbability = m_ui.minimumProbabilityDoubleSpinBox->value();
    textEditText += tr(" Minimum probability: ");
    if (locale)
      textEditText += application->locale().
	toString(m_minimumProbability, 'f', MXP_OLIGOMER_DEC_PLACES);
    else
      textEditText += QString().
	setNum(m_minimumProbability, 'f', MXP_OLIGOMER_DEC_PLACES);;
  
    textEditText += 
      tr("\n\nGaussian simulation of a single peak is performed by "
	  "computing 25 points on each side of the isotopic mass "
	  "according to the following exponential:\n"
	  "f(x) = relativeIntensity * exp(( -(x - mono)^2 ) / c^2 )\n\n"
	  "Full Width At Half Maximum = "
	  "(mono / resolution) = 2 * sqrt(ln(2) * c \n");
  
    double fullWidthHalfMax = m_mono / m_resolution;
    textEditText += tr("FWHM = ");
    if (locale)
      textEditText += application->locale().
	toString(fullWidthHalfMax, 'f', MXP_OLIGOMER_DEC_PLACES);
    else
      textEditText += QString().
	setNum(fullWidthHalfMax, 'f', MXP_OLIGOMER_DEC_PLACES);;
  
    textEditText += "\n";

    double c = fullWidthHalfMax /(2 * sqrt(log(2)));
    double c2 = pow(c, 2);

    textEditText += tr("with c = ");
    if (locale)
      textEditText += application->locale().
	toString(c, 'f', MXP_OLIGOMER_DEC_PLACES);
    else
      textEditText += QString().setNum(c, 'f', MXP_OLIGOMER_DEC_PLACES);;
  
    textEditText += tr(" and c^2 = ");
    if (locale)
      textEditText += application->locale().
	toString(c2, 'f', MXP_OLIGOMER_DEC_PLACES);
    else
      textEditText += QString().setNum(c2, 'f', MXP_OLIGOMER_DEC_PLACES);;
  
    textEditText += "\n============================================\n\n";
  
    // Out put some text in the textEdit:
    m_ui.resultTextEdit->append(textEditText);
    
    return true;
  }


  void 
  IsotopicPatternCalculationDlg::execute()
  {
    Application *application = static_cast<Application *>(qApp);

    bool inserted = false;
    int wholeCount = 0;

    // Clear the textEdit widget so that we can start out-putting text
    // into it.
    m_ui.resultTextEdit->clear();

    if (!fetchValidateInputData())
      return;

    double fullWidthHalfMax = m_mono / m_resolution;


    QList<IsotopicPeak *> oldPeakList;
  
    IsotopicPeak *firstPeak = new IsotopicPeak();
    firstPeak->setMass(0);
    firstPeak->setProbability(1);
  
    oldPeakList.append(firstPeak);
	  
    // After the validation of the formula, it contains a list of
    // allocated AtomCount objects. Get a handy shorthand to it.

    const QList<AtomCount *> &atomCountList = m_formula.atomCountList();
  
    QTime time;
    time.start();
  
    m_ui.feedbackLineEdit->setText(tr("Computing the peaks in the pattern."));  

    for (int iter = 0; iter < atomCountList.size(); ++iter)
      {
	AtomCount *atomCount = atomCountList.at(iter);
	int isotopeCount = atomCount->isotopeList().size();
      
	//       qDebug() << "Current atom:" << atomCount->symbol()
	// 		<< "with" << isotopeCount << "isotope(s)";
      
	for(int jter = 0; jter < atomCount->count(); ++jter)
	  {
	    // 	  qDebug() << "Atom count:"<< atomCount->count()
	    // 		    << "current run:" << jter + 1
	    // 		    << "wholeCount:" << wholeCount;	  
	  
	    m_ui.progressBar->setValue(++wholeCount);
	    qApp->processEvents();
	  
	    if (m_aborted)
	      break;
	  	  
	    QList<IsotopicPeak *> newPeakList;

	    for (int kter = 0; kter < oldPeakList.size(); ++kter)
	      {
		IsotopicPeak *peakListPeak = oldPeakList.at(kter);
	      
		for(int lter = 0; lter < isotopeCount; ++lter)
		  {
		    Isotope *isotope = atomCount->isotopeList().at(lter);
		  
		    // 		  qDebug() << "Current isotope:" << isotope->mass();
		  
		    IsotopicPeak *newPeak = 
		      new IsotopicPeak(peakListPeak->mass() +
					isotope->mass(),
					0,
					peakListPeak->probability() *
					(isotope->abundance() / 100),
					0);
		  
		    // 		  qDebug() << "Newly created peak:"
		    // 			    << newPeak->mass() << "--"
		    // 			    << newPeak->probability();
		  
		    if (!newPeakList.size())
		      {
			newPeakList.append(newPeak);

			// 		      qDebug() << "Appended peak:"
			// 			    << newPeak->mass() << "--"
			// 				<< newPeak->probability();
		      }
		    else
		      {
			for(int pter = 0; pter < newPeakList.size(); ++pter)
			  {
			    IsotopicPeak *iteratedPeak = 
			      newPeakList.at(pter);
			  
			    if (newPeak->mass() > 
				iteratedPeak->mass())
			      continue;
			  
			    newPeakList.insert(pter, newPeak);
	
			    // 			  qDebug() << "Inserted peak:"
			    // 				    << newPeak->mass() << "--"
			    // 				    << newPeak->probability();
			  
			    inserted = true;

			    break;
			  }
		      
			// If we did not insert the peak, then that
			// means that we arrived here because we
			// finished iterating in the list... Thus append
			// the peak.
			if(!inserted)
			  {
			    newPeakList.append(newPeak);

			    // 			  qDebug() << "Appended peak:"
			    // 				    << newPeak->mass() << "--"
			    // 				    << newPeak->probability();	  
			  }
			else
			  {
			    // Reinitialize the bool to false.
			    inserted = false;
			  }
		      
			// 		      qDebug() << "List of peaks:";
		      
			// 		      for (int zter = 0; zter < newPeakList.size(); ++zter)
			// 			qDebug() << newPeakList.at(zter)->mass()
			// 				  << "--" 
			// 				  << newPeakList.at(zter)->probability();
		      }
		  }
	      }
	  
	    qDeleteAll(oldPeakList);
	    oldPeakList.clear();
	  
	    oldPeakList = newPeakList;
	    newPeakList.clear();
	    
	    // We now have to remove peaks that are too near one from
	    // the other.
	  
	    for (int nter = 0; nter < oldPeakList.size() - 1; ++nter)
	      {
		double absDiff = fabs(oldPeakList.at(nter)->mass() -
				       oldPeakList.at(nter + 1)->mass());
	      
		if(absDiff <= fullWidthHalfMax)
		  {
		    // 		  qDebug() << "absDiff is" << absDiff;
		  
		    // Remove the forward peak that is too close to
		    // current one, but do not forget to account for
		    // that one from the abundance point of vue.
		  
		    oldPeakList.at(nter)->
		      incrementProbability(oldPeakList.at(nter + 1)->
					    probability());

		    // 		  qDebug() << "Removed peak"
		    // 			    << oldPeakList.at(nter + 1)->mass() << "--"
		    // 			    << oldPeakList.at(nter + 1)->probability();
		  
		    delete oldPeakList.takeAt(nter + 1);
		    --nter;
		  } 
	      }

	    if (m_maximumPeaks >= 2 && oldPeakList.size() > m_maximumPeaks)
	      {
		// The user wants at most m_maximumPeaks peaks, which means
		// we have to remove all the ones in excess.
	      
		for(int pter = m_maximumPeaks; 
		     pter < oldPeakList.size(); ++pter)
		  {
		    delete oldPeakList.takeAt(pter);
		  }
	      }
	  }
            
	if(m_aborted)
	  break;
      }
    m_aborted = false;
  
    m_ui.progressBar->setRange(0, 100);
    m_ui.progressBar->setValue(0);

    qApp->processEvents();

    // Take into account the charge of the analyte ! Take advantage of
    // this run to compute the sum of all the probabilities, that we'll
    // need later. Also take advantage of this looping to get the value
    // of the most probable peak. We'll need it to compute the relative
    // intensities later.
 
    m_ui.feedbackLineEdit->setText(tr("Accounting for the charge."));
 
    double sumProbabilities = 0;
    double greatestProbability = 0;
    
    for (int iter = 0; iter < oldPeakList.size(); ++iter)
      {
	IsotopicPeak *isotopicPeak = oldPeakList.at(iter);
      
	isotopicPeak->setMass(isotopicPeak->mass() / m_charge);
      
	sumProbabilities += isotopicPeak->probability();
	if(isotopicPeak->probability() > greatestProbability)
	  greatestProbability = isotopicPeak->probability();
      }
  
    // If we have actually divided the mass by m_charge > 1, then the
    // peaks need to be reprocessed so as to remove all the ones too
    // close to each other.
  
    if (m_charge > 1)
      {
	for(int nter = 0; nter < oldPeakList.size() - 1; ++nter)
	  {
	    if (fabs(oldPeakList.at(nter)->mass() -
		      oldPeakList.at(nter + 1)->mass()) <= fullWidthHalfMax)
	      {
		// Remove the forward peak that is too close to
		// current one, but do not forget to account for
		// that one from the abundance point of vue.
	      
		oldPeakList.at(nter)->
		  incrementProbability(oldPeakList.at(nter + 1)->
					probability());
	      
		// qDebug() << "Removed peak"
		// << oldPeakList.at(nter + 1)->mass() << "--"
		// << oldPeakList.at(nter + 1)->probability();
	      
		delete oldPeakList.takeAt(nter + 1);
	      } 
	  }
      }
    
    m_ui.feedbackLineEdit->setText(tr("Computing the relative intensity."));

    // Compute the relative intensity.
  
    for (int iter = 0; iter < oldPeakList.size(); ++iter)
      {
	IsotopicPeak *isotopicPeak = oldPeakList.at(iter);
      
	double relativeIntensity = 
	 ( isotopicPeak->probability() / greatestProbability) * 100;
      
	if(relativeIntensity > std::numeric_limits<double>::min())
	  isotopicPeak->setRelativeIntensity(relativeIntensity);
      }
  
    int elapsedTime = time.elapsed() /(1000 * 60);

    // And now we should do something with the peaks !!!

    m_ui.feedbackLineEdit->setText(tr("Formatting the results"));

    bool locale = m_ui.localeCheckBox->checkState() == Qt::Checked ?
      true : false;
    
    double c = fullWidthHalfMax /(2 * sqrt(log(2)));
    double c2 = pow(c, 2);
    double increment = fullWidthHalfMax / 50;

    QString results(tr("RESULTS for this computation(duration: %1 min)\n"
			 "====================================\n\n").
		     arg(elapsedTime));
  
    for (int oter = 0; oter < oldPeakList.size(); ++oter)
      {
	double mass = oldPeakList.at(oter)->mass();
	double prob = oldPeakList.at(oter)->probability();
	double relInt = oldPeakList.at(oter)->relativeIntensity();
      
	if(prob > m_minimumProbability)
	  {
	    results += tr("Mass: ");
	    if (locale)
	      results += application->locale().
		toString(mass, 'f', MXP_OLIGOMER_DEC_PLACES);
	    else
	      results += QString().setNum(mass, 'f', MXP_OLIGOMER_DEC_PLACES);
	    results += tr(" -- Probability: ");
	    if (locale)
	      results += application->locale().
		toString(prob, 'f', MXP_OLIGOMER_DEC_PLACES);
	    else
	      results += QString().setNum(prob, 'f', MXP_OLIGOMER_DEC_PLACES);
	    results += tr(" -- Rel. intensity: ");
	    if (locale)
	      results += application->locale().
		toString(relInt, 'f', MXP_OLIGOMER_DEC_PLACES);
	    else
	      results += QString().
		setNum(relInt, 'f', MXP_OLIGOMER_DEC_PLACES);
	    results += "\n";
	  }
      }
    
    m_ui.resultTextEdit->append(results);

    //   qDebug() << "fullWidthHalfMax" << fullWidthHalfMax
    // 	    << "x-increment" << increment
    // 	    << "c" << c
    // 	    << "c^2" << c2;

    // For each isotopic peak, compute a gaussian.
    // Write the results to the file, peak after peak.

    QFile file(m_filePath);
    QTextStream stream;
  
    bool ret = file.open(QIODevice::WriteOnly |
			  QIODevice::Truncate | 
			  QIODevice::Text);
    if (ret)
      {
	stream.setDevice(&file);
	stream.setCodec("UTF-8");
      }
    else
      {
	QMessageBox::warning(0, 
			      tr("massXpert - Isotopic Pattern Calculation"),
			      tr("Failed to export "
				  "the data: file could not be opened."),
			      QMessageBox::Ok);
      
	stream.setStatus(QTextStream::ReadCorruptData);
      }
  
    for (int iter = 0; iter < oldPeakList.size(); ++iter)
      {
	double prob = oldPeakList.at(iter)->probability();
      
	// If the formula has a big number of atoms, then the first few
	// peaks will be of almost negligible probability, and we'll
	// want to skip them.
	if(prob < m_minimumProbability)
	  continue;
      
	double mono = oldPeakList.at(iter)->mass();
	double relInt = oldPeakList.at(iter)->relativeIntensity();
      
	double leftPoint = mono -(2 * fullWidthHalfMax);
	double rightPoint = mono +(2 * fullWidthHalfMax);
      
	QString outputString;
	stream << outputString;
	outputString.clear();
      
	//       outputString += "#########";
	//       if (locale)
	// 	outputString += application->locale().toString(mono, 'f', 5);
	//       else
	// 	outputString += QString().setNum(mono, 'f', 5);
	//       outputString += "--> [";
	//       if (locale)
	// 	outputString += application->locale().toString(leftPoint, 'f', 5);
	//       else
	// 	outputString += QString().setNum(leftPoint, 'f', 5);
	//       outputString += "--";
	//       if (locale)
	// 	outputString += application->locale().toString(rightPoint, 'f', 5);
	//       else
	// 	outputString += QString().setNum(rightPoint, 'f', 5);
	//       outputString += "]\n";
      
	for(double x = leftPoint; x <= rightPoint; x += increment)
	  {
	    double y = relInt * exp(( - pow((x - mono), 2) ) / c2);
	  
	    if (locale)
	      outputString += application->locale().
		toString(x, 'f', MXP_OLIGOMER_DEC_PLACES);
	    else
	      outputString += QString().
		setNum(x, 'f', MXP_OLIGOMER_DEC_PLACES);
	    outputString += "  ";
	    if (locale)
	      outputString += application->locale().
		toString(y, 'f', MXP_OLIGOMER_DEC_PLACES);
	    else
	      outputString += QString().
		setNum(y, 'f', MXP_OLIGOMER_DEC_PLACES);
	    outputString += "\n";
	  }

	if(stream.status() == QTextStream::ReadCorruptData)
	  m_ui.resultTextEdit->append(outputString);
	else
	  stream << outputString;       
      
	//       qDebug() << outputString;
	outputString.clear();
      }

    m_ui.feedbackLineEdit->setText(tr(""));

    // Finished sending all the peak data to the file. QFile will close
    // itself.
  }


  void 
  IsotopicPatternCalculationDlg::abort()
  {
    m_aborted = true;
  }


  void 
  IsotopicPatternCalculationDlg::outputFile()
  {
    QString name;
  
    m_filePath = 
      QFileDialog::getSaveFileName(this, tr("Export Raw Text File"),
				    QDir::homePath(),
				    tr("Any file type(*)"));
  }


#if 0
  /***************************************************************************
   *  file                 :  ipc.c                                          *
   *  copyright            :(C) 2001-2005 by Dirk Nolting                   *
   *  email                : nolting@uni-duesseldorf.de                      *
   ***************************************************************************/

  /***************************************************************************
   *                                                                         *
   *   This program is free software; you can redistribute it and/or modify  *
   *   it under the terms of the GNU General Public License as published by  *
   *   the Free Software Foundation; either version~2 of the License, or     *
   *  (at your option) any later version.                                   *
   *                                                                         *
   ***************************************************************************/

#include "global.h"
#include "ipc.h"
#include "pars.h"
#include "element.h"
#include "gp_out.h"
#include <time.h>
#include <signal.h>
#include <math.h>

  //#define SUMMARIZE_LEVEL 0.00000001
  //#define SUMMARIZE_LEVEL 0.0001
  //#define SUMMARIZE_LEVEL 0.001
#define SUMMARIZE_LEVEL 0.01

  compound *verbindung=NULL;
  isotope *m_peakList;
  int fast_calc=0;

  void free_list(isotope *target)
  {
    while(target->next)
      {
	target=target->next;
	free(target->previous);
      }
    free(target);
  }

  void cut_peaks(isotope *spectrum)
  {
    int dummy=1;

    while((spectrum->next) &&(dummy<fast_calc) )
      {
	++dummy;
	spectrum=spectrum->next;
      }

    if(spectrum->next)
      {
	free_list(spectrum->next);
	spectrum->next=NULL;
      }
  }

  void summarize_peaks()
  {  isotope *dummy,*d2;

    for(dummy=m_peakList;dummy;dummy=dummy->next)
      /* Differenz wegen Rundungsfehlern */
      while( dummy->next &&(dummy->next->mass - dummy->mass < SUMMARIZE_LEVEL) )
	{
	  d2=dummy->next;
	  dummy->next=d2->next;
	  if(dummy->next)
	    dummy->next->previous=dummy;
	  dummy->p+=d2->p;
	  free(d2);
	}
  }

  isotope *add_peak(isotope *base,isotope *peak)
  {
    static isotope *IsotopeIterator;

    if(!(base->mass))
      {
	peak->next=NULL;
	peak->previous=NULL;
	IsotopeIterator = peak;
	return peak;
      }

    if( peak->mass >= IsotopeIterator->mass )
      while((IsotopeIterator->next) &&(IsotopeIterator->mass < peak->mass))
	IsotopeIterator=IsotopeIterator->next;
    else
      {
	while((IsotopeIterator->previous) &&(IsotopeIterator->mass > peak->mass) )
	  IsotopeIterator=IsotopeIterator->previous;
	IsotopeIterator=IsotopeIterator->next;
      }

    if((IsotopeIterator->mass) >=(peak->mass) )
      {
	peak->next=IsotopeIterator;
	peak->previous=IsotopeIterator->previous;
	peak->previous->next=peak;
	IsotopeIterator->previous=peak;
	return base;
      }
    else
      {
	IsotopeIterator->next=peak;
	peak->next=NULL;
	peak->previous=IsotopeIterator;
	return base;
      }
    return 0;
  }

  int calculate_peaks(){
    compound *c;
    isotope *newPeakList,*p,*i,*np1;
    int amount;

    if(!(m_peakList=malloc(sizeof(isotope))))
      return 0;
    m_peakList->mass=0;
    m_peakList->p=1;
    m_peakList->previous=NULL;
    m_peakList->next=NULL;

    for (c = verbindung; c; c = c->next)
      {
	for(amount = 0; amount < c->amount; ++amount)
	  {
	    if (!(newPeakList=malloc(sizeof(isotope))))
	      return 0;
	  
	    newPeakList->mass = 0;
	  
	    for (p = m_peakList; p; p = p->next)
	      {
		for(i = c->isotopes; i; i = i->next)
		  {
		    //printf("working on isotope of mass %f\n", i->mass);

		    if (!(np1 = malloc(sizeof(isotope))))
		      return 0;
		  
		    np1->mass=p->mass + i->mass;
		    np1->p=p->p * i->p;
		  
		    if(!(newPeakList = add_peak(newPeakList,np1)))
		      return 0;
		  }
	      }
	  
	    free_list(m_peakList);
	    m_peakList = newPeakList;
	    summarize_peaks();
	  
	    if (fast_calc)
	      cut_peaks(m_peakList);
	  }
      }
  
    return 1;
  }


  void print_result(int digits,int charge){
    isotope *d;
    double maxp=0,relint=0,sump=0;
    int permutationen=0;

    printf("\n");
 
    for(d=m_peakList;d;d=d->next)
      {
	++permutationen;
	sump+=d->p;
	d->mass=d->mass / charge;
	d->mass=(rint( d->mass * pow(10,digits) ) / pow(10,digits) );
      }

    summarize_peaks();
    for(d=m_peakList;d;d=d->next)
      if(d->p > maxp)
	maxp=d->p;

    for(d=m_peakList;d;d=d->next)
      {
	if(( relint=(d->p/maxp)*100) > MIN_INT )
	  printf("M= %f, p= %e, rel. Int.= %f%%\n",
		 d->mass,d->p,relint);
      }
    if(!(fast_calc))
      printf("\nNumber of  permutations: %i\n",permutationen);
    else
      {
	sump=(rint(sump*10000)/100); 
	printf("\nCovered Intensity: %2.2f%% \n",sump);
      }
  }

  int main(int argc,char **argv){
    long seconds;
    int d=1,zeig_summenformel=0,calc_peaks=1,gnuplot=0,use_digits=USE_DIGITS,charge=1;
    char *gnuplotfile=NULL;
  
    if(!argv[d])
      {
	usage();
	return 1;
      }

#ifdef HAVE_SIGNAL
    signal(SIGHUP,SIG_IGN);
#endif

    if(!init_elements()){
      printf("Error in init_elements\n");
      return 1;
    }
    while(argv[d])
      {
	if(!strcmp(argv[d],"-f"))
	  {
	    ++d;
	    if(!argv[d])
	      {
		printf("Missing argument for -f\n");
		return 1;
	      }
	    fast_calc=strtol(argv[d],NULL,10);
	  }

	else if(!strcmp(argv[d],"-z"))
	  {
	    ++d;
	    if(!argv[d])
	      {
		printf("Missing argument for -z\n");
		return 1;
	      }
	    charge=strtol(argv[d],NULL,10);
	  }

	else if(!strcmp(argv[d],"-d"))
	  {
	    ++d;
	    if(!argv[d])
	      {
		printf("Missing argument for -d\n");
		return 1;
	      }
	    use_digits=strtol(argv[d],NULL,10);
	  }

	else if(!strcmp(argv[d],"-c")){
	  ++d;
	  if(!argv[d])
	    {
	      printf("Missing argument for -c\n");
	      return 1;
	    }
	  if(!pars_chem_form(argv[d])){
	    printf("Parser error.\n");
	    return 1;
	  }
	}

	else if(!strcmp(argv[d],"-g")){
	  ++d;
	  if(!argv[d]){
	    printf("Missing argument for -g\n");
	    return 1;
	  }
	  gnuplot=1;
	  if(!(gnuplotfile=strdup(argv[d]))){
	    printf("Not enough memory\n");
	    return 1;
	  }
	}


	else if(!strcmp(argv[d],"-p")){
	  ++d;
	  if(!argv[d]){
	    printf("Missing argument for -p\n");
	    return 1;
	  }
	  if(!(pars_peptid(argv[d]))){
	    printf("Error in peptid parser.\n");
	    return 1;
	  }
	}

	else if(!strcmp(argv[d],"-a")){
	  ++d;
	  if(!argv[d]){
	    printf("Missing argument for -a\n");
	    return 1;
	  }
	  if(!pars_amino_acid(argv[d])){
	    printf("Error in pars_amino_acid.\n");
	    return 1;
	  }
	}

	else if(!strcmp(argv[d],"-s"))
	  zeig_summenformel=1;

	else if(!strcmp(argv[d],"-h"))
	  {
	    usage();
	    return 1;
	  }

	else if(!strcmp(argv[d],"-x"))
	  calc_peaks=0;

	else
	  {
	    printf("Unknown flag: %s\n",argv[d]);
	    usage();
	    return 1;
	  }
	++d;
      }

    if(zeig_summenformel)
      {
	if(!print_sum())
	  {
	    printf("error while showing chemical formula.\n");
	    return 1;
	  }
      }
#ifdef HAVE_TIME
    seconds=time(0);
#endif

    if(calc_peaks)
      if(!calculate_peaks()){
	printf("Error in calculate_peaks\n");
	return 1;
      }

    print_result(use_digits,charge);

#ifdef HAVE_TIME
    seconds=time(0)-seconds;
    printf("Computing time: %li seconds.\n",seconds);
#endif

    if(gnuplot)
      if(!(make_gnuplot_output(gnuplotfile))){
	printf("Error while generating gnuplot file\n");
	return 1;
      }

    return 0;
  }


  void usage()
  {
    printf("\nThis is IPC v %s\nCopyright Dirk Nolting 2001-2005\n\n",VERSION);
    printf("\nSynopsis:\n ipc -d <int> -z <int> -f <int> <-a <amino acid> -c <chemical formula> -p <File> -g <name> -s -x -h\n\n");

    printf("-c <chemical formula> calculates the isotopic pattern of the \n");
    printf("   given chemical formula. Additive with -p\n");

    printf("-p <File> reads peptide sequenz(one letter notation) from the \n");
    printf("   given file and calculates the isotopic pattern. Additive with -c\n");

    printf("-a <amino acids> calculate isotopic pattern from given amino acids\n");
    printf("   in one letter notation\n");

    printf("-s gives the chemical formula. Usefull for -a or -p\n");
    printf("-g creates gnuplot output and stores it in the file <name> and <name>.gnu\n");

    printf("-x no calculation of the isotopic pattern. Usefull for -s\n");
    printf("-f fast calc, calculates only first <int> peaks\n");
    printf("-d <int> digits significant\n");
    printf("-z assume <int> charges on ion \n");
    printf("-h show this text\n\n");
  }
#endif

} // namespace massXpert