/*
 * mdengine.c
 *
 * Summary:  Engine API routines.
 *
 * Author:  David Hardy
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
/* #include <dlfcn.h> */
#include "ihash.h"
#include "mdsim.h"
#if 0
#include "mdclient.h"
#endif


  /* type info for MDAPI general data types */

typedef struct MD_Init_Type_Tag {
  const char *name;
  MD_Int elemsize;
  MD_Int typenum;
  MD_Int subtype[19];
} MD_Init_Type;

static const MD_Init_Type INIT_TYPE_TABLE[] = {
  { "MD_Char", 1, MD_CHAR, { MD_CHAR, -1 } },
  { "MD_Int", 4, MD_INT, { MD_INT, -1 } },
  { "MD_Float", 4, MD_FLOAT, { MD_FLOAT, -1 } },
  { "MD_Double", 8, MD_DOUBLE, { MD_DOUBLE, -1 } },
  { "MD_Fvec", 12, MD_FVEC, { MD_FVEC, -1 } },
  { "MD_Dvec", 24, MD_DVEC, { MD_DVEC, -1 } },
  { "MD_Short_String", 8, MD_SHORT_STRING, { MD_SHORT_STRING, -1 } },
  { "MD_String", 64, MD_STRING, { MD_STRING, -1 } },
  { "MD_Atom_Param", 40, MD_ATOM_PARAM,
    { MD_DOUBLE, MD_DOUBLE, MD_DOUBLE, MD_DOUBLE, MD_SHORT_STRING, -1 } },
  { "MD_Bond_Param", 32, MD_BOND_PARAM,
    { MD_DOUBLE, MD_DOUBLE, MD_SHORT_STRING, MD_SHORT_STRING, -1 } },
  { "MD_Angle_Param", 56, MD_ANGLE_PARAM,
    { MD_DOUBLE, MD_DOUBLE, MD_DOUBLE, MD_DOUBLE,
      MD_SHORT_STRING, MD_SHORT_STRING, MD_SHORT_STRING, -1 } },
  { "MD_Dihed_Param", 120, MD_DIHED_PARAM,
    { MD_DOUBLE, MD_DOUBLE, MD_DOUBLE, MD_DOUBLE,
      MD_DOUBLE, MD_DOUBLE, MD_DOUBLE, MD_DOUBLE,
      MD_INT, MD_INT, MD_INT, MD_INT, MD_INT, MD_INT,
      MD_SHORT_STRING, MD_SHORT_STRING, MD_SHORT_STRING, MD_SHORT_STRING, -1}},
  { "MD_Impr_Param", 120, MD_IMPR_PARAM,
    { MD_DOUBLE, MD_DOUBLE, MD_DOUBLE, MD_DOUBLE,
      MD_DOUBLE, MD_DOUBLE, MD_DOUBLE, MD_DOUBLE,
      MD_INT, MD_INT, MD_INT, MD_INT, MD_INT, MD_INT,
      MD_SHORT_STRING, MD_SHORT_STRING, MD_SHORT_STRING, MD_SHORT_STRING, -1}},
  { "MD_Nbfix_Param", 56, MD_NBFIX_PARAM,
    { MD_DOUBLE, MD_DOUBLE, MD_DOUBLE, MD_DOUBLE,
      MD_SHORT_STRING, MD_SHORT_STRING, MD_INT, MD_INT, -1 } },
  { "MD_Atom", 40, MD_ATOM,
    { MD_DOUBLE, MD_DOUBLE, MD_INT, MD_INT,
      MD_SHORT_STRING, MD_SHORT_STRING, -1 } },
  { "MD_Bond", 12, MD_BOND, { MD_INT, MD_INT, MD_INT, -1 } },
  { "MD_Angle", 16, MD_ANGLE, { MD_INT, MD_INT, MD_INT, MD_INT, -1 } },
  { "MD_Dihed", 20, MD_DIHED, { MD_INT, MD_INT, MD_INT, MD_INT, MD_INT, -1 } },
  { "MD_Impr", 20, MD_IMPR, { MD_INT, MD_INT, MD_INT, MD_INT, MD_INT, -1 } },
  { "MD_Excl", 8, MD_EXCL, { MD_INT, MD_INT, -1 } },
};


  /* prototypes for static functions */

static MD_Int validate_type(MD_Sim *, const MD_Int *subtype, MD_Int len);
static MD_Int add_type(MD_Sim *, const char *name,
                       const MD_Int *subtype, MD_Int len, MD_Int numbytes);


/*****************************************************************************
 *
 * Initialization routines during create_engine.
 *
 *****************************************************************************/

