From fdc784eb54f0036dd8cc17cc3b3dda5d4675715e Mon Sep 17 00:00:00 2001
From: Douglas Clowes <dcl@ansto.gov.au>
Date: Thu, 29 Nov 2012 10:01:30 +1100
Subject: [PATCH] data smoothing plus optional scan object

---
 fitcenter.c | 72 +++++++++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 65 insertions(+), 7 deletions(-)

diff --git a/fitcenter.c b/fitcenter.c
index 0e4aa9ae..a91ae0ed 100644
--- a/fitcenter.c
+++ b/fitcenter.c
@@ -20,6 +20,10 @@
 #include "scan.h"
 #include "fitcenter.h"
 #define THRESHOLD .1
+
+float ggf1[] = { 0.006, 0.061, 0.242, 0.383, 0.242, 0.061, 0.006 };     // Gaussian Smooth Filter
+float ggf2[] = { 0.1429, 0.1429, 0.1429, 0.1429, 0.1429, 0.1429, 0.1429 };      // Mean Smooth Filter
+
 /*--------------------------------------------------------------------------*/
 typedef struct __FitCenter {
   pObjectDescriptor pDes;
@@ -73,6 +77,28 @@ void DeleteFitCenter(void *pData)
   free(self);
 }
 
+/*--------------------------------------------------------------------------*/
+void SmoothScanCounts(long *lData, int iDataLen)
+{
+  int ind = (int) (sizeof(ggf1) / sizeof(float) - 1) / 2;
+  long *pData = (long *) malloc(iDataLen * sizeof(long));
+
+  for (int i = 0; i < iDataLen; i++)
+    pData[i] = lData[i];
+
+  for (int i = 0; i < iDataLen; i++) {
+    lData[i] = pData[i] * ggf1[ind];
+    for (int j = 1; j <= ind; j++) {
+      if ((i - j) > 0)
+        lData[i] += pData[i - j] * ggf1[ind - j];
+      if ((i + j) < iDataLen)
+        lData[i] += pData[i + j] * ggf1[ind + j];
+    }
+  }
+
+  free(pData);
+}
+
 /*--------------------------------------------------------------------------*/
 static int Init(pFit self)
 {
@@ -97,6 +123,9 @@ static int Init(pFit self)
   }
 
   GetScanCounts(self->pScan, self->lCounts, self->iNP);
+
+  SmoothScanCounts(self->lCounts, self->iNP);
+
   GetScanVar(self->pScan, 0, self->fAxis, self->iNP);
   GetScanVarName(self->pScan, 0, self->pName, 131);
 
@@ -433,11 +462,40 @@ int FitWrapper(SConnection * pCon, SicsInterp * pSics, void *pData,
   int iRet;
   char pBueffel[256];
   pDynString buf = NULL;
+  pScanData pScan = NULL;
+  char *cmd = NULL;
 
   self = (pFit) pData;
   assert(self);
-
-  iRet = CalculateFit(self);
+  /*
+   * allow for an optional scan object as first argument
+   * and make cmd point at the command
+   */
+  if (argc > 1) {
+    cmd = argv[1];
+    CommandList *pCom = NULL;
+    pCom = FindCommand(pSics, argv[1]);
+    if (pCom) {
+      pDummy pDum = NULL;
+      pDum = (pDummy) pCom->pData;
+      if (pDum) {
+        if (strcmp(pDum->pDescriptor->name, "ScanObject") == 0) {
+          pScan = (pScanData) pDum;
+          if (argc > 2)
+            cmd = argv[2];
+        }
+      }
+    }
+  }
+  if (pScan) {
+    pScanData pScanSave;
+    pScanSave = self->pScan;
+    self->pScan = pScan;
+    iRet = CalculateFit(self);
+    self->pScan = pScanSave;
+  } else {
+    iRet = CalculateFit(self);
+  }
   switch (iRet) {
   case 0:
     SCWrite(pCon, "ERROR: failure to fit your data!", eError);
@@ -470,14 +528,14 @@ int FitWrapper(SConnection * pCon, SicsInterp * pSics, void *pData,
      This is a little feature to get the peak without rubbish for
      the TAS routines
    */
-  if (argc > 1) {
-    strtolower(argv[1]);
-    if (strcmp(argv[1], "value") == 0) {
-      snprintf(pBueffel,sizeof(pBueffel)-1, "%f", self->fCenter);
+  if (cmd) {
+    strtolower(cmd);
+    if (strcmp(cmd, "value") == 0) {
+      sprintf(pBueffel, "%f", self->fCenter);
       SCWrite(pCon, pBueffel, eValue);
       return 1;
     }
-    if (strcmp(argv[1], "data") == 0) {
+    if (strcmp(cmd, "data") == 0) {
       snprintf(pBueffel, 255, "%f,%f,%ld",
                self->fCenter, self->FWHM, self->lPeak);
       SCWrite(pCon, pBueffel, eValue);