/**
 * This is a new version of the four circle reflection measurement module. Its task is 
 * to provide functionality required for measuring a list of reflections. This is done 
 * with the help of some scripts binding things together. This module provides:
 * - a table of scan parameters for measuring reflections at different angles. 
 * - the functions required to produce the scan data files for single detector 
 *   mode.
 * - a list of reflections to measure. 
 * 
 * In contrast to the old mesure module, this only provides some help to 
 * scripts. Thus, through scripts, a more flexible measurement is possible. 
 * 
 * copyright: see file COPYRIGHT
 * 
 * Mark Koennecke, July 2008
 */
#include <stdio.h>
#include <assert.h>
#include "singlex.h"
#include "sicsobj.h"
#include "fourtable.h"
#include "sicshipadaba.h"
#include "sicsvar.h"
#include <sics.h>
#include "scan.h"
#include "stdscan.h"
#include "evcontroller.h"
#include "integrate.h"
#include "sginfo.h"
#include "matrix/matrix.h"
#include "cell.h"
#include "fourlib.h"
#include "counter.h"
#include "scanvar.h"
#include "scan.i"
#include "sdynar.h"
#include "lld.h"

extern char *trim(char *);

extern void SNXFormatTime(char *pBueffel, int iLen);

#define ABS(x) (x < 0 ? -(x) : (x))
/*---------------------------------------------------------------------------*/
typedef struct {
  FILE *profFile;               /* file with reflection profiles, ccl */
  FILE *hklFile;                /* file with integrated intensities */
  pSICSOBJ stepTable;                /* table with the scan parameters */
  char *currentFileRoot;
  pSICSOBJ messList;
  pHdb currentRefl;             /* the current reflection being measured */
  int count;
  pScanData pScanner;
  int masterCount;              /* the number of master reflection as craeted by indgen */
} FourMess, *pFourMess;
/*---------------------------------------------------------------------------*/
static void KillFourMess(void *data)
{
  pFourMess priv = (pFourMess) data;
  if (priv == NULL) {
    return;
  }
  if (priv->profFile != NULL) {
    fclose(priv->profFile);
  }
  if (priv->hklFile != NULL) {
    fclose(priv->hklFile);
  }
  if (priv->currentFileRoot != NULL) {
    free(priv->currentFileRoot);
  }
  if (priv->stepTable >= 0) {
    DeleteFourCircleTable(priv->stepTable);
  }
  free(priv);
}

/*----------------------------------------------------------------------------*/
static int FourMessAction(SConnection * pCon, SicsInterp * pSics,
                          void *data, int argc, char *argv[])
{
  pFourMess priv = NULL;
  pSICSOBJ self = (pSICSOBJ) data;
  int err, status;

  assert(self != NULL);
  priv = self->pPrivate;
  assert(priv != NULL);

  if (argc < 2) {
    SCWrite(pCon, "ERROR: insufficient number of arguments to fmess",
            eError);
    return 0;
  }

  if (strcmp(argv[1], "table") == 0) {
    return HandleFourCircleCommands(priv->stepTable, pCon,
                                    argc, argv, &err);
  }

  status = InvokeSICSOBJ(pCon, pSics, data, argc, argv);
  if (status < 0) {
    SCPrintf(pCon, eError, "ERROR: %s no subcommand or parameter",
             argv[1]);
    return 0;
  }
  return status;
}

/*-----------------------------------------------------------------------------*/
static int FourMessClose(pSICSOBJ self, SConnection * pCon,
                         pHdb commandNode, pHdb par[], int nPar)
{
  pFourMess priv = NULL;

  priv = self->pPrivate;
  if (priv->hklFile != NULL) {
    fclose(priv->hklFile);
    priv->hklFile = NULL;
  }
  if (priv->profFile != NULL) {
    fclose(priv->profFile);
    priv->profFile = NULL;
  }
  return 1;
}