MD_Errcode MD_set_version(MD_Sim *sim, MD_Int version)
{
  MD_Int len, typenum, i;

  if (sim->version != 0 || version != MD_API_VERSION) {
    ERR(sim, MD_ERR_VERSION, "MD_set_version");
    return MD_FAIL;
  }
  sim->version = version;

  /* use this moment to initialize all MDAPI general data types */
  /* add primitive types first */
  for (i = INDEX(MD_CHAR);  i <= INDEX(MD_STRING);  i++) {
    len = 0;
    while (INIT_TYPE_TABLE[i].subtype[len] != -1)  len++;
    sim->flipflags = INIT_TYPE_TABLE[i].typenum & MD_FLIP_MASK;
    typenum = add_type(sim, INIT_TYPE_TABLE[i].name, 
                       INIT_TYPE_TABLE[i].subtype, len,
                       INIT_TYPE_TABLE[i].elemsize);
    if (typenum != INIT_TYPE_TABLE[i].typenum) {
      if (sim->errnum == 0) {
        ERR(sim, MD_ERR_TYPE, "MD_set_version call to add_type");
      }
      return MD_FAIL;
    }
  }

  /* call MD_new_type to add derived types next */
  for ( ;  i < ARRLEN(INIT_TYPE_TABLE);  i++) {
    len = 0;
    while (INIT_TYPE_TABLE[i].subtype[len] != -1)  len++;
    typenum = MD_new_type(sim, INIT_TYPE_TABLE[i].name,
                          INIT_TYPE_TABLE[i].subtype, len,
                          INIT_TYPE_TABLE[i].elemsize);
    if (typenum != INIT_TYPE_TABLE[i].typenum) {
      if (sim->errnum == 0) {
        ERR(sim, MD_ERR_TYPE, "MD_set_version: MDAPI typenum mismatch");
      }
    }
  }

  return 0;
}


MD_Errcode MD_set_engine(MD_Sim *sim, void *engine)
{
  /* check on incorrect use */
  if (sim->version != MD_API_VERSION) {
    ERR(sim, MD_ERR_VERSION, "MD_set_engine");
    return MD_FAIL;
  }
  else if (sim->iscreated) {
    ERR(sim, MD_ERR_USE, "MD_set_engine");
    return MD_FAIL;
  }
  sim->engine = engine;
  return 0;
}


MD_Int MD_new_data(MD_Sim *sim, void *ppdata, MD_Int *plen, MD_Int typenum,
                   MD_Int maxlen, MD_Int access, MD_Setlen_Fptr setlen,
                   MD_Setdata_Fptr setdata, MD_Getdata_Fptr getdata)
{
  MD_Int idnum;
  MD_Int n, i;
  char buf[20];

  /* check on incorrect use */
  if (sim->version != MD_API_VERSION) {
    ERR(sim, MD_ERR_VERSION, "MD_new_data");
    return MD_FAIL;
  }
  else if (sim->state != STATE_CREATE) {
    ERR(sim, MD_ERR_USE, "MD_new_data");
    return MD_FAIL;
  }

  idnum = sim->datalen;
  if (idnum == sim->datamax) {
    /* need to resize data array */
    void *tmp = realloc(sim->data, 2 * sim->datamax * sizeof(MD_Data));
    if (tmp == NULL) {
      ERR(sim, MD_ERR_MEMALLOC, "MD_new_data call to realloc data array");
      return MD_FAIL;
    }
    sim->data = (MD_Data *) tmp;
    sim->datamax *= 2;
  }

  /* reset corresponding access flag if function pointer is NULL */
  if ((sim->data[idnum].setlen = setlen) == NULL) {
    access &= ~MD_DATA_RESIZE;
  }
  else if ((access & MD_DATA_RESIZE) == 0) {
    sim->data[idnum].setlen = NULL;
  }

  if ((sim->data[idnum].setdata = setdata) == NULL) {
    access &= ~MD_DATA_WRITE;
  }
  else if ((access & MD_DATA_WRITE) == 0) {
    sim->data[idnum].setdata = NULL;
  }

  if ((sim->data[idnum].getdata = getdata) == NULL) {
    access &= ~MD_DATA_READ;
  }
  else if ((access & MD_DATA_READ) == 0) {
    sim->data[idnum].getdata = NULL;
  }

  sim->data[idnum].ppdata = ppdata;
  sim->data[idnum].plen = plen;
  sim->data[idnum].attr.idnum = idnum;
  sim->data[idnum].attr.typenum = typenum;
  sim->data[idnum].attr.maxlen = maxlen;
  sim->data[idnum].attr.access = access;

  /* write plen (address) as a string (ASCII hex) -- check hash table */
  sprintf(buf, "%lx", (unsigned long) plen);
  n = ihash_lookup( &(sim->plentable), buf);
  if (n >= 0) {
    /*
     * already in hash table, n is smallest index of array using
     * this shared length variable
     *
     * need to traverse the circular linked list to see where
     * to insert our own link
     */
    for (i = n;  sim->data[i].link != n;  i = sim->data[i].link) ;
    sim->data[i].link = idnum;
    sim->data[idnum].link = n;
  }
  else {
    /*
     * so this length variable is not (yet) being used by any
     * other data array, make a link to myself
     */
    sim->data[idnum].link = idnum;
    /* store string in plenstr array */
    if (sim->plenstrlen == sim->plenstrmax) {
      /* need to resize plenstr array */
      void *tmp = realloc(sim->plenstr, 2 * sim->plenstrmax * sizeof(char *));
      if (tmp == NULL) {
        ERR(sim, MD_ERR_MEMALLOC, "MD_new_data call to realloc plenstr array");
        return MD_FAIL;
      }
      sim->plenstr = (char **) tmp;
      sim->plenstrmax *= 2;
    }
    /* allocate memory to store string */
    if ((sim->plenstr[sim->plenstrlen] = strdup(buf)) == NULL) {
      ERR(sim, MD_ERR_MEMALLOC, "MD_new_data call to strdup");
      return MD_FAIL;
    }
    /* insert into hash table */
    n = ihash_insert( &(sim->plentable), sim->plenstr[sim->plenstrlen], idnum);
    sim->plenstrlen++;
    if (n != IHASH_FAIL) {
      ERR(sim, MD_ERR_MEMALLOC, "MD_new_data call to ihash_insert");
      return MD_FAIL;
    }
  }
  sim->datalen++;
  return idnum;
}


