/***************************************************************************
 * traceutils.c:
 *
 * Generic routines to handle Traces.
 *
 * Written by Chad Trabant, IRIS Data Management Center
 *
 * modified: 2012.105
 ***************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "libmseed.h"

static int mst_groupsort_cmp ( MSTrace *mst1, MSTrace *mst2, flag quality );


/***************************************************************************
 * mst_init:
 *
 * Initialize and return a MSTrace struct, allocating memory if needed.
 * If the specified MSTrace includes data samples they will be freed.
 *
 * Returns a pointer to a MSTrace struct on success or NULL on error.
 ***************************************************************************/
MSTrace *
mst_init ( MSTrace *mst )
{
  /* Free datasamples, prvtptr and stream state if present */
  if ( mst )
    {
      if ( mst->datasamples )
	free (mst->datasamples);

      if ( mst->prvtptr )
	free (mst->prvtptr);

      if ( mst->ststate )
        free (mst->ststate);
    }
  else
    {
      mst = (MSTrace *) malloc (sizeof(MSTrace));
    }

  if ( mst == NULL )
    {
      ms_log (2, "mst_init(): Cannot allocate memory\n");
      return NULL;
    }

  memset (mst, 0, sizeof (MSTrace));

  return mst;
} /* End of mst_init() */


/***************************************************************************
 * mst_free:
 *
 * Free all memory associated with a MSTrace struct and set the pointer
 * to 0.
 ***************************************************************************/
void
mst_free ( MSTrace **ppmst )
{
  if ( ppmst && *ppmst )
    {
      /* Free datasamples if present */
      if ( (*ppmst)->datasamples )
        free ((*ppmst)->datasamples);

      /* Free private memory if present */
      if ( (*ppmst)->prvtptr )
        free ((*ppmst)->prvtptr);

      /* Free stream processing state if present */
      if ( (*ppmst)->ststate )
        free ((*ppmst)->ststate);

      free (*ppmst);

      *ppmst = 0;
    }
} /* End of mst_free() */


/***************************************************************************
 * mst_initgroup:
 *
 * Initialize and return a MSTraceGroup struct, allocating memory if
 * needed.  If the supplied MSTraceGroup is not NULL any associated
 * memory it will be freed.
 *
 * Returns a pointer to a MSTraceGroup struct on success or NULL on error.
 ***************************************************************************/
MSTraceGroup *
mst_initgroup ( MSTraceGroup *mstg )
{
  MSTrace *mst = 0;
  MSTrace *next = 0;

  if ( mstg )
    {
      mst = mstg->traces;

      while ( mst )
	{
	  next = mst->next;
	  mst_free (&mst);
	  mst = next;
	}
    }
  else
    {
      mstg = (MSTraceGroup *) malloc (sizeof(MSTraceGroup));
    }

  if ( mstg == NULL )
    {
      ms_log (2, "mst_initgroup(): Cannot allocate memory\n");
      return NULL;
    }

  memset (mstg, 0, sizeof (MSTraceGroup));

  return mstg;
} /* End of mst_initgroup() */


/***************************************************************************
 * mst_freegroup:
 *
 * Free all memory associated with a MSTraceGroup struct and set the
 * pointer to 0.
 ***************************************************************************/
void
mst_freegroup ( MSTraceGroup **ppmstg )
{
  MSTrace *mst = 0;
  MSTrace *next = 0;

  if ( *ppmstg )
    {
      mst = (*ppmstg)->traces;

      while ( mst )
	{
	  next = mst->next;
	  mst_free (&mst);
	  mst = next;
	}

      free (*ppmstg);

      *ppmstg = 0;
    }
} /* End of mst_freegroup() */


/***************************************************************************
 * mst_findmatch:
 *
 * Traverse the MSTrace chain starting at 'startmst' until a MSTrace
 * is found that matches the given name identifiers.  If the dataquality
 * byte is not 0 it must also match.
 *
 * Return a pointer a matching MSTrace otherwise 0 if no match found.
 ***************************************************************************/
MSTrace *
mst_findmatch ( MSTrace *startmst, char dataquality,
		char *network, char *station, char *location, char *channel )
{
  int idx;

  if ( ! startmst )
    return 0;

  while ( startmst )
    {
      if ( dataquality && dataquality != startmst->dataquality )
	{
	  startmst = startmst->next;
	  continue;
	}

      /* Compare network */
      idx = 0;
      while ( network[idx] == startmst->network[idx] )
	{
	  if ( network[idx] == '\0' )
	    break;
	  idx++;
	}
      if ( network[idx] != '\0' || startmst->network[idx] != '\0' )
	{
	  startmst = startmst->next;
	  continue;
	}
      /* Compare station */
      idx = 0;
      while ( station[idx] == startmst->station[idx] )
	{
	  if ( station[idx] == '\0' )
	    break;
	  idx++;
	}
      if ( station[idx] != '\0' || startmst->station[idx] != '\0' )
	{
	  startmst = startmst->next;
	  continue;
	}
      /* Compare location */
      idx = 0;
      while ( location[idx] == startmst->location[idx] )
	{
	  if ( location[idx] == '\0' )
	    break;
	  idx++;
	}
      if ( location[idx] != '\0' || startmst->location[idx] != '\0' )
	{
	  startmst = startmst->next;
	  continue;
	}
      /* Compare channel */
      idx = 0;
      while ( channel[idx] == startmst->channel[idx] )
	{
	  if ( channel[idx] == '\0' )
	    break;
	  idx++;
	}
      if ( channel[idx] != '\0' || startmst->channel[idx] != '\0' )
	{
	  startmst = startmst->next;
	  continue;
	}

      /* A match was found if we made it this far */
      break;
    }

  return startmst;
} /* End of mst_findmatch() */


/***************************************************************************
 * mst_findadjacent:
 *
 * Find a MSTrace in a MSTraceGroup matching the given name
 * identifiers, samplerate and is adjacent with a time span.  If the
 * dataquality byte is not 0 it must also match.
 *
 * The time tolerance and sample rate tolerance are used to determine
 * if traces abut.  If timetol is -1.0 the default tolerance of 1/2
 * the sample period will be used.  If samprratetol is -1.0 the
 * default tolerance check of abs(1-sr1/sr2) < 0.0001 is used (defined
 * in libmseed.h).  If timetol or sampratetol is -2.0 the respective
 * tolerance check will not be performed.
 *
 * The 'whence' flag will be set, when a matching MSTrace is found, to
 * indicate where the indicated time span is adjacent to the MSTrace
 * using the following values:
 * 1: time span fits at the end of the MSTrace
 * 2: time span fits at the beginning of the MSTrace
 *
 * Return a pointer a matching MSTrace and set the 'whence' flag
 * otherwise 0 if no match found.
 ***************************************************************************/