/*-----------------------------------------------------------------------------*/
static int FourMessStart(pSICSOBJ self, SConnection * pCon,
                         pHdb commandNode, pHdb par[], int nPar)
{
  pFourMess priv = NULL;
  char pFilename[1024], pRoot[512];
  char pBueffel[1024];
  double lambda;
  const double *dUB;
  pSicsVariable pVar = NULL;
  char *pFile = NULL, *pPtr;
  FILE *temp = NULL;
  pHdb node;
  hdbValue v;

  assert(self);
  assert(pCon);

  priv = self->pPrivate;

  if (nPar < 1) {
    SCWrite(pCon, "ERROR: need file name parameter to start", eError);
    return 0;
  }

  /* close open files if so */
  if (priv->hklFile != NULL) {
    FourMessClose(self, pCon, commandNode, par, nPar);
  }

  /* create filename root */
  GetHipadabaPar(par[0], &v, pCon);
  pFile = strdup(v.v.text);
  pPtr = strrchr(pFile, (int) '.');
  pPtr++;
  *pPtr = '\0';
  if (priv->currentFileRoot != NULL) {
    free(priv->currentFileRoot);
  }
  priv->currentFileRoot = strdup(pFile);
  free(pFile);
  strlcpy(pRoot, priv->currentFileRoot, 511);


  /* open the reflection file */
  sprintf(pBueffel, "Writing to %s.ccl, .rfl", pRoot);
  SCWrite(pCon, pBueffel, eLog);
  strcpy(pFilename, pRoot);
  strcat(pFilename, "ccl");
  priv->profFile = fopen(pFilename, "w");
  if (!priv->profFile) {
    sprintf(pBueffel, "ERROR: SERIOUS TROUBLE: cannot open %s!",
            pFilename);
    SCWrite(pCon, pBueffel, eError);
    SCSetInterrupt(pCon, eAbortBatch);
    return 0;
  }

  node = GetHipadabaNode(self->objectNode, "template");
  assert(node != NULL);
  GetHipadabaPar(node, &v, pCon);
  temp = fopen(v.v.text, "r");
  if (temp == NULL) {
    SCWrite(pCon, "ERROR: failed to open header template", eError);
  }
  if (temp != NULL) {
    WriteTemplate(priv->profFile, temp, pFilename, NULL,
                  pCon, pServ->pSics);
    fclose(temp);
  }

  /* open hkl-data file */
  strcpy(pFilename, pRoot);
  strcat(pFilename, "rfl");
  priv->hklFile = fopen(pFilename, "w");
  if (!priv->hklFile) {
    sprintf(pBueffel, "ERROR: SERIOUS TROUBLE: cannot open %s!",
            pFilename);
    SCWrite(pCon, pBueffel, eError);
    SCSetInterrupt(pCon, eAbortBatch);
    return 0;
  }
  fputs(pFilename, priv->hklFile);
  fputs("\n", priv->hklFile);

  /* write some header data */
  SNXFormatTime(pBueffel, 1024);
  fprintf(priv->hklFile, "filetime = %s\n", pBueffel);
  lambda = SXGetLambda();
  fprintf(priv->hklFile, "lambda = %f Angstroem\n", lambda);
  dUB = SXGetUB();
  fprintf(priv->hklFile,
          "UB = %7.6f %7.6f %7.6f %7.6f %7.6f %7.6f %7.6f %7.6f %7.6f\n",
          dUB[0], dUB[1], dUB[2], dUB[3], dUB[4], dUB[5], dUB[6], dUB[7],
          dUB[8]);

  /* write sample & user info */
  strcpy(pBueffel, "CCL, Instr=TRICS, ");
  pVar = FindVariable(pServ->pSics, "sample");
  if (pVar) {
    fprintf(priv->hklFile, "sample = %s\n", pVar->text);
  }
  pVar = FindVariable(pServ->pSics, "user");
  if (pVar) {
    fprintf(priv->hklFile, "user = %s \n", pVar->text);
  }
  priv->count = 0;
  return 1;
}
/*-----------------------------------------------------------------------------*/
static int FourMessReopen(pSICSOBJ self, SConnection * pCon,
                         pHdb commandNode, pHdb par[], int nPar)
{
	pFourMess priv = NULL;
	char pFilename[1024];
	
	priv = self->pPrivate;
	if(priv->hklFile != NULL){
		return 1; /* still open! */
	}
	
	if(priv->currentFileRoot == NULL){
		SCWrite(pCon,"ERROR: nothing to append too", eLogError);
		return 0;
	}
	
	strcpy(pFilename, priv->currentFileRoot);
	strcat(pFilename, "ccl");
	priv->profFile = fopen(pFilename, "a");
	if(priv->profFile == NULL){
		SCPrintf(pCon,eLogError,"ERROR: failed to reopen %s", pFilename);
		return 0;
	}
	strcpy(pFilename, priv->currentFileRoot);
	strcat(pFilename, "rfl");
	priv->hklFile = fopen(pFilename, "a");
	if(priv->hklFile == NULL){
		SCPrintf(pCon,eLogError,"ERROR: failed to reopen %s", pFilename);
		return 0;
	}
	SCPrintf(pCon,eValue,"Reopened %s.ccl,.rfl", priv->currentFileRoot);
	
	return 1;
}
/*----------------------------------------------------------------------------*/
static int FourMessScanPar(pSICSOBJ self, SConnection * pCon,
                           pHdb commandNode, pHdb par[], int nPar)
{
  pFourMess priv = self->pPrivate;
  char *scanvar;
  double stt, dVal;
  int np, preset;
  pDynString data;

  assert(priv != NULL);
  if (nPar < 1) {
    SCWrite(pCon, "ERROR: need two theta parameter to scanpar", eError);
    return 0;
  }
  stt = par[0]->value.v.doubleValue;

  scanvar = GetFourCircleScanVar(priv->stepTable, stt);
  dVal = GetFourCircleStep(priv->stepTable, stt);
  np = GetFourCircleScanNP(priv->stepTable, stt);
  preset = GetFourCirclePreset(priv->stepTable, stt);
  if (strcmp(scanvar, "Not found") == 0) {
    SCPrintf(pCon, eValue, "%s,%f,%d,%d", "om", dVal, np, preset);
  } else {
    SCPrintf(pCon, eValue, "%s,%f,%d,%d", scanvar, dVal, np, preset);
  }
  return 1;
}
/*--------------------------------------------------------------------------*/
static int FourMessPrepareScan(pScanData self)
{
  pVarEntry pVar = NULL;
  void *pDings;
  int i, iRet;
  float fVal;
  char pBueffel[512];
  char pMessage[1024];

  assert(self);
  assert(self->pCon);

  /* check boundaries of scan variables and allocate storage */
  for (i = 0; i < self->iScanVar; i++) {
    DynarGet(self->pScanVar, i, &pDings);
    pVar = (pVarEntry) pDings;
    if (pVar) {
      iRet = CheckScanVar(pVar, self->pCon, self->iNP - 1);
      if (!iRet) {
        return 0;
      }
      InitScanVar(pVar);
    } else {
      SCWrite(self->pCon,
              "WARNING: Internal error, no scan variable, I try to continue",
              eLog);
    }
    pVar = NULL;
  }                             /* end for */

  /* configure counter */
  SetCounterMode((pCounter) self->pCounterData, self->iMode);
  SetCounterPreset((pCounter) self->pCounterData, self->fPreset);
  self->iCounts = 0;

  return 1;
}
/*-----------------------------------------------------------------------------*/
static int FourMessPrepareCmd(pSICSOBJ self, SConnection * pCon,
                         pHdb commandNode, pHdb par[], int nPar)
{
	  pFourMess priv = self->pPrivate;

	  return FourMessPrepareScan(priv->pScanner); 
}
/*----------------------------------------------------------------------------*/
static hdbCallbackReturn SetScannerCB(pHdb node, void *userData,
                                      pHdbMessage mm)
{
  pHdbDataMessage set = NULL;
  pFourMess priv = (pFourMess) userData;
  SConnection *pCon = NULL;
  pScanData old;

  if ((set = GetHdbSetMessage(mm)) != NULL) {
    old = priv->pScanner;
    priv->pScanner =
        FindCommandData(pServ->pSics, set->v->v.text, "ScanObject");
    if (priv->pScanner == NULL) {
      priv->pScanner = old;
      pCon = set->callData;
      if (pCon != NULL) {
        SCWrite(pCon, "ERROR: scan object not found", eError);
      }
      return hdbAbort;
    }
  }
  return hdbContinue;
}