MD_Errcode MD_new_name(MD_Sim *sim, const char *name, MD_Int idnum)
{
  MD_Int retval;
  char *str, *pch;

  /* check on incorrect use */
  if (sim->version != MD_API_VERSION) {
    ERR(sim, MD_ERR_VERSION, "MD_new_name");
    return MD_FAIL;
  }
  else if (sim->state != STATE_CREATE) {
    ERR(sim, MD_ERR_USE, "MD_new_name");
    return MD_FAIL;
  }

  /* check for legal idnum value */
  if (idnum < 0 || idnum >= sim->datalen) {
    ERR(sim, MD_ERR_IDNUM, "MD_new_name");
    return MD_FAIL;
  }

  /* place name in hash table for fast lookup of idnum */
  /* make sure to convert hash key to all lower case */
  if (sim->hashkeylen == sim->hashkeymax) {
    /* need to resize hashkey array */
    void *tmp = realloc(sim->hashkey, 2 * sim->hashkeymax * sizeof(char *));
    if (tmp == NULL) {
      ERR(sim, MD_ERR_MEMALLOC, "MD_new_name call to realloc");
      return MD_FAIL;
    }
    sim->hashkey = (char **) tmp;
    sim->hashkeymax *= 2;
  }
  if ((str = sim->hashkey[sim->hashkeylen] = strdup(name)) == NULL) {
    ERR(sim, MD_ERR_MEMALLOC, "MD_new_name call to strdup");
    return MD_FAIL;
  }
  for (pch = str;  *pch != '\0';  pch++) {
    if (isupper(*pch))  *pch = tolower(*pch);
  }
  sim->hashkeylen++;

  retval = ihash_insert( &(sim->nametable), str, idnum);
  if (retval >= 0) {
    /* name is already in table */
    free(str);
    sim->hashkeylen--;
    ERR(sim, MD_ERR_NAME, "MD_new_name call to ihash_insert");
    return MD_FAIL;
  }
  else if (retval != IHASH_FAIL) {
    /* in this case, there was a memory allocation error */
    ERR(sim, MD_ERR_MEMALLOC, "MD_new_name call to ihash_insert");
    return MD_FAIL;
  }

  /* if access is not invisible, publish name in namelist */
  if ((sim->data[idnum].attr.access & MD_DATA_INVISIBLE) == 0) {
    if (sim->namelen == sim->namemax - 1) {
      /* need to resize name array */
      void *tmp = realloc(sim->name, 2 * sim->namemax * sizeof(const char *));
      if (tmp == NULL) {
        ERR(sim, MD_ERR_MEMALLOC, "MD_new_name call to realloc");
        return MD_FAIL;
      }
      sim->name = (const char **) tmp;
      sim->namemax *= 2;
    }
    sim->name[sim->namelen++] = name;
    sim->name[sim->namelen] = NULL;
  }
  return 0;
}


MD_Int MD_new_type(MD_Sim *sim, const char *name,
                   const MD_Int *subtype, MD_Int subtypelen, MD_Int size)
{
  MD_Int numbytes;

  /* check on incorrect use */
  if (sim->version != MD_API_VERSION) {
    ERR(sim, MD_ERR_VERSION, "MD_new_type");
    return MD_FAIL;
  }
  else if (sim->state != STATE_CREATE) {
    ERR(sim, MD_ERR_USE, "MD_new_type");
    return MD_FAIL;
  }

  /* check validity of length and name */
  if (subtypelen <= 0) {
    ERR(sim, MD_ERR_TYPE, "MD_new_type: illegal subtype array length");
    return MD_FAIL;
  }
  else if (ihash_lookup( &(sim->typetable), name) >= 0) {
    ERR(sim, MD_ERR_TYPE, "MD_new_type: name conflict");
    return MD_FAIL;
  }
  else if ((numbytes = validate_type(sim, subtype, subtypelen)) <= 0) {
    return MD_FAIL;
  }
  else if (sim->typelen == 0x0FFF || numbytes > 0xFFFF) {
    /* these range checks probably aren't necessary in practice */
    ERR(sim, MD_ERR_TYPE, "MD_new_type: range check failed");
    return MD_FAIL;
  }
  else if (numbytes != size) {
    ERR(sim, MD_ERR_TYPE, "MD_new_type: inconsistent size");
    return MD_FAIL;
  }
  return add_type(sim, name, subtype, subtypelen, numbytes);
}