MSTrace *
mst_findadjacent ( MSTraceGroup *mstg, flag *whence, char dataquality,
		   char *network, char *station, char *location, char *channel,
		   double samprate, double sampratetol,
		   hptime_t starttime, hptime_t endtime, double timetol )
{
  MSTrace *mst = 0;
  hptime_t pregap;
  hptime_t postgap;
  hptime_t hpdelta;
  hptime_t hptimetol = 0;
  hptime_t nhptimetol = 0;
  int idx;

  if ( ! mstg )
    return 0;

  *whence = 0;

  /* Calculate high-precision sample period */
  hpdelta = (hptime_t)(( samprate ) ? (HPTMODULUS / samprate) : 0.0);

  /* Calculate high-precision time tolerance */
  if ( timetol == -1.0 )
    hptimetol = (hptime_t) (0.5 * hpdelta);   /* Default time tolerance is 1/2 sample period */
  else if ( timetol >= 0.0 )
    hptimetol = (hptime_t) (timetol * HPTMODULUS);

  nhptimetol = ( hptimetol ) ? -hptimetol : 0;

  mst = mstg->traces;

  while ( mst )
    {
      /* post/pregap are negative when the record overlaps the trace
       * segment and positive when there is a time gap. */
      postgap = starttime - mst->endtime - hpdelta;

      pregap = mst->starttime - endtime - hpdelta;

      /* If not checking the time tolerance decide if beginning or end is a better fit */
      if ( timetol == -2.0 )
	{
	  if ( ms_dabs((double)postgap) < ms_dabs((double)pregap) )
	    *whence = 1;
	  else
	    *whence = 2;
	}
      else
	{
	  if ( postgap <= hptimetol && postgap >= nhptimetol )
	    {
              /* Span fits right at the end of the trace */
	      *whence = 1;
	    }
	  else if ( pregap <= hptimetol && pregap >= nhptimetol )
	    {
              /* Span fits right at the beginning of the trace */
	      *whence = 2;
	    }
          else
            {
	      /* Span does not fit with this Trace */
	      mst = mst->next;
	      continue;
            }
	}

      /* Perform samprate tolerance check if requested */
      if ( sampratetol != -2.0 )
	{
	  /* Perform default samprate tolerance check if requested */
	  if ( sampratetol == -1.0 )
	    {
	      if ( ! MS_ISRATETOLERABLE (samprate, mst->samprate) )
		{
		  mst = mst->next;
		  continue;
		}
	    }
	  /* Otherwise check against the specified sample rate tolerance */
	  else if ( ms_dabs(samprate - mst->samprate) > sampratetol )
	    {
	      mst = mst->next;
	      continue;
	    }
	}

       /* Compare data qualities */
       if ( dataquality && dataquality != mst->dataquality )
	{
	  mst = mst->next;
	  continue;
	}

      /* Compare network */
      idx = 0;
      while ( network[idx] == mst->network[idx] )
	{
	  if ( network[idx] == '\0' )
	    break;
	  idx++;
	}
      if ( network[idx] != '\0' || mst->network[idx] != '\0' )
	{
	  mst = mst->next;
	  continue;
	}
      /* Compare station */
      idx = 0;
      while ( station[idx] == mst->station[idx] )
	{
	  if ( station[idx] == '\0' )
	    break;
	  idx++;
	}
      if ( station[idx] != '\0' || mst->station[idx] != '\0' )
	{
	  mst = mst->next;
	  continue;
	}
      /* Compare location */
      idx = 0;
      while ( location[idx] == mst->location[idx] )
	{
	  if ( location[idx] == '\0' )
	    break;
	  idx++;
	}
      if ( location[idx] != '\0' || mst->location[idx] != '\0' )
	{
	  mst = mst->next;
	  continue;
	}
      /* Compare channel */
      idx = 0;
      while ( channel[idx] == mst->channel[idx] )
	{
	  if ( channel[idx] == '\0' )
	    break;
	  idx++;
	}
      if ( channel[idx] != '\0' || mst->channel[idx] != '\0' )
	{
	  mst = mst->next;
	  continue;
	}

      /* A match was found if we made it this far */
      break;
    }

  return mst;
} /* End of mst_findadjacent() */


/***************************************************************************
 * mst_addmsr:
 *
 * Add MSRecord time coverage to a MSTrace.  The start or end time will
 * be updated and samples will be copied if they exist.  No checking
 * is done to verify that the record matches the trace in any way.
 *
 * If whence is 1 the coverage will be added at the end of the trace,
 * whereas if whence is 2 the coverage will be added at the beginning
 * of the trace.
 *
 * Return 0 on success and -1 on error.
 ***************************************************************************/