/*---------------------------------------------------------------------------*/
static double getProtonAverage(pFourMess self)
{
  int np, i;
  long *lData = NULL, lSum = 0;

  np = GetScanNP(self->pScanner);
  lData = (long *) malloc((np + 1) * sizeof(long));
  if (lData == NULL || np == 0) {
    return 0.;
  }
  memset(lData, 0, (np + 1) * sizeof(long));
  GetScanMonitor(self->pScanner, 2, lData, np);
  for (i = 0; i < np; i++) {
    lSum += lData[i];
  }
  free(lData);
  return (double) lSum / (double) np;
}

/*---------------------------------------------------------------------------*/
static int FourMessStoreIntern(pSICSOBJ self, SConnection * pCon,
                               double fHkl[3], double fPosition[4], char *extra)
{
  pFourMess priv = self->pPrivate;
  float fSum, fSigma, fTemp, fStep = .0, fPreset =.0, fMF;
  int i, iLF, iRet, iNP, ii;
  long *lCounts = NULL;
  pEVControl pEva = NULL;
  pDummy pPtr = NULL;
  pIDrivable pDriv = NULL;
  char pBueffel[512], pTime[512], pNum[10];
  double prot;

  if (priv->pScanner == NULL) {
    SCWrite(pCon, "ERROR: store: scan not configured", eLogError);
    return 0;
  }

  if (priv->hklFile == NULL) {
    SCWrite(pCon, "ERROR: store: no files open", eLogError);
    return 0;
  }
  priv->count++;

  /* get necessary data */
  fSum = 0.;
  fSigma = 0.;
  iRet = ScanIntegrate(priv->pScanner, &fSum, &fSigma);
  if (iRet != 1) {
    switch (iRet) {
    case INTEGLEFT:
      sprintf(pBueffel,
              "WARNING: integration failed --> no left side to: %f %f %f",
              fHkl[0], fHkl[1], fHkl[2]);
      break;
    case INTEGRIGHT:
      sprintf(pBueffel,
              "WARNING: integration failed -->no right side to: %f %f %f",
              fHkl[0], fHkl[1], fHkl[2]);
      break;
    case INTEGNOPEAK:
      sprintf(pBueffel,
              "WARNING: integration failed -->no peak found: %f %f %f",
              fHkl[0], fHkl[1], fHkl[2]);
      break;
    case INTEGFUNNYBACK:
      sprintf(pBueffel,
              "WARNING: integration problem, asymmetric background: %f %f %f",
              fHkl[0], fHkl[1], fHkl[2]);
      break;
    }
    SCWrite(pCon, pBueffel, eLog);
  }
  iNP = GetScanNP(priv->pScanner);
  lCounts = malloc(iNP * sizeof(long));
  if (lCounts == NULL) {
    SCWrite(pCon, "ERROR: out of memory in store", eLogError);
    return 0;
  }
  GetScanCounts(priv->pScanner, lCounts, iNP);

  /* write it */
  if (priv->profFile) {
    fprintf(priv->profFile,
            "%4d %7.3f %7.3f %7.3f %7.2f %7.2f %7.2f %7.2f %7.0f %7.2f\n",
            priv->count, fHkl[0], fHkl[1], fHkl[2], fPosition[0],
            fPosition[1], fPosition[2], fPosition[3], fSum, fSigma);
  }
  if (priv->hklFile) {
    fprintf(priv->hklFile,
            "%5d %6.2f %6.2f %6.2f %7.2f %7.2f %7.2f %7.2f %7.0f %7.2f\n",
            priv->count, fHkl[0], fHkl[1], fHkl[2], fPosition[0],
            fPosition[1], fPosition[2], fPosition[3], fSum, fSigma);
  }
  sprintf(pBueffel,
          "%5d %6.2f %6.2f %6.2f %7.2f %7.2f %7.2f %7.2f %7.0f %7.2f\n",
          priv->count, fHkl[0], fHkl[1], fHkl[2], fPosition[0],
          fPosition[1], fPosition[2], fPosition[3], fSum, fSigma);
  SCWrite(pCon, pBueffel, eLog);

  /* get temperature */
  fTemp = -777.77;
  pEva = (pEVControl) FindCommandData(pServ->pSics, "temperature",
                                      "Environment Controller");
  if (pEva == NULL) {
    pPtr = (pDummy) FindCommandData(pServ->pSics, "temperature",
                                    "RemObject");
    if (pPtr != NULL) {
      pDriv = pPtr->pDescriptor->GetInterface(pPtr, DRIVEID);
      if (pDriv != NULL) {
        fTemp = pDriv->GetValue(pPtr, pCon);
      }
    }
  } else {
    iRet = EVCGetPos(pEva, pCon, &fTemp);
  }

  /* get mf */
  fMF = -0.00;
  pEva = (pEVControl) FindCommandData(pServ->pSics, "mf",
                                      "Environment Controller");
  if (pEva == NULL) {
    pPtr = (pDummy) FindCommandData(pServ->pSics, "mf",
                                    "RemObject");
    if (pPtr != NULL) {
      pDriv = pPtr->pDescriptor->GetInterface(pPtr, DRIVEID);
      if (pDriv != NULL) {
        fMF = pDriv->GetValue(pPtr, pCon);
      }
    }
  } else {
    iRet = EVCGetPos(pEva, pCon, &fMF);
  }

  /* write profile */
  if (priv->profFile) {
    /* collect data */
    SNXFormatTime(pBueffel, 512);
    GetScanVarStep(priv->pScanner, 0, &fStep);
    fPreset = GetScanPreset(priv->pScanner);
    prot = getProtonAverage(priv);
    /*
      They rather wanted fMF in place of the proton average which fell victim to the 
      monitor assignement chaos anyway.
    */
    if(extra == NULL){
    	fprintf(priv->profFile, "%3d %7.4f %9.0f %7.3f %12f %s\n", iNP, fStep,
    			fPreset, fTemp, fMF, pBueffel);
    } else {
    	fprintf(priv->profFile, "%3d %7.4f %9.0f %7.3f %12f %s %s\n", iNP, fStep,
		fPreset, fTemp, fMF, pBueffel,extra);
    }
    for (i = 0; i < iNP; i++) {
      for (ii = 0; ii < 10 && i < iNP; ii++) {
        fprintf(priv->profFile, " %7ld", lCounts[i]);
        iLF = 1;
        i++;
      }
      fprintf(priv->profFile, "\n");
      i--;
      iLF = 0;
    }
    if (iLF) {
      fprintf(priv->profFile, "\n");
    }
    fflush(priv->profFile);
  }

  strcpy(pTime, pBueffel);
  sprintf(pBueffel, "%3d%8.4f%10.0f%8.3f %s\n", iNP, fStep,
          fPreset, fTemp, pTime);
  SCWrite(pCon, pBueffel, eLog);
  pBueffel[0] = '\0';
  for (i = 0; i < iNP; i++) {
    for (ii = 0; ii < 10 && i < iNP; ii++) {
      sprintf(pNum, " %6ld", lCounts[i]);
      strcat(pBueffel, pNum);
      iLF = 1;
      i++;
    }
    SCWrite(pCon, pBueffel, eLog);
    pBueffel[0] = '\0';
    i--;
    iLF = 0;
  }
  if (iLF) {
    SCWrite(pCon, pBueffel, eLog);
  }
  free(lCounts);
  return 1;

}