MD_Int validate_type(MD_Sim *sim, const MD_Int *subtype, MD_Int len)
{
  MD_Int numbytes = 0;
  MD_Int stacksize = 0;
  MD_Int typenum, index, newstackmax, i;
  MD_Int lastflipbyte = 0;
  MD_Int alignment = 1;

  /*
   * check validity of type definition
   *
   * validity includes: previously defined types + correct word alignment
   *   assume that word alignment requires 4 byte quantities be aligned
   *   on every 4-byte word boundary; same with 8 byte quantities
   */

  /* push subtype array on stack */
  newstackmax = sim->typestackmax;
  while (len > newstackmax)  newstackmax *= 2;
  if (newstackmax > sim->typestackmax) {
    void *tmp = realloc(sim->typestack, newstackmax * sizeof(MD_Int));
    if (tmp == NULL) {
      ERR(sim, MD_ERR_MEMALLOC, "validate_type call to realloc");
      return MD_FAIL;
    }
    sim->typestack = (MD_Int *) tmp;
    sim->typestackmax = newstackmax;
  }
  for (i = len - 1;  i >= 0;  i--) {
    sim->typestack[stacksize++] = subtype[i];
  }
  sim->flipqueuelen = 0;
  while (stacksize > 0) {
    typenum = sim->typestack[--stacksize];
    index = INDEX(typenum);
    if (index < 0 || index >= sim->typelen
        || typenum != sim->type[index].typenum) {
      ERR(sim, MD_ERR_TYPE, "validate_type");
      return MD_FAIL;
    }
    switch (typenum) {
    case MD_INT:
    case MD_FLOAT:
    case MD_FVEC:
      if (numbytes % 4 != 0) {
        ERR(sim, MD_ERR_TYPE, "validate_type: illegal data alignment");
        return MD_FAIL;
      }
      if (sim->flipqueuelen == sim->flipqueuemax - 1) {
        void *tmp = realloc(sim->flipqueue, 2*sim->flipqueuemax*sizeof(MD_Int));
        if (tmp == NULL) {
          ERR(sim, MD_ERR_MEMALLOC, "validate_type call to realloc");
          return MD_FAIL;
        }
        sim->flipqueue = (MD_Int *) tmp;
        sim->flipqueuemax *= 2;
      }
      sim->flipqueue[sim->flipqueuelen] = numbytes - lastflipbyte;
      lastflipbyte = numbytes;
      numbytes += MD_SIZEOF(typenum);
      if (alignment < 4)  alignment = 4;
      if (typenum == MD_FVEC) {
        sim->flipqueue[sim->flipqueuelen++] |= MD_FLIP_12_BYTE;
      }
      else {
        sim->flipqueue[sim->flipqueuelen++] |= MD_FLIP_4_BYTE;
      }
      break;
    case MD_DOUBLE:
    case MD_DVEC:
      if (numbytes % 8 != 0) {
        ERR(sim, MD_ERR_TYPE, "validate_type: illegal data alignment");
        return MD_FAIL;
      }
      if (sim->flipqueuelen == sim->flipqueuemax - 1) {
        void *tmp = realloc(sim->flipqueue, 2*sim->flipqueuemax*sizeof(MD_Int));
        if (tmp == NULL) {
          ERR(sim, MD_ERR_MEMALLOC, "validate_type call to realloc");
          return MD_FAIL;
        }
        sim->flipqueue = (MD_Int *) tmp;
        sim->flipqueuemax *= 2;
      }
      sim->flipqueue[sim->flipqueuelen] = numbytes - lastflipbyte;
      lastflipbyte = numbytes;
      numbytes += MD_SIZEOF(typenum);
      if (alignment < 8)  alignment = 8;
      if (typenum == MD_DVEC) {
        sim->flipqueue[sim->flipqueuelen++] |= MD_FLIP_24_BYTE;
      }
      else {
        sim->flipqueue[sim->flipqueuelen++] |= MD_FLIP_8_BYTE;
      }
      break;
    case MD_CHAR:
    case MD_SHORT_STRING:
    case MD_STRING:
      numbytes += MD_SIZEOF(typenum);
      break;
    default:
      /* not a primitive type, must push its subtype array onto stack */
      newstackmax = sim->typestackmax;
      while (sim->type[index].subtypelen + stacksize > newstackmax) {
        newstackmax *= 2;
      }
      if (newstackmax > sim->typestackmax) {
        void *tmp = realloc(sim->typestack, newstackmax * sizeof(MD_Int));
        if (tmp == NULL) {
          ERR(sim, MD_ERR_MEMALLOC, "validate_type call to realloc");
          return MD_FAIL;
        }
        sim->typestack = (MD_Int *) tmp;
        sim->typestackmax = newstackmax;
      }
      for (i = sim->type[index].subtypelen - 1;  i >= 0;  i--) {
        sim->typestack[stacksize++] = sim->type[index].subtype[i];
      }
    }
  }
  if (sim->flipqueuelen == 0) {
    /* nothing to flip in type */
    sim->flipflags = MD_FLIP_NONE;
  }
  else if (len == 1) {
    /* single type that requires flipping */
    sim->flipflags = (MD_FLIP_MASK & sim->flipqueue[sim->flipqueuelen - 1]);
  }
  else {
    /* composite type that requires flipping */
    sim->flipflags = MD_FLIP_COMPOSITE;
  }
  /* left at least 1 element at end of queue */
  sim->flipqueue[sim->flipqueuelen] = numbytes - lastflipbyte;
  if (numbytes % alignment != 0) {
    ERR(sim, MD_ERR_TYPE, "validate_type: illegal data alignment");
    return MD_FAIL;
  }
  return numbytes;
}