int
mst_addmsr ( MSTrace *mst, MSRecord *msr, flag whence )
{
  int samplesize = 0;

  if ( ! mst || ! msr )
    return -1;

  /* Reallocate data sample buffer if samples are present */
  if ( msr->datasamples && msr->numsamples >= 0 )
    {
      /* Check that the entire record was decompressed */
      if ( msr->samplecnt != msr->numsamples )
	{
	  ms_log (2, "mst_addmsr(): Sample counts do not match, record not fully decompressed?\n");
	  ms_log (2, "  The sample buffer will likely contain a discontinuity.\n");
	}

      if ( (samplesize = ms_samplesize(msr->sampletype)) == 0 )
	{
	  ms_log (2, "mst_addmsr(): Unrecognized sample type: '%c'\n",
		  msr->sampletype);
	  return -1;
	}

      if ( msr->sampletype != mst->sampletype )
	{
	  ms_log (2, "mst_addmsr(): Mismatched sample type, '%c' and '%c'\n",
		  msr->sampletype, mst->sampletype);
	  return -1;
	}

      mst->datasamples = realloc (mst->datasamples,
				  (size_t) (mst->numsamples * samplesize + msr->numsamples * samplesize) );

      if ( mst->datasamples == NULL )
	{
	  ms_log (2, "mst_addmsr(): Cannot allocate memory\n");
	  return -1;
	}
    }

  /* Add samples at end of trace */
  if ( whence == 1 )
    {
      if ( msr->datasamples && msr->numsamples >= 0 )
	{
	  memcpy ((char *)mst->datasamples + (mst->numsamples * samplesize),
		  msr->datasamples,
		  (size_t) (msr->numsamples * samplesize));

	  mst->numsamples += msr->numsamples;
	}

      mst->endtime = msr_endtime (msr);

      if ( mst->endtime == HPTERROR )
	{
	  ms_log (2, "mst_addmsr(): Error calculating record end time\n");
	  return -1;
	}
    }

  /* Add samples at the beginning of trace */
  else if ( whence == 2 )
    {
      if ( msr->datasamples && msr->numsamples >= 0 )
	{
	  /* Move any samples to end of buffer */
	  if ( mst->numsamples > 0 )
	    {
	      memmove ((char *)mst->datasamples + (msr->numsamples * samplesize),
		       mst->datasamples,
		       (size_t) (mst->numsamples * samplesize));
	    }

	  memcpy (mst->datasamples,
		  msr->datasamples,
		  (size_t) (msr->numsamples * samplesize));

	  mst->numsamples += msr->numsamples;
	}

      mst->starttime = msr->starttime;
    }

  /* If two different data qualities reset the MSTrace.dataquality to 0 */
  if ( mst->dataquality && msr->dataquality && mst->dataquality != msr->dataquality )
    mst->dataquality = 0;

  /* Update MSTrace sample count */
  mst->samplecnt += msr->samplecnt;

  return 0;
} /* End of mst_addmsr() */


/***************************************************************************
 * mst_addspan:
 *
 * Add a time span to a MSTrace.  The start or end time will be updated
 * and samples will be copied if they are provided.  No checking is done to
 * verify that the record matches the trace in any way.
 *
 * If whence is 1 the coverage will be added at the end of the trace,
 * whereas if whence is 2 the coverage will be added at the beginning
 * of the trace.
 *
 * Return 0 on success and -1 on error.
 ***************************************************************************/
int
mst_addspan ( MSTrace *mst, hptime_t starttime, hptime_t endtime,
	      void *datasamples, int64_t numsamples, char sampletype,
	      flag whence )
{
  int samplesize = 0;

  if ( ! mst )
    return -1;

  if ( datasamples && numsamples > 0 )
    {
      if ( (samplesize = ms_samplesize(sampletype)) == 0 )
	{
	  ms_log (2, "mst_addspan(): Unrecognized sample type: '%c'\n",
		  sampletype);
	  return -1;
	}

      if ( sampletype != mst->sampletype )
	{
	  ms_log (2, "mst_addspan(): Mismatched sample type, '%c' and '%c'\n",
		  sampletype, mst->sampletype);
	  return -1;
	}

      mst->datasamples = realloc (mst->datasamples,
				  (size_t) (mst->numsamples * samplesize + numsamples * samplesize));

      if ( mst->datasamples == NULL )
	{
	  ms_log (2, "mst_addspan(): Cannot allocate memory\n");
	  return -1;
	}
    }

  /* Add samples at end of trace */
  if ( whence == 1 )
    {
      if ( datasamples && numsamples > 0 )
	{
	  memcpy ((char *)mst->datasamples + (mst->numsamples * samplesize),
		  datasamples,
		  (size_t) (numsamples * samplesize));

	  mst->numsamples += numsamples;
	}

      mst->endtime = endtime;
    }

  /* Add samples at the beginning of trace */
  else if ( whence == 2 )
    {
      if ( datasamples && numsamples > 0 )
	{
	  /* Move any samples to end of buffer */
	  if ( mst->numsamples > 0 )
	    {
	      memmove ((char *)mst->datasamples + (numsamples * samplesize),
		       mst->datasamples,
		       (size_t) (mst->numsamples * samplesize));
	    }

	  memcpy (mst->datasamples,
		  datasamples,
		  (size_t) (numsamples * samplesize));

	  mst->numsamples += numsamples;
	}

      mst->starttime = starttime;
    }

  /* Update MSTrace sample count */
  if ( numsamples > 0 )
    mst->samplecnt += numsamples;

  return 0;
} /* End of mst_addspan() */


/***************************************************************************
 * mst_addmsrtogroup:
 *
 * Add data samples from a MSRecord to a MSTrace in a MSTraceGroup by
 * searching the group for the approriate MSTrace and either adding data
 * to it or creating a new MSTrace if no match found.
 *
 * Matching traces are found using the mst_findadjacent() routine.  If
 * the dataquality flag is true the data quality bytes must also match
 * otherwise they are ignored.
 *
 * Return a pointer to the MSTrace updated or 0 on error.
 ***************************************************************************/
MSTrace *
mst_addmsrtogroup ( MSTraceGroup *mstg, MSRecord *msr, flag dataquality,
		    double timetol, double sampratetol )
{
  MSTrace *mst = 0;
  hptime_t endtime;
  flag whence;
  char dq;

  if ( ! mstg || ! msr )
    return 0;

  dq = ( dataquality ) ? msr->dataquality : 0;

  endtime = msr_endtime (msr);

  if ( endtime == HPTERROR )
    {
      ms_log (2, "mst_addmsrtogroup(): Error calculating record end time\n");
      return 0;
    }

  /* Find matching, time adjacent MSTrace */
  mst = mst_findadjacent (mstg, &whence, dq,
			  msr->network, msr->station, msr->location, msr->channel,
			  msr->samprate, sampratetol,
			  msr->starttime, endtime, timetol);

  /* If a match was found update it otherwise create a new MSTrace and
     add to end of MSTrace chain */
  if ( mst )
    {
      /* Records with no time coverage do not contribute to a trace */
      if ( msr->samplecnt <= 0 || msr->samprate <= 0.0 )
	return mst;

      if ( mst_addmsr (mst, msr, whence) )
	{
	  return 0;
	}
    }
  else
    {
      mst = mst_init (NULL);

      mst->dataquality = dq;

      strncpy (mst->network, msr->network, sizeof(mst->network));
      strncpy (mst->station, msr->station, sizeof(mst->station));
      strncpy (mst->location, msr->location, sizeof(mst->location));
      strncpy (mst->channel, msr->channel, sizeof(mst->channel));

      mst->starttime = msr->starttime;
      mst->samprate = msr->samprate;
      mst->sampletype = msr->sampletype;

      if ( mst_addmsr (mst, msr, 1) )
	{
	  mst_free (&mst);
	  return 0;
	}

      /* Link new MSTrace into the end of the chain */
      if ( ! mstg->traces )
	{
	  mstg->traces = mst;
	}
      else
	{
	  MSTrace *lasttrace = mstg->traces;

	  while ( lasttrace->next )
	    lasttrace = lasttrace->next;

	  lasttrace->next = mst;
	}

      mstg->numtraces++;
    }

  return mst;
}  /* End of mst_addmsrtogroup() */