/*----------------------------------------------------------------------------*/
static int FourMessStore(pSICSOBJ self, SConnection * pCon,
                         pHdb commandNode, pHdb par[], int nPar)
{
  double fHkl[3], fPosition[4];
  int i;

  if (nPar < 7) {
    SCWrite(pCon, "ERROR: not enough arguments for store", eError);
    return 0;
  }

  /* load hkl */
  for (i = 0; i < 3; i++) {
    fHkl[i] = par[i]->value.v.doubleValue;
  }

  /* load positions */
  for (i = 0; i < 4; i++) {
    fPosition[i] = par[i + 3]->value.v.doubleValue;
  }

  return FourMessStoreIntern(self, pCon, fHkl, fPosition, NULL);
}

/*----------------------------------------------------------------------------*/
static int FourMessStoreExtra(pSICSOBJ self, SConnection * pCon,
                         pHdb commandNode, pHdb par[], int nPar)
{
  double fHkl[3], fPosition[4];
  int i;

  if (nPar < 8) {
    SCWrite(pCon, "ERROR: not enough arguments for storeextra", eError);
    return 0;
  }

  /* load hkl */
  for (i = 0; i < 3; i++) {
    fHkl[i] = par[i]->value.v.doubleValue;
  }

  /* load positions */
  for (i = 0; i < 4; i++) {
    fPosition[i] = par[i + 3]->value.v.doubleValue;
  }

  return FourMessStoreIntern(self, pCon, fHkl, fPosition, par[7]->value.v.text);
}

/*------------------------------------------------------------------*/
static int weakScan(pSICSOBJ self, SConnection * pCon)
{
  int i, np;
  long low = 99999, high = -99999, *lCounts = NULL;
  pFourMess priv = self->pPrivate;
  hdbValue v;

  /*
     the scan is always OK if we do not test for weak conditions or we are in psd mode
   */
  SICSHdbGetPar(self, pCon, "weak", &v);
  if (v.v.intValue == 0) {
    return 0;
  }

  np = GetScanNP(priv->pScanner);
  lCounts = malloc(np * sizeof(long));
  if (lCounts == NULL) {
    SCWrite(pCon, "ERROR: out of memory in weakScan test", eLogError);
    return 0;
  }
  GetScanCounts(priv->pScanner, lCounts, np);
  for (i = 0; i < np; i++) {
    if (lCounts[i] < low) {
      low = lCounts[i];
    }
    if (lCounts[i] > high) {
      high = lCounts[i];
    }
  }
  /*
     I am using the weakest point here as a rough estimate of
     the background
   */
  SICSHdbGetPar(self, pCon, "weakthreshold", &v);
  if (high - 2 * low > v.v.intValue) {
    return 0;
  } else {
    return 1;
  }
}