/*
 * must first call validate_type()
 */
MD_Int add_type(MD_Sim *sim, const char *name,
                const MD_Int *subtype, MD_Int len, MD_Int numbytes)
{
  MD_Int typenum;

  /* put all type information into type array */
  if (sim->typelen == sim->typemax) {
    void *tmp = realloc(sim->type, 2 * sim->typemax * sizeof(MD_Type));
    if (tmp == NULL) {
      ERR(sim, MD_ERR_MEMALLOC, "add_type call to realloc");
      return MD_FAIL;
    }
    sim->type = (MD_Type *) tmp;
    sim->typemax *= 2;
  }
  memset( &(sim->type[sim->typelen]), 0, sizeof(MD_Type));
  sim->type[sim->typelen].typenum = typenum 
      = (sim->typelen << 19) | sim->flipflags | numbytes;
  sim->type[sim->typelen].name = name;

  /* store the typelist */
  if ((sim->type[sim->typelen].subtype
        = (MD_Int *) malloc(len * sizeof(MD_Int))) == NULL) {
    ERR(sim, MD_ERR_MEMALLOC, "add_type call to malloc");
    return MD_FAIL;
  }
  memcpy(sim->type[sim->typelen].subtype, subtype, len * sizeof(MD_Int));
  sim->type[sim->typelen].subtypelen = len;

  /* store the flip offsets if this is a flip-able composite type */
  /* this information is determined by validate_type() */
  if (sim->flipflags == MD_FLIP_COMPOSITE) {
    if ((sim->type[sim->typelen].flipoffset
          = (MD_Int *) malloc((sim->flipqueuelen+1)*sizeof(MD_Int))) == NULL) {
      ERR(sim, MD_ERR_MEMALLOC, "add_type call to malloc");
      return MD_FAIL;
    }
    memcpy(sim->type[sim->typelen].flipoffset, sim->flipqueue, 
           (sim->flipqueuelen + 1) * sizeof(MD_Int));
    sim->type[sim->typelen].flipoffsetlen = sim->flipqueuelen;
  }
  sim->typelen++;

  /* insert type name into hash table */
  if (ihash_insert( &(sim->typetable), name, typenum) != IHASH_FAIL) {
    ERR(sim, MD_ERR_MEMALLOC, "add_type call to ihash_insert");
    return MD_FAIL;
  }

  /* put name into type namelist */
  if (sim->typenamelen == sim->typenamemax - 1) {
    void *tmp = realloc(sim->typename, 2*sim->typenamemax*sizeof(const char *));
    if (tmp == NULL) {
      ERR(sim, MD_ERR_MEMALLOC, "add_type call to realloc");
      return MD_FAIL;
    }
    sim->typename = (const char **) tmp;
    sim->typenamemax *= 2;
  }
  sim->typename[sim->typenamelen++] = name;
  sim->typename[sim->typenamelen] = NULL;
  return typenum;
}


MD_Int MD_new_err(MD_Sim *sim, const char *errmsg, MD_Int isfatal)
{
  /* check on incorrect use */
  if (sim->version != MD_API_VERSION) {
    ERR(sim, MD_ERR_VERSION, "MD_new_err");
    return MD_FAIL;
  }
  else if (sim->state != STATE_CREATE) {
    ERR(sim, MD_ERR_USE, "MD_new_err");
    return MD_FAIL;
  }

  if (sim->errlen == sim->errmax) {
    void *tmp = realloc(sim->err, 2 * sim->errmax * sizeof(MD_Err));
    if (tmp == NULL) {
      ERR(sim, MD_ERR_MEMALLOC, "MD_new_err call to realloc");
      return MD_FAIL;
    }
    sim->err = (MD_Err *) tmp;
    sim->errmax *= 2;
  }
  sim->err[sim->errlen].msg = errmsg;
  sim->err[sim->errlen].isfatal = -(isfatal != 0);  /* nonzero ==> -1 */
  return sim->errlen++;
}


MD_Errcode MD_set_run(MD_Sim *sim, MD_Errcode (*run)(MD_Sim *, MD_Int))
{
  if (sim->version != MD_API_VERSION) {
    ERR(sim, MD_ERR_VERSION, "MD_set_run");
    return MD_FAIL;
  }
  else if (sim->iscreated) {
    ERR(sim, MD_ERR_USE, "MD_set_run");
    return MD_FAIL;
  }
  sim->run = run;
  return 0;
}


MD_Errcode MD_set_init(MD_Sim *sim, MD_Errcode (*init)(MD_Sim *))
{
  if (sim->version != MD_API_VERSION) {
    ERR(sim, MD_ERR_VERSION, "MD_set_init");
    return MD_FAIL;
  }
  else if (sim->iscreated) {
    ERR(sim, MD_ERR_USE, "MD_set_init");
    return MD_FAIL;
  }
  sim->init = init;
  return 0;
}


/*****************************************************************************
 *
 * Call this either during create_engine or init.
 *
 *****************************************************************************/

MD_Errcode MD_set_stepnum(MD_Sim *sim, MD_Int stepnum)
{
  if (sim->version != MD_API_VERSION) {
    ERR(sim, MD_ERR_VERSION, "MD_set_stepnum");
    return MD_FAIL;
  }
  else if (sim->state != STATE_CREATE && sim->state != STATE_INIT) {
    ERR(sim, MD_ERR_USE, "MD_set_stepnum");
    return MD_FAIL;
  }
  sim->stepnum = stepnum;
  return 0;
}