/***************************************************************************
 * mst_addtracetogroup:
 *
 * Add a MSTrace to a MSTraceGroup at the end of the MSTrace chain.
 *
 * Return a pointer to the MSTrace added or 0 on error.
 ***************************************************************************/
MSTrace *
mst_addtracetogroup ( MSTraceGroup *mstg, MSTrace *mst )
{
  MSTrace *lasttrace;

  if ( ! mstg || ! mst )
    return 0;

  if ( ! mstg->traces )
    {
      mstg->traces = mst;
    }
  else
    {
      lasttrace = mstg->traces;

      while ( lasttrace->next )
	lasttrace = lasttrace->next;

      lasttrace->next = mst;
    }

  mst->next = 0;

  mstg->numtraces++;

  return mst;
} /* End of mst_addtracetogroup() */


/***************************************************************************
 * mst_groupheal:
 *
 * Check if traces in MSTraceGroup can be healed, if contiguous segments
 * belong together they will be merged.  This routine is only useful
 * if the trace group was assembled from segments out of time order
 * (e.g. a file of Mini-SEED records not in time order) but forming
 * contiguous time coverage.  The MSTraceGroup will be sorted using
 * mst_groupsort() before healing.
 *
 * The time tolerance and sample rate tolerance are used to determine
 * if the traces are indeed the same.  If timetol is -1.0 the default
 * tolerance of 1/2 the sample period will be used.  If samprratetol
 * is -1.0 the default tolerance check of abs(1-sr1/sr2) < 0.0001 is
 * used (defined in libmseed.h).
 *
 * Return number of trace mergings on success otherwise -1 on error.
 ***************************************************************************/
int
mst_groupheal ( MSTraceGroup *mstg, double timetol, double sampratetol )
{
  int mergings = 0;
  MSTrace *curtrace = 0;
  MSTrace *nexttrace = 0;
  MSTrace *searchtrace = 0;
  MSTrace *prevtrace = 0;
  int8_t merged = 0;
  double postgap, pregap, delta;

  if ( ! mstg )
    return -1;

  /* Sort MSTraceGroup before any healing */
  if ( mst_groupsort (mstg, 1) )
    return -1;

  curtrace = mstg->traces;

  while ( curtrace )
    {
      nexttrace = mstg->traces;
      prevtrace = mstg->traces;

      while ( nexttrace )
	{
	  searchtrace = nexttrace;
	  nexttrace = searchtrace->next;

	  /* Do not process the same MSTrace we are trying to match */
	  if ( searchtrace == curtrace )
	    {
	      prevtrace = searchtrace;
	      continue;
	    }

	  /* Check if this trace matches the curtrace */
	  if ( strcmp (searchtrace->network, curtrace->network) ||
	       strcmp (searchtrace->station, curtrace->station) ||
	       strcmp (searchtrace->location, curtrace->location) ||
	       strcmp (searchtrace->channel, curtrace->channel) )
	    {
	      prevtrace = searchtrace;
	      continue;
	    }

	  /* Perform default samprate tolerance check if requested */
	  if ( sampratetol == -1.0 )
	    {
	      if ( ! MS_ISRATETOLERABLE (searchtrace->samprate, curtrace->samprate) )
		{
		  prevtrace = searchtrace;
		  continue;
		}
	    }
	  /* Otherwise check against the specified sample rates tolerance */
	  else if ( ms_dabs(searchtrace->samprate - curtrace->samprate) > sampratetol )
	    {
	      prevtrace = searchtrace;
	      continue;
	    }

	  merged = 0;

	  /* post/pregap are negative when searchtrace overlaps curtrace
	   * segment and positive when there is a time gap.
	   */
	  delta = ( curtrace->samprate ) ? (1.0 / curtrace->samprate) : 0.0;

	  postgap = ((double)(searchtrace->starttime - curtrace->endtime)/HPTMODULUS) - delta;

	  pregap = ((double)(curtrace->starttime - searchtrace->endtime)/HPTMODULUS) - delta;

	  /* Calculate default time tolerance (1/2 sample period) if needed */
	  if ( timetol == -1.0 )
	    timetol = 0.5 * delta;

	  /* Fits right at the end of curtrace */
	  if ( ms_dabs(postgap) <= timetol )
	    {
	      /* Merge searchtrace with curtrace */
	      mst_addspan (curtrace, searchtrace->starttime, searchtrace->endtime,
			   searchtrace->datasamples, searchtrace->numsamples,
			   searchtrace->sampletype, 1);

	      /* If no data is present, make sure sample count is updated */
	      if ( searchtrace->numsamples <= 0 )
		curtrace->samplecnt += searchtrace->samplecnt;

	      /* If qualities do not match reset the indicator */
	      if (curtrace->dataquality != searchtrace->dataquality)
		curtrace->dataquality = 0;

	      merged = 1;
	    }

	  /* Fits right at the beginning of curtrace */
	  else if ( ms_dabs(pregap) <= timetol )
	    {
	      /* Merge searchtrace with curtrace */
	      mst_addspan (curtrace, searchtrace->starttime, searchtrace->endtime,
			   searchtrace->datasamples, searchtrace->numsamples,
			   searchtrace->sampletype, 2);

	      /* If no data is present, make sure sample count is updated */
	      if ( searchtrace->numsamples <= 0 )
		curtrace->samplecnt += searchtrace->samplecnt;

	      /* If qualities do not match reset the indicator */
	      if (curtrace->dataquality != searchtrace->dataquality)
		curtrace->dataquality = 0;

	      merged = 1;
	    }

	  /* If searchtrace was merged with curtrace remove it from the chain */
	  if ( merged )
	    {
	      /* Re-link trace chain and free searchtrace */
	      if ( searchtrace == mstg->traces )
		mstg->traces = nexttrace;
	      else
		prevtrace->next = nexttrace;

	      mst_free (&searchtrace);

	      mstg->numtraces--;
	      mergings++;
	    }
	  else
	    {
	      prevtrace = searchtrace;
	    }
	}

      curtrace = curtrace->next;
    }

  return mergings;
}  /* End of mst_groupheal() */