/*----------------------------------------------------------------------------*/
static int FourMessWeak(pSICSOBJ self, SConnection * pCon,
                        pHdb commandNode, pHdb par[], int nPar)
{
  int weak;

  weak = weakScan(self, pCon);
  SCPrintf(pCon, eLog, "weak = %d", weak);
  return 1;
}

/*----------------------------------------------------------------------------
 * Usage of the suppress flag:
 *  0 do notthing
 *  1 suppress symmetriqually equivalent reflections
 *  2 opposite: add only reflections NOT allowed by the spacegroup
 * ---------------------------------------------------------------------------*/
static int GenIndex(pSICSOBJ self, SConnection * pCon, pHdb commandNode,
                    pHdb par[], int nPar)
{
  double lambda = SXGetLambda(), d, om, hkl[3];
  int h, k, l, minh, mink, minl, suppress;
  hdbValue hkllim, sttlim;
  T_SgInfo *sginfo = SXGetSpaceGroup();
  pFourMess priv = self->pPrivate;
  int count = 0;
  MATRIX B, H, Z1;
  const double *cell;
  lattice direct;

  if (nPar < 1) {
    SCWrite(pCon, "ERROR: need a suppression flag for indgen", eError);
    return 0;
  }
  SICSHdbGetPar(self, pCon, "hkllim", &hkllim);
  SICSHdbGetPar(self, pCon, "sttlim", &sttlim);
  suppress = par[0]->value.v.intValue;

  minh = hkllim.v.intArray[3];
  mink = hkllim.v.intArray[4];
  minl = hkllim.v.intArray[5];

  SetListMin_hkl(sginfo, hkllim.v.intArray[1], hkllim.v.intArray[2],
                 &minh, &mink, &minl);
  ClearReflectionList(priv->messList);
  cell = SXGetCell();
  direct.a = cell[0];
  direct.b = cell[1];
  direct.c = cell[2];
  direct.alpha = cell[3];
  direct.beta = cell[4];
  direct.gamma = cell[5];
  B = mat_creat(3, 3, UNIT_MATRIX);
  if (!calculateBMatrix(direct, B)) {
    SCWrite(pCon, "ERROR: invalid cell", eError);
    return 0;
  }
  H = mat_creat(3, 1, ZERO_MATRIX);


  for (h = hkllim.v.intArray[0]; h <= hkllim.v.intArray[3]; h++) {
    for (k = hkllim.v.intArray[1]; k <= hkllim.v.intArray[4]; k++) {
      for (l = hkllim.v.intArray[2]; l <= hkllim.v.intArray[5]; l++) {
    	  /* SCPrintf(pCon,eLog, "Testing %d, %d, %d", h,k,l); */
        /* first test: extinct */
    	  
    	if(suppress != 2){  
    		if (IsSysAbsent_hkl(sginfo, h, k, l, NULL) != 0) {
    			/* SCPrintf(pCon,eLog, "%d, %d, %d rejected for sys absent", h, k, l); */
    			continue;
    		}
    	} else {
    		if (!IsSysAbsent_hkl(sginfo, h, k, l, NULL) != 0) {
    			/* SCPrintf(pCon,eLog, "%d, %d, %d rejected for not sys absent", h, k, l); */
    			continue;
    		}
    	}
        /* second test: a symmetrically equivalent already seen */
        if ((suppress > 0) && IsSuppressed_hkl(sginfo, minh, mink, minl,
                                                hkllim.v.intArray[1],
                                                hkllim.v.intArray[2], h, k,
                                                l) != 0) {
        	/* SCPrintf(pCon,eLog, "%d, %d, %d rejected for symmetrical eqiv", h, k, l); */
          continue;
        }
        /* third test: within stt limits */
        H[0][0] = (double) h;
        H[1][0] = (double) k;
        H[2][0] = (double) l;
        Z1 = mat_mul(B, H);
        calcTheta(lambda, Z1, &d, &om);
        om *= 2.;
        mat_free(Z1);
        if (om > sttlim.v.floatArray[0] && om < sttlim.v.floatArray[1]) {
          hkl[0] = (double) h;
          hkl[1] = (double) k;
          hkl[2] = (double) l;
          AddRefIdx(priv->messList, hkl);
          count++;
        } 
        if(SCGetInterrupt(pCon) != eContinue){
        	SCWrite(pCon,"ERROR: reflection generation aborted", eError);
        	return 0;
        }
      }
    }
  }
  mat_free(B);
  mat_free(H);
  priv->masterCount = count;
  SCPrintf(pCon, eValue, "%d reflections generated", count);
  return 1;
}

/*-----------------------------------------------------------------------------*/
static int GenInconsumerate(pSICSOBJ self, SConnection * pCon,
                            pHdb commandNode, pHdb par[], int nPar)
{
  double hkl[3], qvec[3];
  pFourMess priv = self->pPrivate;
  int i, j, iGen = 0, startCount;

  if (nPar < 3) {
    SCWrite(pCon,
            "ERROR: need q displacement vector with three components",
            eError);
    return 0;
  }
  qvec[0] = par[0]->value.v.doubleValue;
  qvec[1] = par[1]->value.v.doubleValue;
  qvec[2] = par[2]->value.v.doubleValue;

  startCount = priv->masterCount;
  for (i = 0; i < startCount; i++) {
    GetRefIndex(priv->messList, i, hkl);
    for (j = 0; j < 3; j++) {
      hkl[j] += qvec[j];
    }
    AddRefIdx(priv->messList, hkl);
    iGen++;
    GetRefIndex(priv->messList, i, hkl);
    for (j = 0; j < 3; j++) {
      hkl[j] -= qvec[j];
    }
    if(FindHKL(priv->messList, hkl[0], hkl[1], hkl[2]) == NULL){
      AddRefIdx(priv->messList, hkl);
      iGen++;
    }
    if(SCGetInterrupt(pCon) != eContinue){
    	SCWrite(pCon,"ERROR: generating incommensurate reflections aborted", eError);
    	return 0;
    }
    if( (i % 50) == 0 ){
      SCPrintf(pCon,eLog, "%d of %d input reflections processed", i, startCount);
    }
  }
  /*
    add satellites of 0,0,0,
  */
  for(i = 0; i < 3; i++){
    hkl[i] = 0. + qvec[i];
  }
  AddRefIdx(priv->messList, hkl);
  iGen++;
  for(i = 0; i < 3; i++){
    hkl[i] = 0. - qvec[i];
  }
  if(FindHKL(priv->messList, hkl[0], hkl[1], hkl[2]) == NULL){
    AddRefIdx(priv->messList, hkl);
    iGen++;
  }
  SCPrintf(pCon, eValue,
           "%d additional inconsumerate reflections generated", iGen);
  return 1;
}