/*****************************************************************************
 *
 * Default setlen, setdata, and getdata routines.
 *
 *****************************************************************************/

MD_Errcode MD_def_setlen(MD_Sim *sim, MD_Data *data, MD_Int len)
{
  if (*(data->plen) != len) {
    MD_Int i, idnum;
    void *tmp;

    /* need to resize all arrays in my circular link */
    i = idnum = data->attr.idnum;
    while (1) {
      tmp = realloc(*((void **) (sim->data[i].ppdata)),
                    len * MD_SIZEOF(sim->data[i].attr.typenum));
      if (tmp == NULL) {
        ERR(sim, MD_ERR_MEMALLOC, "MD_def_setlen call to realloc");
        return MD_FAIL;
      }
      *((void **) (sim->data[i].ppdata)) = tmp;
      *(sim->data[i].plen) = len;
      if (sim->data[i].link == idnum)  break;
      i = sim->data[i].link;
    }
  }
  return 0;
}


MD_Errcode MD_def_setdata(MD_Sim *sim, MD_Data *data, const void *buf,
                          MD_Int nelems, MD_Int first)
{
  void *dest = *((char **) (data->ppdata))
               + (first * MD_SIZEOF(data->attr.typenum));

  /* assume that memory locations are either identical or they don't overlap */
  if (nelems > 0 && buf != dest) {
    memcpy(dest, buf, nelems * MD_SIZEOF(data->attr.typenum));
  }
  return 0;
}


MD_Errcode MD_def_getdata(MD_Sim *sim, MD_Data *data, void *buf,
                          MD_Int nelems, MD_Int first)
{
  void *src = *((char **) (data->ppdata))
              + (first * MD_SIZEOF(data->attr.typenum));

  /* assume that memory locations are either identical or they don't overlap */
  if (nelems > 0 && buf != src) {
    memcpy(buf, src, nelems * MD_SIZEOF(data->attr.typenum));
  }
  return 0;
}


/*****************************************************************************
 *
 * Obtain the engine handle.
 *
 *****************************************************************************/

void *MD_engine(MD_Sim *sim)
{
  return sim->engine;
}


/*****************************************************************************
 *
 * Call back at end of each step and before each force evaluation.
 *
 *****************************************************************************/

MD_Errcode MD_callback(MD_Sim *sim, MD_Int stepnum)
{
  if (sim->state != STATE_RUN && sim->state != STATE_INIT) {
    ERR(sim, MD_ERR_USE, "MD_callback");
    return MD_FAIL;
  }
  else if (stepnum != sim->stepnum) {
    ERR(sim, MD_ERR_CALLBACK, "MD_callback, illegal stepnum");
    return MD_FAIL;
  }
  sim->prev_state = sim->state;  /* save current state */
  sim->state = STATE_CALLBACK;   /* set call state */
  sim->stepnum = stepnum;
  sim->numer = 0;
  sim->denom = 0;
  sim->cbptr = sim->cb;
  sim->cbcounter = 0;
  sim->cbtotal = sim->cblen;
  sim->cbnumcalls = 0;
  sim->iscb = 0;
  return 0;
}


MD_Errcode MD_fcallback(MD_Sim *sim, MD_Int stepnum,
                        MD_Int numer, MD_Int denom)
{
  if (sim->state != STATE_RUN && sim->state != STATE_INIT) {
    ERR(sim, MD_ERR_USE, "MD_fcallback");
    return MD_FAIL;
  }
  else if (stepnum != sim->stepnum && stepnum != sim->stepnum + 1) {
/*
 * should we also test numer and denom?
 */
    ERR(sim, MD_ERR_CALLBACK, "MD_fcallback, illegal stepnum");
    return MD_FAIL;
  }
  sim->prev_state = sim->state;  /* save current state */
  sim->state = STATE_FCALLBACK;  /* set call state */
  sim->stepnum = stepnum;
  sim->numer = numer;
  sim->denom = denom;
  sim->cbptr = sim->fcb;
  sim->cbcounter = 0;
  sim->cbtotal = sim->fcblen;
  sim->cbnumcalls = 0;
  sim->iscb = 0;
  return 0;
}


