#include        <stdio.h>
#include        "spline.h"

#include        "utility.h"

void spline_print(FILE *fp, spline *s)
{
  int        i;
  fprintf(fp,"Spline: ");
  if ( s == NULL )
  {
    fprintf(fp,"NULL\n");
    return;
  }
  fprintf(fp,"length: %d\n", s->length);
  fprintf(fp,"xlist: ");
  for ( i = 0; i < s->length; i++ )
    fprintf(fp,"%20.14g ", s->xlist[i]);
  fprintf(fp,"\n");
  fprintf(fp,"ylist: ");
  for ( i = 0; i < s->length; i++ )
    fprintf(fp,"%20.14g ", s->ylist[i]);
  fprintf(fp,"\n");
  fprintf(fp,"Mlist: ");
  for ( i = 0; i < s->length; i++ )
    fprintf(fp,"%20.14g ", s->Mlist[i]);
  fprintf(fp,"\n");
}


void spline_dump(FILE *fp, spline *s)
{
  int        i;
  fprintf(fp,"Spline: ");
  if ( s == NULL )
  {
    fprintf(fp,"NULL\n");
    return;
  }
  fprintf(fp,"@ %p ", s);
  fprintf(fp,"length: %d\n", s->length);
  fprintf(fp,"xlist @ %p: ", s->xlist);
  for ( i = 0; i < s->length; i++ )
    fprintf(fp,"%20.14g ", s->xlist[i]);
  fprintf(fp,"\n");
  fprintf(fp,"ylist @ %p: ", s->ylist);
  for ( i = 0; i < s->length; i++ )
    fprintf(fp,"%20.14g ", s->ylist[i]);
  fprintf(fp,"\n");
  fprintf(fp,"Mlist @ %p: ", s->Mlist);
  for ( i = 0; i < s->length; i++ )
    fprintf(fp,"%20.14g ", s->Mlist[i]);
  fprintf(fp,"\n");
}


spline *spline_read(FILE *fp)
{
  int        i, length;
  spline *s;

  /* read header */
  if ( fscanf(fp," Spline: length:%u",&length) < 1)
    return NULL;
  s = spline_create(length);

  /* read xlist */
  fscanf(fp," xlist:");
  for ( i = 0; i < length; i++ )
    if ( fscanf(fp,"%lf",&(s->xlist[i])) < 1 )
      spline_error("spline_read: fscanf error");
  /* read ylist */
  fscanf(fp," ylist:");
  for ( i = 0; i < length; i++ )
    if ( fscanf(fp,"%lf",&(s->ylist[i])) < 1 )
      spline_error("spline_read: fscanf error");
  /* read Mlist */
  fscanf(fp," Mlist:");
  for ( i = 0; i < length; i++ )
    if ( fscanf(fp,"%lf",&(s->Mlist[i])) < 1 )
      spline_error("spline_read: fscanf error");

  return s;
}


spline        *spline_create(int length)
{
  spline        *s;

  if ( length <= 0 )
    spline_error("spline_create: length <= 0");

  if ( (s=NEW(spline)) == NULL )
    spline_error("spline_create: allocation error");
  if ( (s->xlist = NEW_ARRAY(length,double)) == NULL )
    spline_error("spline_create: allocation error");
  if ( (s->ylist = NEW_ARRAY(length,double)) == NULL )
    spline_error("spline_create: allocation error");
  if ( (s->Mlist = NEW_ARRAY(length,double)) == NULL )
    spline_error("spline_create: allocation error");

  s->length = length;

  return s;
}


void spline_destroy(spline *s)
{
  if ( s == NULL )
    return;    /* nothing to do */

  if ( s->Mlist != NULL )
    free(s->Mlist);
  if ( s->ylist != NULL )
    free(s->ylist);
  if ( s->xlist != NULL )
    free(s->xlist);

  free(s);
}


spline        *spline_setxy(spline *s, 
                      double x[], double y[], int length)
{
  int        i;

  if ( s == NULL )
    spline_error("spline_setxy: NULL pointer passed for spline");
  if ( s->length != length )
    spline_error("spline_setxy: length mismatch");

  /* check that x[] is in strictly increaseing order */
  for ( i = 0; i < length-1; i++ )
    if ( x[i+1] <= x[i] )
      spline_error("spline_setxy: x input not strictly increasing");

  for ( i = 0; i < length; i++ )
    {
      s->xlist[i] = x[i];
      s->ylist[i] = y[i];
    }

  return s;
}