/*----------------------------------------------------------------------------*/
static int hklCompare(const void *h1, const void *h2)
{
  const double *hkl1 = h1, *hkl2 = h2;

  if (hkl1[3] < hkl2[3]) {
    return -1;
  } else if (hkl1[3] == hkl2[3]) {
    return 0;
  } else {
    return 1;
  }
}

/*----------------------------------------------------------------------------
 * This sorts on two theta only.
 -----------------------------------------------------------------------------*/
static int SortRef(pSICSOBJ self, SConnection * pCon, pHdb commandNode,
                   pHdb par[], int nPar)
{
  double *sortlist, d, lambda, om, hkl[4], ang[4];
  const double *cell;
  double *hklc;
  int nRefl, i, j;
  MATRIX B, H, Z1;
  lattice direct;
  pFourMess priv = self->pPrivate;

  lambda = SXGetLambda();
  cell = SXGetCell();
  direct.a = cell[0];
  direct.b = cell[1];
  direct.c = cell[2];
  direct.alpha = cell[3];
  direct.beta = cell[4];
  direct.gamma = cell[5];
  B = mat_creat(3, 3, UNIT_MATRIX);
  if (!calculateBMatrix(direct, B)) {
    SCWrite(pCon, "ERROR: invalid cell", eError);
    return 0;
  }
  H = mat_creat(3, 1, ZERO_MATRIX);

  nRefl = ReflectionListCount(priv->messList);
  sortlist = malloc(nRefl * 4 * sizeof(double));
  if (sortlist == NULL) {
    SCWrite(pCon, "ERROR: out of memory in SortRef", eError);
    return 0;
  }
  /*
   * I am using hkl[3] for storing the theta value to sort for! */
  for (i = 0; i < nRefl; i++) {
    GetRefIndex(priv->messList, i, hkl);
    for (j = 0; j < 3; j++) {
      H[j][0] = hkl[j];
    }
    Z1 = mat_mul(B, H);
    calcTheta(lambda, Z1, &d, &om);
    mat_free(Z1);
    hkl[3] = om;
    memcpy(sortlist + i * 4, hkl, 4 * sizeof(double));
  }

  qsort(sortlist, nRefl, 4 * sizeof(double), hklCompare);

  ClearReflectionList(priv->messList);
  for (i = 0; i < nRefl; i++) {
    ang[1] = sortlist[i * 4 + 3];
    ang[0] = 2. * ang[1];
    ang[2] = .0;
    ang[3] = .0;
    /* AddRefIdxAng(priv->messList, sortlist+i*4,ang); */
    AddRefIdx(priv->messList, sortlist + i * 4);
  }
  free(sortlist);
  mat_free(B);
  mat_free(H);
  SCSendOK(pCon);
  return 1;
}
/*----------------------------------------------------------------------------
 * Sorting reflection is actually a variant of the traveling salesman
 * problem. Optimal solutions for this are computationally expensive.
 * Optimise for the shortest path between reflections. Not all motors
 * are equal though. Om and phi move quick whereas stt and chi go slow.
 * As a compromise, what I try here is sorting by a weighted addition of
 * stt and chi.
 ----------------------------------------------------------------------------*/