/***************************************************************************
 * mst_groupsort:
 *
 * Sort a MSTraceGroup using a mergesort algorithm.  MSTrace entries
 * are compared using the mst_groupsort_cmp() function.
 *
 * The mergesort implementation was inspired by the listsort function
 * published and copyright 2001 by Simon Tatham.
 *
 * Return 0 on success and -1 on error.
 ***************************************************************************/
int
mst_groupsort ( MSTraceGroup *mstg, flag quality )
{
  MSTrace *p, *q, *e, *top, *tail;
  int nmerges;
  int insize, psize, qsize, i;

  if ( ! mstg )
    return -1;

  if ( ! mstg->traces )
    return 0;

  top = mstg->traces;
  insize = 1;

  for (;;)
    {
      p = top;
      top = NULL;
      tail = NULL;

      nmerges = 0;  /* count number of merges we do in this pass */

      while ( p )
        {
          nmerges++;  /* there exists a merge to be done */

          /* step `insize' places along from p */
          q = p;
          psize = 0;
          for (i = 0; i < insize; i++)
            {
              psize++;
              q = q->next;
              if ( ! q )
                break;
            }

          /* if q hasn't fallen off end, we have two lists to merge */
          qsize = insize;

          /* now we have two lists; merge them */
          while ( psize > 0 || (qsize > 0 && q) )
            {
              /* decide whether next element of merge comes from p or q */
              if ( psize == 0 )
                {  /* p is empty; e must come from q. */
                  e = q; q = q->next; qsize--;
                }
              else if ( qsize == 0 || ! q )
                {  /* q is empty; e must come from p. */
                  e = p; p = p->next; psize--;
                }
              else if ( mst_groupsort_cmp (p, q, quality) <= 0 )
                {  /* First element of p is lower (or same), e must come from p. */
                  e = p; p = p->next; psize--;
                }
              else
                {  /* First element of q is lower; e must come from q. */
                  e = q; q = q->next; qsize--;
                }

              /* add the next element to the merged list */
              if ( tail )
                tail->next = e;
              else
                top = e;

              tail = e;
            }

          /* now p has stepped `insize' places along, and q has too */
          p = q;
        }

      tail->next = NULL;

      /* If we have done only one merge, we're finished. */
      if ( nmerges <= 1 )   /* allow for nmerges==0, the empty list case */
        {
          mstg->traces = top;

	  return 0;
        }

      /* Otherwise repeat, merging lists twice the size */
      insize *= 2;
    }
} /* End of mst_groupsort() */


/***************************************************************************
 * mst_groupsort_cmp:
 *
 * Compare two MSTrace entities for the purposes of sorting a
 * MSTraceGroup.  Criteria for MSTrace comparison are (in order of
 * testing): source name, start time, descending endtime (longest
 * trace first) and sample rate.
 *
 * Return 1 if mst1 is "greater" than mst2, otherwise return 0.
 ***************************************************************************/
static int
mst_groupsort_cmp ( MSTrace *mst1, MSTrace *mst2, flag quality )
{
  char src1[50], src2[50];
  int strcmpval;

  if ( ! mst1 || ! mst2 )
    return -1;

  mst_srcname (mst1, src1, quality);
  mst_srcname (mst2, src2, quality);

  strcmpval = strcmp (src1, src2);

  /* If the source names do not match make sure the "greater" string is 2nd,
   * otherwise, if source names do match, make sure the later start time is 2nd
   * otherwise, if start times match, make sure the earlier end time is 2nd
   * otherwise, if end times match, make sure the highest sample rate is 2nd
   */
  if ( strcmpval > 0 )
    {
      return 1;
    }
  else if ( strcmpval == 0 )
    {
      if ( mst1->starttime > mst2->starttime )
	{
	  return 1;
	}
      else if ( mst1->starttime == mst2->starttime )
	{
	  if ( mst1->endtime < mst2->endtime )
	    {
	      return 1;
	    }
	  else if ( mst1->endtime == mst2->endtime )
	    {
	      if ( ! MS_ISRATETOLERABLE (mst1->samprate, mst2->samprate) &&
		   mst1->samprate > mst2->samprate )
		{
		  return 1;
		}
	    }
	}
    }

  return 0;
} /* End of mst_groupsort_cmp() */


/***************************************************************************
 * mst_srcname:
 *
 * Generate a source name string for a specified MSTrace in the
 * format: 'NET_STA_LOC_CHAN[_QUAL]'.  The quality is added to the
 * srcname if the quality flag argument is 1 and mst->dataquality is
 * not zero.  The passed srcname must have enough room for the
 * resulting string.
 *
 * Returns a pointer to the resulting string or NULL on error.
 ***************************************************************************/
char *
mst_srcname (MSTrace *mst, char *srcname, flag quality)
{
  char *src = srcname;
  char *cp = srcname;

  if ( ! mst || ! srcname )
    return NULL;

  /* Build the source name string */
  cp = mst->network;
  while ( *cp ) { *src++ = *cp++; }
  *src++ = '_';
  cp = mst->station;
  while ( *cp ) { *src++ = *cp++; }
  *src++ = '_';
  cp = mst->location;
  while ( *cp ) { *src++ = *cp++; }
  *src++ = '_';
  cp = mst->channel;
  while ( *cp ) { *src++ = *cp++; }

  if ( quality && mst->dataquality )
    {
      *src++ = '_';
      *src++ = mst->dataquality;
    }

  *src = '\0';

  return srcname;
} /* End of mst_srcname() */