MD_Int MD_isdonecb(MD_Sim *sim)
{
  MD_Callback_Data *pcbdata;
  MD_Int datanum, i, idnum, len;
  void *ppdata;

  switch (sim->state) {

  case STATE_CALLBACK:
  case STATE_FCALLBACK:
    if (sim->errnum != 0)  return 1;
    else if (sim->iscb && sim->cbnumcalls > sim->cbmaxcalls) {
      ERR(sim, MD_ERR_CALLBACK,
          "MD_isdone, exceeded maximum consecutive calls");
      return 1;
    }
    else if (sim->iscb) {
      if (sim->isdoneprcb(sim)) {
        if (sim->waitprcb(sim) != 0) {
          /* make sure errnum is set */
          if (sim->errnum == 0) {
            ERR(sim, MD_ERR_CALLBACK, "MD_isdonecb call to waitprcb");
          }
          else {
            sim->err[sim->errnum].isfatal = MD_FAIL;
          }
          return MD_FAIL;
        }
        sim->iscb = 0;
        sim->cbnumcalls = 0;
        sim->cbcounter++;
        sim->cbptr++;
      }
      else {
        sim->cbnumcalls++;
      }
    }

    if (!sim->iscb) {
      /* done with previous callback, find next one */
      while (sim->cbcounter < sim->cbtotal
             && (sim->stepnum % sim->cbptr->stepnum != 0
                 || ((sim->cbptr->numer != 0 || sim->cbptr->denom != 0)
                     && (sim->cbptr->numer != sim->numer
                         || sim->cbptr->denom != sim->denom)))) {
        sim->cbcounter++;
        sim->cbptr++;
      }
      if (sim->cbcounter == sim->cbtotal)  return 1;

      /* save the callback information */
      sim->cbsave = *(sim->cbptr);

      /* make ranges for callback data arrays valid */
      /* set pointers to engine data */
      pcbdata = sim->cbsave.data;
      datanum = sim->cbsave.datanum;
      for (i = 0;  i < datanum;  i++, pcbdata++) {
        idnum = pcbdata->attr.idnum;
        len = *(sim->data[idnum].plen);
        ppdata = sim->data[idnum].ppdata;
        if (ppdata != NULL && pcbdata->first < len) {
          if (pcbdata->nelems > len - pcbdata->first) {
            pcbdata->nelems = len - pcbdata->first;
          }
          /* set callback data pointer to data array */
          pcbdata->data = *((char **) ppdata)
                          + (pcbdata->first * MD_SIZEOF(pcbdata->attr.typenum));
        }
        else {
          pcbdata->data = NULL;
          pcbdata->first = len;
          pcbdata->nelems = 0;
        }
      }

      sim->cbsave.stepnum = sim->stepnum;
      sim->cbsave.numer = sim->numer;
      sim->cbsave.denom = sim->denom;
      if (sim->prcb(sim, &(sim->cbsave)) != 0) {
        /* make sure errnum is set */
        if (sim->errnum == 0) {
          ERR(sim, MD_ERR_CALLBACK, "MD_isdonecb call to prcb");
        }
        else {
          sim->err[sim->errnum].isfatal = MD_FAIL;
        }
        return 1;
      }
      sim->iscb = 1;
      sim->cbnumcalls++;
    }
    break;

  default:
    ERR(sim, MD_ERR_USE, "MD_isdonecb");
    return 1;
  }
  return 0;
}


MD_Errcode MD_waitcb(MD_Sim *sim)
{
  MD_Callback_Data *pcbdata;
  MD_Int datanum, i, idnum, len;
  void *ppdata;
  MD_Int notdone = 1;

  switch (sim->state) {

  case STATE_CALLBACK:
  case STATE_FCALLBACK:
    sim->state = sim->prev_state;  /* restore previous state */
    if (sim->errnum != 0) {
      /* make sure error is fatal */
      sim->err[sim->errnum].isfatal = MD_FAIL;
      return MD_FAIL;
    }
    else if (sim->iscb && sim->cbnumcalls > sim->cbmaxcalls) {
      ERR(sim, MD_ERR_CALLBACK,
          "MD_waitcb, exceeded maximum consecutive calls");
      return MD_FAIL;
    }
    else if (sim->iscb) {
      while (!sim->isdoneprcb(sim)) {
        sim->cbnumcalls++;
        if (sim->cbnumcalls > sim->cbmaxcalls) {
          ERR(sim, MD_ERR_CALLBACK,
              "MD_waitcb, exceeded maximum consecutive calls");
          return MD_FAIL;
        }
      }
      sim->iscb = 0;
      sim->cbnumcalls = 0;
      if (sim->waitprcb(sim) != 0) {
        /* make sure errnum is set */
        if (sim->errnum == 0) {
          ERR(sim, MD_ERR_CALLBACK, "MD_waitcb call to waitprcb");
        }
        else {
          sim->err[sim->errnum].isfatal = MD_FAIL;
        }
        return MD_FAIL;
      }
      sim->cbcounter++;
      sim->cbptr++;
    }

    while (notdone) {
      /* process all remaining callbacks */
      /* done with previous callback, find next one */
      while (sim->cbcounter < sim->cbtotal
             && (sim->stepnum % sim->cbptr->stepnum != 0
                 || ((sim->cbptr->numer != 0 || sim->cbptr->denom != 0)
                     && (sim->cbptr->numer != sim->numer
                         || sim->cbptr->denom != sim->denom)))) {
        sim->cbcounter++;
        sim->cbptr++;
      }
      if (sim->cbcounter == sim->cbtotal)  notdone = 0;
      else {
        /* save the callback information */
        sim->cbsave = *(sim->cbptr);

        /* make ranges for callback data arrays valid */
        /* set pointers to engine data */
        pcbdata = sim->cbsave.data;
        datanum = sim->cbsave.datanum;
        for (i = 0;  i < datanum;  i++, pcbdata++) {
          idnum = pcbdata->attr.idnum;
          len = *(sim->data[idnum].plen);
          ppdata = sim->data[idnum].ppdata;
          if (ppdata != NULL && pcbdata->first < len) {
            if (pcbdata->nelems > len - pcbdata->first) {
              pcbdata->nelems = len - pcbdata->first;
            }
            /* set callback data pointer to data array */
            pcbdata->data = *((char **) ppdata)
                        + (pcbdata->first * MD_SIZEOF(pcbdata->attr.typenum));
          }
          else {
            pcbdata->data = NULL;
            pcbdata->first = len;
            pcbdata->nelems = 0;
          }
        }

        sim->cbsave.stepnum = sim->stepnum;
        sim->cbsave.numer = sim->numer;
        sim->cbsave.denom = sim->denom;
        if (sim->prcb(sim, &(sim->cbsave)) != 0) {
          /* make sure errnum is set */
          if (sim->errnum == 0) {
            ERR(sim, MD_ERR_CALLBACK, "MD_waitcb call to prcb");
          }
          else {
            sim->err[sim->errnum].isfatal = MD_FAIL;
          }
          return MD_FAIL;
        }
        sim->cbnumcalls++;
        while (sim->cbnumcalls <= sim->cbmaxcalls
               && !sim->isdoneprcb(sim)) {
          sim->cbnumcalls++;
        }
        if (sim->cbnumcalls > sim->cbmaxcalls) {
          ERR(sim, MD_ERR_CALLBACK,
              "MD_waitcb, exceeded maximum consecutive calls");
          return MD_FAIL;
        }
        else if (sim->errnum != 0) {
          /* make sure error is fatal */
          sim->err[sim->errnum].isfatal = MD_FAIL;
          return MD_FAIL;
        }
        else if (sim->waitprcb(sim) != 0) {
          /* make sure errnum is set */
          if (sim->errnum == 0) {
            ERR(sim, MD_ERR_CALLBACK, "MD_waitcb call to waitprcb");
          }
          else {
            sim->err[sim->errnum].isfatal = MD_FAIL;
          }
          return MD_FAIL;
        }
        sim->cbnumcalls = 0;
        sim->cbcounter++;
        sim->cbptr++;
      }
    }
    break;

  default:
    ERR(sim, MD_ERR_USE, "MD_waitcb");
    return 1;
  }
  return 0;
}