static int SortRefNew(pSICSOBJ self, SConnection * pCon, pHdb commandNode,
                   pHdb par[], int nPar)
{
  double *sortlist, d, lambda, om, hkl[4], ang[5];
  const double *cell;
  double *hklc;
  int nRefl, i, j, status, count;
  SingleXModes mode;
  pFourMess priv = self->pPrivate;
  pSingleDiff diffi = NULL;

  diffi = SXGetDiffractometer();
  mode = SXGetMode();
  if(diffi == NULL){
	  SCWrite(pCon,"ERROR: not initialized, no diffractometer", eError);
	  return 0;
  }

  nRefl = ReflectionListCount(priv->messList);
  sortlist = malloc(nRefl * 4 * sizeof(double));
  if (sortlist == NULL) {
    SCWrite(pCon, "ERROR: out of memory in SortRef", eError);
    return 0;
  }

  /*
   * I am using hkl[3] for storing the sort value to sort for! */
  for (i = 0, count = 0; i < nRefl; i++) {
    GetRefIndex(priv->messList, i, hkl);
    status = diffi->calculateSettings(diffi,hkl,ang);
    if(status == 1){
    	if(mode == Bisecting){
    		hkl[3] = ang[0] + 0.5* ang[2];
    		if(ang[2] < 0){
    			SCWrite(pCon,"WARNING: chi < 0!", eWarning);
    		}
    	} else if(mode == NB){
    		/*
    		 * sort for nu in the first place. In order to cope with negativity, add
    		 * 10.
    		 */
    		hkl[3] = ang[2]+10 + 0.1 * ang[0];
    	} else if(mode == BiNB){
    		/*
    		 * sort for nu in the first place. In order to cope with negativity, add
    		 * 10.
    		 */
    		hkl[3] = ang[4]+10 + 0.1 * ang[0];
    	} else {
    		hkl[3] = ang[0];
    	}
    	memcpy(sortlist + count * 4, hkl, 4 * sizeof(double));
    	count++;
    }
  }

  qsort(sortlist, count, 4 * sizeof(double), hklCompare);

  ClearReflectionList(priv->messList);
  for (i = 0; i < count; i++) {
	  hklc = sortlist+i*4;
	  hklc[3] = .0; /* otherwise this will be interpreted as psi! */
	  status = diffi->calculateSettings(diffi,hklc,ang);
	  AddRefIdxAng(priv->messList, sortlist+i*4,ang);
	  /* AddRefIdx(priv->messList, sortlist + i * 4); */
  }
  free(sortlist);
  SCSendOK(pCon);
  return 1;
}
/*----------------------------------------------------------------------------
 * Write a reflection list in the TSPLIB format for running a traveling
 * salesman algorithm against it.
----------------------------------------------------------------------------- */
static int WriteTSP(pSICSOBJ self, SConnection * pCon, pHdb commandNode,
                   pHdb par[], int nPar)
{
  double hkl[4], ang[4];
  int nRefl, i, j, status, count;
  SingleXModes mode;
  pFourMess priv = self->pPrivate;
  pSingleDiff diffi = NULL;
  FILE *fd = NULL;

  if(nPar < 1) {
	  SCWrite(pCon,"ERROR: need a filename parameter", eError);
	  return 0;
  }
  fd = fopen(trim(par[0]->value.v.text),"w");
  if(fd == NULL){
	SCPrintf(pCon,eError,"ERROR: failed to open %s for writing", par[0]->value.v.text);
	return 0;
  }

  diffi = SXGetDiffractometer();
  mode = SXGetMode();
  if(diffi == NULL){
	  SCWrite(pCon,"ERROR: not initialized, no diffractometer", eError);
	  return 0;
  }

  nRefl = ReflectionListCount(priv->messList);
  /*
   * smear header...
   */
  fprintf(fd,"NAME : %s\n", par[0]->value.v.text);
  fprintf(fd,"TYPE : TSP\n");
  fprintf(fd,"DIMENSION : %d\n", nRefl);
  fprintf(fd,"EDGE_WEIGHT_TYPE : EUC_3D\n");
  fprintf(fd,"NODE_COORD_SECTION\n");


  for (i = 0, count = 0; i < nRefl; i++) {
    GetRefIndex(priv->messList, i, hkl);
    status = diffi->calculateSettings(diffi,hkl,ang);
    if(status == 1){
    	if(mode == Bisecting){
    		fprintf(fd,"%d %e %e %e\n",i+1,ang[0]*3*100,ang[2]*2*100,ang[3]*100);
    	} else if(mode == NB){
    		fprintf(fd,"%d %e %e %e\n",i,ang[2]*4*100,ang[0]*3*100,ang[1]*100);
    	} else {
    		fprintf(fd,"%d %e %e %e\n",i,ang[0]*3*100,ang[2]*2*100,ang[3]*2*100);
    	}
    }
  }
  fclose(fd);

  SCSendOK(pCon);
  return 1;
}
/*----------------------------------------------------------------------------
 * read a tour file as generated by LKH or other TSP solvers and
 * reorder the reflection list accordingly.
 -----------------------------------------------------------------------------*/
static int ReadTour(pSICSOBJ self, SConnection * pCon, pHdb commandNode,
                   pHdb par[], int nPar)
{
	FILE *fd = NULL;
	char buffer[132];
	double hkl[4], ang[5];
	int idx, tour;
	pFourMess priv = self->pPrivate;
	pSingleDiff diffi = NULL;


	if(nPar < 1) {
		SCWrite(pCon,"ERROR: need name of a tour file to read", eError);
		return 0;
	}
	fd = fopen(par[0]->value.v.text,"r");
	if(fd == NULL){
		SCPrintf(pCon,eError,"ERROR: failed to open tour file %s for reading",
				par[0]->value.v.text);
		return 0;
	}
	tour = LLDcreate(sizeof(hkl));
	diffi = SXGetDiffractometer();
	assert(diffi != NULL);

	/* skip header */
	while(fgets(buffer,sizeof(buffer),fd) != NULL){
		if(strstr(buffer,"TOUR_SECTION") != NULL){
			break;
		}
	}
	while(fscanf(fd,"%d",&idx) == 1){
		if(idx < 0){
			break;
		}
		GetRefIndex(priv->messList,idx-1,hkl);
		LLDnodeAppendFrom(tour,hkl);
	}
	fclose(fd);

	idx = LLDnodePtr2First(tour);
	ClearReflectionList(priv->messList);
	while(idx == 1){
		LLDnodeDataTo(tour,hkl);
		hkl[3] = 0;
		diffi->calculateSettings(diffi,hkl,ang);
		AddRefIdxAng(priv->messList, hkl,ang);
		idx = LLDnodePtr2Next(tour);
	}
	LLDdelete(tour);
	SCSendOK(pCon);
	return 1;
}
/*----------------------------------------------------------------------------*/
static int FourMessSave(void *data, char *name, FILE * fd)
{
  pSICSOBJ self = data;
  pFourMess priv = self->pPrivate;

  SaveSICSOBJ(data, name, fd);
  priv->stepTable->pDes->SaveStatus(priv->stepTable,"fmess table", fd);
  return 1;
}