/***************************************************************************
 * mst_printtracelist:
 *
 * Print trace list summary information for the specified MSTraceGroup.
 *
 * By default only print the srcname, starttime and endtime for each
 * trace.  If details is greater than 0 include the sample rate,
 * number of samples and a total trace count.  If gaps is greater than
 * 0 and the previous trace matches (srcname & samprate) include the
 * gap between the endtime of the last trace and the starttime of the
 * current trace.
 *
 * The timeformat flag can either be:
 * 0 : SEED time format (year, day-of-year, hour, min, sec)
 * 1 : ISO time format (year, month, day, hour, min, sec)
 * 2 : Epoch time, seconds since the epoch
 ***************************************************************************/
void
mst_printtracelist ( MSTraceGroup *mstg, flag timeformat,
		     flag details, flag gaps )
{
  MSTrace *mst = 0;
  char srcname[50];
  char prevsrcname[50];
  char stime[30];
  char etime[30];
  char gapstr[20];
  flag nogap;
  double gap;
  double delta;
  double prevsamprate;
  hptime_t prevendtime;
  int tracecnt = 0;

  if ( ! mstg )
    {
      return;
    }

  mst = mstg->traces;

  /* Print out the appropriate header */
  if ( details > 0 && gaps > 0 )
    ms_log (0, "   Source                Start sample             End sample        Gap  Hz  Samples\n");
  else if ( details <= 0 && gaps > 0 )
    ms_log (0, "   Source                Start sample             End sample        Gap\n");
  else if ( details > 0 && gaps <= 0 )
    ms_log (0, "   Source                Start sample             End sample        Hz  Samples\n");
  else
    ms_log (0, "   Source                Start sample             End sample\n");

  prevsrcname[0] = '\0';
  prevsamprate = -1.0;
  prevendtime = 0;

  while ( mst )
    {
      mst_srcname (mst, srcname, 1);

      /* Create formatted time strings */
      if ( timeformat == 2 )
	{
	  snprintf (stime, sizeof(stime), "%.6f", (double) MS_HPTIME2EPOCH(mst->starttime) );
	  snprintf (etime, sizeof(etime), "%.6f", (double) MS_HPTIME2EPOCH(mst->endtime) );
	}
      else if ( timeformat == 1 )
	{
	  if ( ms_hptime2isotimestr (mst->starttime, stime, 1) == NULL )
	    ms_log (2, "Cannot convert trace start time for %s\n", srcname);

	  if ( ms_hptime2isotimestr (mst->endtime, etime, 1) == NULL )
	    ms_log (2, "Cannot convert trace end time for %s\n", srcname);
	}
      else
	{
	  if ( ms_hptime2seedtimestr (mst->starttime, stime, 1) == NULL )
	    ms_log (2, "Cannot convert trace start time for %s\n", srcname);

	  if ( ms_hptime2seedtimestr (mst->endtime, etime, 1) == NULL )
	    ms_log (2, "Cannot convert trace end time for %s\n", srcname);
	}

      /* Print trace info at varying levels */
      if ( gaps > 0 )
	{
	  gap = 0.0;
	  nogap = 0;

	  if ( ! strcmp (prevsrcname, srcname) && prevsamprate != -1.0 &&
	       MS_ISRATETOLERABLE (prevsamprate, mst->samprate) )
	    gap = (double) (mst->starttime - prevendtime) / HPTMODULUS;
	  else
	    nogap = 1;

	  /* Check that any overlap is not larger than the trace coverage */
	  if ( gap < 0.0 )
	    {
	      delta = ( mst->samprate ) ? (1.0 / mst->samprate) : 0.0;

	      if ( (gap * -1.0) > (((double)(mst->endtime - mst->starttime)/HPTMODULUS) + delta) )
		gap = -(((double)(mst->endtime - mst->starttime)/HPTMODULUS) + delta);
	    }

	  /* Fix up gap display */
	  if ( nogap )
	    snprintf (gapstr, sizeof(gapstr), " == ");
	  else if ( gap >= 86400.0 || gap <= -86400.0 )
	    snprintf (gapstr, sizeof(gapstr), "%-3.1fd", (gap / 86400));
	  else if ( gap >= 3600.0 || gap <= -3600.0 )
	    snprintf (gapstr, sizeof(gapstr), "%-3.1fh", (gap / 3600));
	  else if ( gap == 0.0 )
	    snprintf (gapstr, sizeof(gapstr), "-0  ");
	  else
	    snprintf (gapstr, sizeof(gapstr), "%-4.4g", gap);

	  if ( details <= 0 )
	    ms_log (0, "%-17s %-24s %-24s %-4s\n",
		    srcname, stime, etime, gapstr);
	  else
	    ms_log (0, "%-17s %-24s %-24s %-s %-3.3g %-lld\n",
		    srcname, stime, etime, gapstr, mst->samprate, (long long int)mst->samplecnt);
	}
      else if ( details > 0 && gaps <= 0 )
	ms_log (0, "%-17s %-24s %-24s %-3.3g %-lld\n",
		srcname, stime, etime, mst->samprate, (long long int)mst->samplecnt);
      else
	ms_log (0, "%-17s %-24s %-24s\n", srcname, stime, etime);

      if ( gaps > 0 )
	{
	  strcpy (prevsrcname, srcname);
	  prevsamprate = mst->samprate;
	  prevendtime = mst->endtime;
	}

      tracecnt++;
      mst = mst->next;
    }

  if ( tracecnt != mstg->numtraces )
    ms_log (2, "mst_printtracelist(): number of traces in trace group is inconsistent\n");

  if ( details > 0 )
    ms_log (0, "Total: %d trace segment(s)\n", tracecnt);

}  /* End of mst_printtracelist() */


/***************************************************************************
 * mst_printsynclist:
 *
 * Print SYNC trace list summary information for the specified MSTraceGroup.
 *
 * The SYNC header line will be created using the supplied dccid, if
 * the pointer is NULL the string "DCC" will be used instead.
 *
 * If the subsecond flag is true the segment start and end times will
 * include subsecond precision, otherwise they will be truncated to
 * integer seconds.
 *
 ***************************************************************************/