/*****************************************************************************
 *
 * Raise an error condition.
 *
 *****************************************************************************/

void MD_raiserr(MD_Sim *sim, const char *myname, MD_Int errnum,
                const char *optmsg, const char *filename, MD_Int linenum)
{
  char format_buf[511];
  char empty_str[1] = "";
  const char *msg;

  if (errnum < 0 || errnum >= sim->errlen)  return;
  sim->errnum = errnum;

  /* space for myname */
  /* if a simulation is currently running, print the step number */
  if (myname != NULL) {
    if (sim->state == STATE_RUN
        || sim->state == STATE_CALLBACK || sim->state == STATE_FCALLBACK) {
      if (sim->numer != 0 && sim->denom != 0) {
        sprintf(format_buf, "%%.%ds, step number %d %d/%d, ",
                MD_STRBUF_SIZE / 4, sim->stepnum, sim->numer, sim->denom);
      }
      else {
        sprintf(format_buf, "%%.%ds, step number %d, ", MD_STRBUF_SIZE / 4,
                sim->stepnum);
      }
    }
    else {
      sprintf(format_buf, "%%.%ds ", MD_STRBUF_SIZE / 4);
    }
  }
  else {
    myname = empty_str;
    if (sim->state == STATE_RUN
        || sim->state == STATE_CALLBACK || sim->state == STATE_FCALLBACK) {
      if (sim->numer != 0 && sim->denom != 0) {
        sprintf(format_buf, "%%sstep number %d %d/%d, ",
                sim->stepnum, sim->numer, sim->denom);
      }
      else {
        sprintf(format_buf, "%%sstep number %d, ", sim->stepnum);
      }
    }
    else {
      sprintf(format_buf, "%%s");
    }
  }

  /* is error fatal? */
  if (sim->err[errnum].isfatal != 0) {
    sprintf(format_buf + strlen(format_buf), "fatal error %d: ", errnum);
  }
  else {
    sprintf(format_buf + strlen(format_buf), "error %d: ", errnum);
  }

  /* space for msg */
  msg = sim->err[errnum].msg;
  if (msg != NULL) {
    sprintf(format_buf + strlen(format_buf), "%%.%ds: ", MD_STRBUF_SIZE / 8);
  }
  else {
    msg = empty_str;
    sprintf(format_buf + strlen(format_buf), "%%s");
  }

  /* space for optmsg */
  if (optmsg != NULL) {
    sprintf(format_buf + strlen(format_buf), "%%.%ds, ",
            MD_STRBUF_SIZE / 4);
  }
  else {
    optmsg = empty_str;
    sprintf(format_buf + strlen(format_buf), "%%s");
  }

  /* space for filename */
  if (filename != NULL) {
    sprintf(format_buf + strlen(format_buf), "file %%.%ds, ",
            MD_STRBUF_SIZE / 8);
  }
  else {
    filename = empty_str;
    sprintf(format_buf + strlen(format_buf), "%%s");
  }

  /* print linenum into format string */
  if (linenum > 0) {
    sprintf(format_buf + strlen(format_buf), "line %d", linenum);
  }

  /* use format string to print into error message buffer */
  sprintf(sim->errmsg_buf, format_buf, myname, msg, optmsg, filename);
}