/*----------------------------------------------------------------------------*/
void InstallFourMess(SConnection * pCon, SicsInterp * pSics)
{
  pFourMess priv = NULL;
  pSICSOBJ pNew = NULL;
  pHdb cmd = NULL;
  int hkl[] = { -10, -10, 10, 10, 10, 10 };
  double sttlim[] = { 5.0, 180 };

  pNew = MakeSICSOBJ("fmess", "FourMess");
  priv = calloc(1, sizeof(FourMess));
  if (pNew == NULL || priv == NULL) {
    SCWrite(pCon, "ERROR: out of memory creating fourmess", eError);
  }
  pNew->pDes->SaveStatus = FourMessSave;
  pNew->pPrivate = priv;
  pNew->KillPrivate = KillFourMess;

  priv->stepTable = MakeFourCircleTable();
  AddCommand(pSics, "fmesstable", InterInvokeSICSOBJ, NULL, priv->stepTable);

  priv->messList = CreateReflectionList(pCon, pSics, "messref");
  if (priv->stepTable < 0 || priv->messList == NULL) {
    SCWrite(pCon, "ERROR: out of memory creating fourmess", eError);
  }

  cmd = AddSICSHdbPar(pNew->objectNode, "weak", usUser, MakeHdbInt(0));
  SetHdbProperty(cmd, "__save", "true");
  cmd = AddSICSHdbPar(pNew->objectNode, "weakthreshold", usUser,
                      MakeHdbInt(20));
  SetHdbProperty(cmd, "__save", "true");
  
  cmd = AddSICSHdbPar(pNew->objectNode, "fast", usUser, MakeHdbInt(0));
  SetHdbProperty(cmd, "__save", "true");
  
  cmd =
      AddSICSHdbPar(pNew->objectNode, "mode", usUser,
                    MakeHdbText("Monitor"));

  cmd =
      AddSICSHdbPar(pNew->objectNode, "hkllim", usUser,
                    MakeHdbIntArray(6, hkl));
  SetHdbProperty(cmd, "__save", "true");
  cmd =
      AddSICSHdbPar(pNew->objectNode, "sttlim", usUser,
                    MakeHdbFloatArray(2, sttlim));
  SetHdbProperty(cmd, "__save", "true");

  cmd =
      AddSICSHdbPar(pNew->objectNode, "start", usUser,
                    MakeSICSFunc(FourMessStart));
  cmd = AddSICSHdbPar(cmd, "filename", usUser, MakeHdbText("Unknown"));
  cmd =
      AddSICSHdbPar(pNew->objectNode, "close", usUser,
                    MakeSICSFunc(FourMessClose));
  cmd =
      AddSICSHdbPar(pNew->objectNode, "reopen", usUser,
                    MakeSICSFunc(FourMessReopen));

  cmd =
      AddSICSHdbPar(pNew->objectNode, "prepare", usUser,
                    MakeSICSFunc(FourMessPrepareCmd));

  cmd =
      AddSICSHdbPar(pNew->objectNode, "scanpar", usUser,
                    MakeSICSFunc(FourMessScanPar));
  cmd = AddSICSHdbPar(cmd, "twotheta", usUser, MakeHdbFloat(.0));

  cmd =
      AddSICSHdbPar(pNew->objectNode, "store", usUser,
                    MakeSICSFunc(FourMessStore));
  AddSICSHdbPar(cmd, "h", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "k", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "l", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "stt", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "om", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "chi", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "phi", usUser, MakeHdbFloat(.0));

  cmd =
      AddSICSHdbPar(pNew->objectNode, "storeextra", usUser,
                    MakeSICSFunc(FourMessStoreExtra));
  AddSICSHdbPar(cmd, "h", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "k", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "l", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "stt", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "om", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "chi", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "phi", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "extra", usUser, MakeHdbText(""));
  
  cmd =
      AddSICSHdbPar(pNew->objectNode, "weak", usUser,
                    MakeSICSFunc(FourMessWeak));

  cmd =
      AddSICSHdbPar(pNew->objectNode, "indgen", usUser,
                    MakeSICSFunc(GenIndex));
  AddSICSHdbPar(cmd, "sup", usUser, MakeHdbInt(1));

  cmd =
      AddSICSHdbPar(pNew->objectNode, "genw", usUser,
                    MakeSICSFunc(GenInconsumerate));
  AddSICSHdbPar(cmd, "hw", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "kw", usUser, MakeHdbFloat(.0));
  AddSICSHdbPar(cmd, "lw", usUser, MakeHdbFloat(.0));

  cmd =
      AddSICSHdbPar(pNew->objectNode, "indsort", usUser,
                    MakeSICSFunc(SortRef));
  cmd =
      AddSICSHdbPar(pNew->objectNode, "writetsp", usUser,
                    MakeSICSFunc(WriteTSP));
  AddSICSHdbPar(cmd, "filename", usUser, MakeHdbText("Unknown"));

  cmd =
      AddSICSHdbPar(pNew->objectNode, "readtour", usUser,
                    MakeSICSFunc(ReadTour));
  AddSICSHdbPar(cmd, "filename", usUser, MakeHdbText("Unknown"));

  cmd =
      AddSICSHdbPar(pNew->objectNode, "template", usMugger,
                    MakeHdbText("Unknown"));
  cmd =
      AddSICSHdbPar(pNew->objectNode, "scanobj", usMugger,
                    MakeHdbText("xxxscan"));
  PrependHipadabaCallback(cmd,
                          MakeHipadabaCallback(SetScannerCB, priv, NULL));

  priv->pScanner = FindCommandData(pSics, "xxxscan", "ScanObject");
  /* AddHipadabaChild(pNew->objectNode, priv->stepTable->objectNode,pCon); */

  AddCommand(pSics, "fmess", FourMessAction, KillSICSOBJ, pNew);
}