void
mst_printsynclist ( MSTraceGroup *mstg, char *dccid, flag subsecond )
{
  MSTrace *mst = 0;
  char stime[30];
  char etime[30];
  char yearday[10];
  time_t now;
  struct tm *nt;

  if ( ! mstg )
    {
      return;
    }

  /* Generate current time stamp */
  now = time (NULL);
  nt = localtime ( &now ); nt->tm_year += 1900; nt->tm_yday += 1;
  snprintf ( yearday, sizeof(yearday), "%04d,%03d", nt->tm_year, nt->tm_yday);

  /* Print SYNC header line */
  ms_log (0, "%s|%s\n", (dccid)?dccid:"DCC", yearday);

  /* Loope through trace list */
  mst = mstg->traces;
  while ( mst )
    {
      ms_hptime2seedtimestr (mst->starttime, stime, subsecond);
      ms_hptime2seedtimestr (mst->endtime, etime, subsecond);

      /* Print SYNC line */
      ms_log (0, "%s|%s|%s|%s|%s|%s||%.10g|%lld|||||||%s\n",
	      mst->network, mst->station, mst->location, mst->channel,
	      stime, etime, mst->samprate, (long long int)mst->samplecnt,
	      yearday);

      mst = mst->next;
    }
}  /* End of mst_printsynclist() */


/***************************************************************************
 * mst_printgaplist:
 *
 * Print gap/overlap list summary information for the specified
 * MSTraceGroup.  Overlaps are printed as negative gaps.  The trace
 * summary information in the MSTraceGroup is logically inverted so gaps
 * for like channels are identified.
 *
 * If mingap and maxgap are not NULL their values will be enforced and
 * only gaps/overlaps matching their implied criteria will be printed.
 *
 * The timeformat flag can either be:
 * 0 : SEED time format (year, day-of-year, hour, min, sec)
 * 1 : ISO time format (year, month, day, hour, min, sec)
 * 2 : Epoch time, seconds since the epoch
 ***************************************************************************/
void
mst_printgaplist (MSTraceGroup *mstg, flag timeformat,
		  double *mingap, double *maxgap)
{
  MSTrace *mst;
  char src1[50], src2[50];
  char time1[30], time2[30];
  char gapstr[30];
  double gap;
  double delta;
  double nsamples;
  flag printflag;
  int gapcnt = 0;

  if ( ! mstg )
    return;

  if ( ! mstg->traces )
    return;

  mst = mstg->traces;

  ms_log (0, "   Source                Last Sample              Next Sample       Gap  Samples\n");

  while ( mst->next )
    {
      mst_srcname (mst, src1, 1);
      mst_srcname (mst->next, src2, 1);

      if ( ! strcmp (src1, src2) )
	{
	  /* Skip MSTraces with 0 sample rate, usually from SOH records */
	  if ( mst->samprate == 0.0 )
	    {
	      mst = mst->next;
	      continue;
	    }

	  /* Check that sample rates match using default tolerance */
	  if ( ! MS_ISRATETOLERABLE (mst->samprate, mst->next->samprate) )
	    {
	      ms_log (2, "%s Sample rate changed! %.10g -> %.10g\n",
		      src1, mst->samprate, mst->next->samprate );
	    }

	  gap = (double) (mst->next->starttime - mst->endtime) / HPTMODULUS;

	  /* Check that any overlap is not larger than the trace coverage */
	  if ( gap < 0.0 )
	    {
	      delta = ( mst->next->samprate ) ? (1.0 / mst->next->samprate) : 0.0;

	      if ( (gap * -1.0) > (((double)(mst->next->endtime - mst->next->starttime)/HPTMODULUS) + delta) )
		gap = -(((double)(mst->next->endtime - mst->next->starttime)/HPTMODULUS) + delta);
	    }

	  printflag = 1;

	  /* Check gap/overlap criteria */
	  if ( mingap )
	    if ( gap < *mingap )
	      printflag = 0;

	  if ( maxgap )
	    if ( gap > *maxgap )
	      printflag = 0;

	  if ( printflag )
	    {
	      nsamples = ms_dabs(gap) * mst->samprate;

	      if ( gap > 0.0 )
		nsamples -= 1.0;
	      else
		nsamples += 1.0;

	      /* Fix up gap display */
	      if ( gap >= 86400.0 || gap <= -86400.0 )
		snprintf (gapstr, sizeof(gapstr), "%-3.1fd", (gap / 86400));
	      else if ( gap >= 3600.0 || gap <= -3600.0 )
		snprintf (gapstr, sizeof(gapstr), "%-3.1fh", (gap / 3600));
	      else if ( gap == 0.0 )
		snprintf (gapstr, sizeof(gapstr), "-0  ");
	      else
		snprintf (gapstr, sizeof(gapstr), "%-4.4g", gap);

	      /* Create formatted time strings */
	      if ( timeformat == 2 )
		{
		  snprintf (time1, sizeof(time1), "%.6f", (double) MS_HPTIME2EPOCH(mst->endtime) );
		  snprintf (time2, sizeof(time2), "%.6f", (double) MS_HPTIME2EPOCH(mst->next->starttime) );
		}
	      else if ( timeformat == 1 )
		{
		  if ( ms_hptime2isotimestr (mst->endtime, time1, 1) == NULL )
		    ms_log (2, "Cannot convert trace end time for %s\n", src1);

		  if ( ms_hptime2isotimestr (mst->next->starttime, time2, 1) == NULL )
		    ms_log (2, "Cannot convert next trace start time for %s\n", src1);
		}
	      else
		{
		  if ( ms_hptime2seedtimestr (mst->endtime, time1, 1) == NULL )
		    ms_log (2, "Cannot convert trace end time for %s\n", src1);

		  if ( ms_hptime2seedtimestr (mst->next->starttime, time2, 1) == NULL )
		    ms_log (2, "Cannot convert next trace start time for %s\n", src1);
		}

	      ms_log (0, "%-17s %-24s %-24s %-4s %-.8g\n",
		      src1, time1, time2, gapstr, nsamples);

	      gapcnt++;
	    }
	}

      mst = mst->next;
    }

  ms_log (0, "Total: %d gap(s)\n", gapcnt);

}  /* End of mst_printgaplist() */


/***************************************************************************
 * mst_pack:
 *
 * Pack MSTrace data into Mini-SEED records using the specified record
 * length, encoding format and byte order.  The datasamples array and
 * numsamples field will be adjusted (reduced) based on how many
 * samples were packed.
 *
 * As each record is filled and finished they are passed to
 * record_handler which expects 1) a char * to the record, 2) the
 * length of the record and 3) a pointer supplied by the original
 * caller containing optional private data (handlerdata).  It is the
 * responsibility of record_handler to process the record, the memory
 * will be re-used or freed when record_handler returns.
 *
 * If the flush flag is > 0 all of the data will be packed into data
 * records even though the last one will probably not be filled.
 *
 * If the mstemplate argument is not NULL it will be used as the
 * template for the packed Mini-SEED records.  Otherwise a new
 * MSRecord will be initialized and populated from values in the
 * MSTrace.  The reclen, encoding and byteorder arguments take
 * precedence over those in the template.  The start time, sample
 * rate, datasamples, numsamples and sampletype values from the
 * template will be preserved.
 *
 * Returns the number of records created on success and -1 on error.
 ***************************************************************************/
int
mst_pack ( MSTrace *mst, void (*record_handler) (char *, int, void *),
	   void *handlerdata, int reclen, flag encoding, flag byteorder,
	   int64_t *packedsamples, flag flush, flag verbose,
	   MSRecord *mstemplate )
{
  MSRecord *msr;
  char srcname[50];
  int trpackedrecords;
  int64_t trpackedsamples;
  int samplesize;
  int64_t bufsize;

  hptime_t preservestarttime = 0;
  double preservesamprate = 0.0;
  void *preservedatasamples = 0;
  int64_t preservenumsamples = 0;
  char preservesampletype = 0;
  StreamState *preserveststate = 0;

  if ( packedsamples )
    *packedsamples = 0;

  /* Allocate stream processing state space if needed */
  if ( ! mst->ststate )
    {
      mst->ststate = (StreamState *) malloc (sizeof(StreamState));
      if ( ! mst->ststate )
        {
          ms_log (2, "mst_pack(): Could not allocate memory for StreamState\n");
          return -1;
        }
      memset (mst->ststate, 0, sizeof(StreamState));
    }

  if ( mstemplate )
    {
      msr = mstemplate;

      preservestarttime = msr->starttime;
      preservesamprate = msr->samprate;
      preservedatasamples = msr->datasamples;
      preservenumsamples = msr->numsamples;
      preservesampletype = msr->sampletype;
      preserveststate = msr->ststate;
    }
  else
    {
      msr = msr_init (NULL);

      if ( msr == NULL )
	{
	  ms_log (2, "mst_pack(): Error initializing msr\n");
	  return -1;
	}

      msr->dataquality = 'D';
      strcpy (msr->network, mst->network);
      strcpy (msr->station, mst->station);
      strcpy (msr->location, mst->location);
      strcpy (msr->channel, mst->channel);
    }

  /* Setup MSRecord template for packing */
  msr->reclen = reclen;
  msr->encoding = encoding;
  msr->byteorder = byteorder;

  msr->starttime = mst->starttime;
  msr->samprate = mst->samprate;
  msr->datasamples = mst->datasamples;
  msr->numsamples = mst->numsamples;
  msr->sampletype = mst->sampletype;
  msr->ststate = mst->ststate;

  /* Sample count sanity check */
  if ( mst->samplecnt != mst->numsamples )
    {
      ms_log (2, "mst_pack(): Sample counts do not match, abort\n");
      return -1;
    }

  /* Pack data */
  trpackedrecords = msr_pack (msr, record_handler, handlerdata, &trpackedsamples, flush, verbose);

  if ( verbose > 1 )
    {
      ms_log (1, "Packed %d records for %s trace\n", trpackedrecords, mst_srcname (mst, srcname, 1));
    }

  /* Adjust MSTrace start time, data array and sample count */
  if ( trpackedsamples > 0 )
    {
      /* The new start time was calculated my msr_pack */
      mst->starttime = msr->starttime;

      samplesize = ms_samplesize (mst->sampletype);
      bufsize = (mst->numsamples - trpackedsamples) * samplesize;

      if ( bufsize )
	{
	  memmove (mst->datasamples,
		   (char *) mst->datasamples + (trpackedsamples * samplesize),
		   (size_t) bufsize);

	  mst->datasamples = realloc (mst->datasamples, (size_t) bufsize);

	  if ( mst->datasamples == NULL )
	    {
	      ms_log (2, "mst_pack(): Cannot (re)allocate datasamples buffer\n");
	      return -1;
	    }
	}
      else
	{
	  if ( mst->datasamples )
	    free (mst->datasamples);
	  mst->datasamples = 0;
	}

      mst->samplecnt -= trpackedsamples;
      mst->numsamples -= trpackedsamples;
    }

  /* Reinstate preserved values if a template was used */
  if ( mstemplate )
    {
      msr->starttime = preservestarttime;
      msr->samprate = preservesamprate;
      msr->datasamples = preservedatasamples;
      msr->numsamples = preservenumsamples;
      msr->sampletype = preservesampletype;
      msr->ststate = preserveststate;
    }
  else
    {
      msr->datasamples = 0;
      msr->ststate = 0;
      msr_free (&msr);
    }

  if ( packedsamples )
    *packedsamples = trpackedsamples;

  return trpackedrecords;
}  /* End of mst_pack() */


/***************************************************************************
 * mst_packgroup:
 *
 * Pack MSTraceGroup data into Mini-SEED records by calling mst_pack()
 * for each MSTrace in the group.
 *
 * Returns the number of records created on success and -1 on error.
 ***************************************************************************/
int
mst_packgroup ( MSTraceGroup *mstg, void (*record_handler) (char *, int, void *),
		void *handlerdata, int reclen, flag encoding, flag byteorder,
		int64_t *packedsamples, flag flush, flag verbose,
		MSRecord *mstemplate )
{
  MSTrace *mst;
  int trpackedrecords = 0;
  int64_t trpackedsamples = 0;
  char srcname[50];

  if ( ! mstg )
    {
      return -1;
    }

  if ( packedsamples )
    *packedsamples = 0;

  mst = mstg->traces;

  while ( mst )
    {
      if ( mst->numsamples <= 0 )
	{
	  if ( verbose > 1 )
	    {
	      mst_srcname (mst, srcname, 1);
	      ms_log (1, "No data samples for %s, skipping\n", srcname);
	    }
	}
      else
	{
	  trpackedrecords += mst_pack (mst, record_handler, handlerdata, reclen,
				     encoding, byteorder, &trpackedsamples, flush,
				     verbose, mstemplate);

	  if ( trpackedrecords == -1 )
	    break;

	  if ( packedsamples )
	    *packedsamples += trpackedsamples;
	}

      mst = mst->next;
    }

  return trpackedrecords;
}  /* End of mst_packgroup() */
