text: /**
   * @file yapp_dedisperse.c
   * Program to dedisperse the input signal for the given value of DM.
   *
   * @verbatim
   * Usage: yapp_dedisperse [options] <data-file>
   *     -h  --help                           Display this usage information
   *     -s  --skip <time>                    The length of data in seconds, to
   *                                          be skipped
   *                                          (default is 0 s)
   *     -p  --proc <time>                    The length of data in seconds, to
   *                                          be processed
   *                                          (default is all)
   *     -n  --nsamp <samples>                Number of samples read in one block
   *                                          (default is 4096 samples)
   *     -d  --dm <dm>                        DM at which to de-disperse
   *                                          (default is 10.0)
   *     -l  --law <law>                      Dispersion law
   *                                          (default is 2.0)
   *     -b  --nsubband <nsubband>            Number of sub-bands
   *                                          (must be < number of channels)
   *     -u  --subband <subband>              Sub-band to dedisperse
   *                                          (within [0, nsubband))
   *     -o  --out-format <format>            Output format - 'dds', 'tim', or
   *                                          'fil'
   *                                          (default is 'tim')
   *     -f  --f-decimate <factor>            Decimation-in-frequency factor,
   *                                          for 'fil' output
   *     -g  --graphics                       Turn on plotting
   *     -m  --colour-map <name>              Colour map for plotting
   *                                          (default is 'jet')
   *     -i  --invert                         Invert the background and
   *                                          foreground colours in plots
   *     -e  --non-interactive                Run in non-interactive mode
   *     -v  --version                        Display the version @endverbatim
   *
   * @author Jayanth Chennamangalam
   * @date 2008.11.14
   */
  
  #include "yapp.h"
  #include "yapp_sigproc.h"   /* for SIGPROC filterbank file format support */
  #include "colourmap.h"
  
  /**
   * The build version string, maintained in the file version.c, which is
   * generated by makever.c.
   */
  extern const char *g_pcVersion;
  
  /* PGPLOT device ID */
  extern int g_iPGDev;
  
  /* data file */
  extern FILE *g_pFData;
  
  /* the following are global only to enable cleaning up in case of abnormal
     termination, such as those triggered by SIGINT or SIGTERM */
  int *g_piOffsetTab = NULL;
  int *g_piOutOffsetTab = NULL;
  char *g_pcIsTimeGood = NULL;
  float *g_pfBFTimeSectMean = NULL;
  float *g_pfBFGain = NULL;
  double (*g_padBadTimes)[][NUM_BAD_BOUNDS] = NULL;
  float *g_pfBuf0 = NULL;
  float *g_pfBuf1 = NULL;
  float *g_pfPlotBuf = NULL;
  float *g_pfDedispData = NULL;
  float *g_pfDecDedispData = NULL;
  float *g_pfXAxis = NULL;
  float *g_pfYAxis = NULL;
  
  int main(int argc, char *argv[])
  {
      FILE *pFDedispData = NULL;
      char *pcFileSpec = NULL;
      char acFileDedisp[LEN_GENSTRING] = {0};
      int iFormat = DEF_FORMAT;
      int iOutputFormat = DEF_OUT_FORMAT;
      char cIsPlotToFile = YAPP_FALSE;
      double dDataSkipTime = 0.0;
      double dDataProcTime = 0.0;
      YUM_t stYUM = {{0}};
      YUM_t stYUMOut = {{0}};
      char cIsDMGiven = YAPP_FALSE;
      double dDM = 0.0;
      float fLaw = DEF_LAW;
      float fChanBW = 0.0;
      int iMaxOffset = 0;
      int iOutMaxOffset = 0;
      int iNumChans = 0;
      int iOutNumChans = 0;
      float fSampSize = 0.0;      /* number of bits that make a sample */
      int iTotSampsPerBlock = 0;  /* iNumChans * iBlockSize */
      int iEffcNumGoodChans = 0;
      float fStatBW = 0.0;
      float fNoiseRMS = 0.0;
      double dTNextBF = 0.0;
      double dTSamp = 0.0;        /* holds sampling time in ms */
      double dTSampInSec = 0.0;   /* holds sampling time in s */
      double dTNow = 0.0;
      int iTimeSect = 0;
      int iBadTimeSect = 0;
      char cIsInBadTimeRange = YAPP_FALSE;
      float *pfTimeSectGain = NULL;
      float *pfPriBuf = NULL;
      float *pfSecBuf = NULL;
      float *pfSpectrum = NULL;
      float *pfOutSpec = NULL;
      float *pfOffsetSpec = NULL;
      int iPrimaryBuf = BUF_0;
      int iOffset = 0;
      long int lBytesToSkip = 0;
      long int lBytesToProc = 0;
      int iTimeSamps = 0;
      int iTimeSampsToSkip = 0;
      int iTimeSampsToProc = 0;
      int iBlockSize = DEF_SIZE_BLOCK;
      int iNumReads = 0;
      int iReadBlockCount = 0;
      int iReadSmpCount = 0;         /* (iReadBlockCount - 1) * iBlockSize */
      char cIsLastBlock = YAPP_FALSE;
      long int lDataSizeTotal = 0;
      int iRet = YAPP_RET_SUCCESS;
      float fDataMin = 0.0;
      float fDataMax = 0.0;
      int iReadItems = 0;
      char acDev[LEN_GENSTRING] = {0};
      char *pcFilename = NULL;
      int iNumSamps = 0;
      int iDiff = 0;
      int iNumSubBands = 0;
      int iSubBand = 0;
      int iChansPerSubBand = 0;
      int iStartChan = 0;
      int iEndChan = 0;
      float fButX = 0.0;
      float fButY = 0.0;
      char cCurChar = 0;
      int i = 0;
      int j = 0;
      int k = 0;
      int l = 0;
      int m = 0;
      int n = 0;
      double dDelay = 0.0;
      int iStartOffset = 0;
      float fStartOffset = 0.0;
      int iDecFac = 0;
      char cHasGraphics = YAPP_FALSE;
      int iColourMap = DEF_CMAP;
      int iInvCols = YAPP_FALSE;
      char cIsNonInteractive = YAPP_FALSE;
      const char *pcProgName = NULL;
      int iNextOpt = 0;
      /* valid short options */
      const char* const pcOptsShort = "hs:p:n:d:l:b:u:o:f:gm:iev";
      /* valid long options */
      const struct option stOptsLong[] = {
          { "help",                   0, NULL, 'h' },
          { "skip",                   1, NULL, 's' },
          { "proc",                   1, NULL, 'p' },
          { "nsamp",                  1, NULL, 'n' },
          { "dm",                     1, NULL, 'd' },
          { "law",                    1, NULL, 'l' },
          { "nsubband",               1, NULL, 'b' },
          { "subband",                1, NULL, 'u' },
          { "out-format",             1, NULL, 'o' },
          { "f-decimate",             1, NULL, 'f' },
          { "graphics",               0, NULL, 'g' },
          { "colour-map",             1, NULL, 'm' },
          { "invert",                 0, NULL, 'i' },
          { "non-interactive",        0, NULL, 'e' },
          { "version",                0, NULL, 'v' },
          { NULL,                     0, NULL, 0   }
      };
  
      /* get the filename of the program from the argument list */
      pcProgName = argv[0];
  
      /* parse the input */
      do
      {
          iNextOpt = getopt_long(argc, argv, pcOptsShort, stOptsLong, NULL);
          switch (iNextOpt)
          {
              case 'h':   /* -h or --help */
                  /* print usage info and terminate */
                  PrintUsage(pcProgName);
                  return YAPP_RET_SUCCESS;
  
              case 's':   /* -s or --skip */
                  /* set option */
                  dDataSkipTime = atof(optarg);
                  break;
  
              case 'p':   /* -p or --proc */
                  /* set option */
                  dDataProcTime = atof(optarg);
                  break;
  
              case 'n':   /* -n or --nsamp */
                  /* set option */
                  iBlockSize = atoi(optarg);
                  /* validate - PGPLOT does not like iBlockSize = 1 */
                  if (iBlockSize < 2)
                  {
                      (void) fprintf(stderr,
                                     "ERROR: Number of samples must be > 1!\n");
                      PrintUsage(pcProgName);
                      return YAPP_RET_ERROR;
                  }
                  break;
  
              case 'd':   /* -d or --dm */
                  /* set option */
                  dDM = atof(optarg);
                  cIsDMGiven = YAPP_TRUE;
                  break;
  
              case 'l':   /* -l or --law */
                  /* set option */
                  fLaw = atof(optarg);
                  break;
  
              case 'b':   /* -b or --nsubband */
                  /* set option */
                  iNumSubBands = atoi(optarg);
                  break;
  
              case 'u':   /* -u or --subband */
                  /* set option */
                  iSubBand = atoi(optarg);
                  break;
  
              case 'o':   /* -o or --out-format */
                  /* set option */
                  if (0 == strcmp(optarg, YAPP_FORMATSTR_DTS_DDS))
                  {
                      iOutputFormat = YAPP_FORMAT_DTS_DDS;
                  }
                  else if (0 == strcmp(optarg, YAPP_FORMATSTR_DTS_TIM))
                  {
                      iOutputFormat = YAPP_FORMAT_DTS_TIM;
                  }
                  else if (0 == strcmp(optarg, YAPP_FORMATSTR_FIL))
                  {
                      iOutputFormat = YAPP_FORMAT_FIL;
                  }
                  else
                  {
                      (void) fprintf(stderr,
                                     "ERROR: Format should be 'dds' or "
                                     "'tim' or 'fil'!\n");
                      PrintUsage(pcProgName);
                      return YAPP_RET_ERROR;
                  }
                  break;
  
              case 'f':   /* -f or --f-decimate */
                  /* set option */
                  iDecFac = atoi(optarg);
                  break;
  
              case 'g':   /* -g or --graphics */
                  /* set option */
                  cHasGraphics = YAPP_TRUE;
                  break;
  
              case 'm':   /* -m or --colour-map */
                  /* set option */
                  iColourMap = GetColourMapFromName(optarg);
                  break;
  
              case 'i':  /* -i or --invert */
                  /* set option */
                  iInvCols = YAPP_TRUE;
                  break;
  
              case 'e':  /* -e or --non-interactive */
                  /* set option */
                  cIsNonInteractive = YAPP_TRUE;
                  break;
  
              case 'v':   /* -v or --version */
                  /* display the version */
                  (void) printf("%s\n", g_pcVersion);
                  return YAPP_RET_SUCCESS;
  
              case '?':   /* user specified an invalid option */
                  /* print usage info and terminate with error */
                  (void) fprintf(stderr, "ERROR: Invalid option!\n");
                  PrintUsage(pcProgName);
                  return YAPP_RET_ERROR;
  
              case -1:    /* done with options */
                  break;
  
              default:    /* unexpected */
                  assert(0);
          }
      } while (iNextOpt != -1);
  
      /* no arguments */
      if (argc <= optind)
      {
          (void) fprintf(stderr, "ERROR: Input file not specified!\n");
          PrintUsage(pcProgName);
          return YAPP_RET_ERROR;
      }
  
      if (!(cIsDMGiven))
      {
          (void) fprintf(stderr,
                         "ERROR: Required option not given! DM is required.\n");
          PrintUsage(pcProgName);
          return YAPP_RET_ERROR;
      }
  
      if (iNumSubBands > 0)
      {
          if (iSubBand >= iNumSubBands)
          {
              (void) fprintf(stderr,
                             "ERROR: Sub-band number inconsistent with number "
                             "of sub-bands!\n");
              PrintUsage(pcProgName);
              return YAPP_RET_ERROR;
          }
          if (iOutputFormat != YAPP_FORMAT_DTS_TIM)
          {
              (void) fprintf(stderr,
                             "ERROR: Sub-band dedispersion support only for "
                             ".tim output!\n");
              PrintUsage(pcProgName);
              return YAPP_RET_ERROR;
          }
      }
  
      if (iDecFac !=0 && iOutputFormat != YAPP_FORMAT_FIL)
      {
          (void) fprintf(stderr,
                         "WARNING: Ignoring decimation factor!\n");
          /* reset so that later checks evaluate to false */
          iDecFac = 0;
      }
  
      /* register the signal-handling function */
      iRet = YAPP_RegisterSignalHandlers();
      if (iRet != YAPP_RET_SUCCESS)
      {
          (void) fprintf(stderr,
                         "ERROR: Handler registration failed!\n");
          return YAPP_RET_ERROR;
      }
  
      /* get the input filename */
      pcFileSpec = argv[optind];
  
      /* determine the file type */
      iFormat = YAPP_GetFileType(pcFileSpec);
      if (YAPP_RET_ERROR == iFormat)
      {
          (void) fprintf(stderr,
                         "ERROR: File type determination failed!\n");
          return YAPP_RET_ERROR;
      }
      if (!((YAPP_FORMAT_FIL == iFormat) || (YAPP_FORMAT_SPEC == iFormat)))
      {
          (void) fprintf(stderr,
                         "ERROR: Invalid file type!\n");
          return YAPP_RET_ERROR;
      }
  
      /* read metadata */
      iRet = YAPP_ReadMetadata(pcFileSpec, iFormat, &stYUM);
      if (iRet != YAPP_RET_SUCCESS)
      {
          (void) fprintf(stderr,
                         "ERROR: Reading metadata failed for file %s!\n",
                         pcFileSpec);
          return YAPP_RET_ERROR;
      }
  
      dTSamp = stYUM.dTSamp;
      /* convert sampling interval to seconds */
      dTSampInSec = stYUM.dTSamp / 1e3;
      fChanBW = stYUM.fChanBW;
      iTimeSamps = stYUM.iTimeSamps; 
      fSampSize = stYUM.fSampSize;
      lDataSizeTotal = stYUM.lDataSizeTotal;
      pfTimeSectGain = stYUM.pfBFGain;    /* for .spec */
      dTNextBF = stYUM.dTNextBF;          /* for .spec */
      iNumChans = stYUM.iNumChans;
  
      /* some more user-input validation */
      if (iNumSubBands > 0)
      {
          if (iNumSubBands > iNumChans)
          {
              (void) fprintf(stderr,
                             "ERROR: Invalid number of sub-bands! Number of "
                             "sub-bands must be less than or equal to number of "
                             "channels!\n");
              PrintUsage(pcProgName);
              return YAPP_RET_ERROR;
          }
  
          /* there should be an integer number of sub-bands within the band */
          if ((iNumChans % iNumSubBands) != 0)
          {
              (void) fprintf(stderr,
                             "ERROR: Invalid number of sub-bands! Number of "
                             "sub-bands must be a factor of number of "
                             "channels!\n");
              PrintUsage(pcProgName);
              return YAPP_RET_ERROR;
          }
      }
      if (iDecFac != 0)
      {
          if (iNumChans % iDecFac != 0)
          {
              (void) fprintf(stderr,
                             "ERROR: Invalid decimation factor! Decimation factor "
                             "must be a factor of number of channels!\n");
              PrintUsage(pcProgName);
              return YAPP_RET_ERROR;
          }
  
          /* calculate the number of channels in the output */
          iOutNumChans = iNumChans / iDecFac;
      }
  
      /* allocate memory for the delay table */
      g_piOffsetTab = (int *) YAPP_Malloc((size_t) stYUM.iNumChans,
                                          sizeof(int),
                                          YAPP_FALSE);
      if (NULL == g_piOffsetTab)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s!\n",
                         strerror(errno));
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
      iRet = YAPP_CalcDelays(dDM, stYUM, fLaw, &g_piOffsetTab, &iMaxOffset);
      if (iRet != YAPP_RET_SUCCESS)
      {
          (void) fprintf(stderr,
                         "ERROR: Calculating delays failed!\n");
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
  
      if (0 == iDecFac)
      {
          /* calculate the corrected start time */
          dDelay = (double) -4.148741601e6
                   * (((double) 1.0 / pow(INFINITY, fLaw))
                      - ((double) 1.0 / pow(stYUM.fFMax, fLaw)))
                   * dDM;    /* in ms */
          iStartOffset = (int) (dDelay / stYUM.dTSamp);
          fStartOffset = iStartOffset * dTSampInSec;
      }
  
      /* ensure that the block size is at least equivalent to the maximum offset,
         because we don't read beyond the second buffer */
      if (iBlockSize < iMaxOffset)
      {
          (void) printf("WARNING: Block size is less than the calculated maximum"
                        " offset! Changing block size to %d.\n",
                        iMaxOffset);
          iBlockSize = iMaxOffset;
      }
  
      /* calculate bytes to skip and read */
      if (0.0 == dDataProcTime)
      {
          dDataProcTime = (iTimeSamps * dTSampInSec) - dDataSkipTime;
      }
      /* check if the input time duration is less than the length of the
         data */
      else if (dDataProcTime > (iTimeSamps * dTSampInSec))
      {
          (void) fprintf(stderr,
                         "WARNING: Input time is longer than length of "
                         "data!\n");
      }
  
      lBytesToSkip = (long) floor(dDataSkipTime / dTSampInSec)
                                                      /* number of samples */
                            * iNumChans
                            * fSampSize;
      if (fSampSize >= 2.0)
      {
          /* will be true for fSampSize = 0.5 and 1.0 */
          assert(0 == lBytesToSkip % (int) fSampSize);
      }
      lBytesToProc = (long) floor(dDataProcTime / dTSampInSec)
                                                      /* number of samples */
                            * iNumChans
                            * fSampSize;
      if (fSampSize >= 2.0)
      {
          /* will be true for fSampSize = 0.5 and 1.0 */
          assert(0 == lBytesToProc % (int) fSampSize);
      }
  
      if (lBytesToSkip >= lDataSizeTotal)
      {
          (void) fprintf(stderr,
                         "ERROR: Data to be skipped is greater than or equal to "
                         "the size of the file!\n");
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
  
      if ((lBytesToSkip + lBytesToProc) > lDataSizeTotal)
      {
          (void) printf("WARNING: Total data to be read (skipped and processed) "
                        "is more than the size of the file! ");
          lBytesToProc = lDataSizeTotal - lBytesToSkip;
          dDataProcTime = ((double) lBytesToProc * dTSampInSec)
                          / (iNumChans * fSampSize);
          (void) printf("Newly calculated size of data to be processed: %ld "
                        "bytes\n",
                        lBytesToProc);
      }
  
      if (iBlockSize == iMaxOffset)
      {
          if (lBytesToProc < (iBlockSize * iNumChans * fSampSize))
          {
              /* if the block size is equivalent to the maximum delay that is to
                 be applied, and if the number of bytes to be processed is less
                 than the block size, de-dispersion will be affected, as we don't
                 have more than two block buffers. if both conditions are true,
                 force the number of bytes to be processed to be equivalent to
                 the block size/maximum offset */
              (void) printf("WARNING: Amount of data to be processed is less "
                            "than the calculated maximum offset! Will process "
                            "more data than what was requested.\n");
              lBytesToProc = iBlockSize * iNumChans * fSampSize;
          }
      }
      else
      {
          /* here, iBlockSize > iMaxOffset */
          assert(iBlockSize > iMaxOffset);
  
          if (lBytesToProc < (iMaxOffset * iNumChans * fSampSize))
          {
              /* if the number of bytes to be processed is less than the maximum
                 offset, de-dispersion will be affected, as we don't have more
                 than two block buffers. in this situation, force the number of
                 bytes to be processed to be equivalent to the maximum offset */
              (void) printf("WARNING: Amount of data to be processed is less "
                            "than the calculated maximum offset! Will process "
                            "more data than what was requested.\n");
              lBytesToProc = iMaxOffset * iNumChans * fSampSize;
              (void) printf("WARNING: Amount of data to be processed is less "
                            "than the block size! Adjusting block size "
                            "accordingly.\n");
              iBlockSize = iMaxOffset;
          }
          else
          {
              if (lBytesToProc < (iBlockSize * iNumChans * fSampSize))
              {
                  /* here, iMaxOffset <=(eqv) lBytesToProc <(eqv) iBlockSize */
                  (void) printf("WARNING: Amount of data to be processed is "
                                "less than the block size! Adjusting block size "
                                "accordingly.\n");
                  iBlockSize = lBytesToProc / (iNumChans * fSampSize);
              }
          }
      }
  
      /* since we may have adjusted the number of bytes to be processed, correct
         the number of bytes to be skipped, too */
      if ((lBytesToSkip + lBytesToProc) > lDataSizeTotal)
      {
          (void) printf("WARNING: Total data to be read (skipped and processed) "
                        "is more than the size of the file! ");
          lBytesToSkip = lDataSizeTotal - lBytesToProc;
          (void) printf("Newly calculated size of data to be skipped: %ld "
                        "bytes\n",
                        lBytesToSkip);
      }
  
      iTimeSampsToSkip = (int) (lBytesToSkip / (iNumChans * fSampSize));
      (void) printf("Skipping\n"
                    "    %ld of %ld bytes\n"
                    "    %d of %d time samples\n"
                    "    %.10g of %.10g seconds\n",
                    lBytesToSkip,
                    lDataSizeTotal,
                    iTimeSampsToSkip,
                    stYUM.iTimeSamps,
                    (iTimeSampsToSkip * dTSampInSec),
                    (stYUM.iTimeSamps * dTSampInSec));
  
      iTimeSampsToProc = (int) (lBytesToProc / (iNumChans * fSampSize));
      iNumReads = (int) ceilf((float) iTimeSampsToProc / iBlockSize);
  
      /* optimisation - store some commonly used values in variables */
      iTotSampsPerBlock = iNumChans * iBlockSize;
  
      (void) printf("Processing\n"
                    "    %ld of %ld bytes\n"
                    "    %d of %d time samples\n"
                    "    %.10g of %.10g seconds\n"
                    "in %d reads with block size %d time samples...\n",
                    lBytesToProc,
                    lDataSizeTotal,
                    iTimeSampsToProc,
                    iTimeSamps,
                    (iTimeSampsToProc * dTSampInSec),
                    (iTimeSamps * dTSampInSec),
                    iNumReads,
                    iBlockSize);
  
      /* allocate memory for the time sample goodness flag array */
      g_pcIsTimeGood = (char *) YAPP_Malloc((size_t) iTimeSampsToProc,
                                            sizeof(char),
                                            YAPP_FALSE);
      if (NULL == g_pcIsTimeGood)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s!\n",
                         strerror(errno));
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
      /* set all elements to 'YAPP_TRUE' */
      (void) memset(g_pcIsTimeGood, YAPP_TRUE, iTimeSampsToProc);
  
      /* add header for .fil or .tim file format */
      if ((YAPP_FORMAT_FIL == iOutputFormat)
          || (YAPP_FORMAT_DTS_TIM == iOutputFormat))
      {
          /* populate the YUM structure */
          (void) memcpy(&stYUMOut, &stYUM, sizeof(YUM_t));
      }
  
      /* kludge: simulate sub-band dedispersion by marking other channels as bad
         so that they don't go into the summation */
      if (iNumSubBands > 0)
      {
          /* compute the number of channels per sub-band. iNumChans is a multiple
             of iNumSubBands */
          iChansPerSubBand = iNumChans / iNumSubBands;
          if (stYUM.cIsBandFlipped)
          {
              iStartChan = ((iNumSubBands - iSubBand) * iChansPerSubBand) - 1;
              iEndChan = (iNumSubBands - (iSubBand + 1)) * iChansPerSubBand;
              for (i = 0; i < iNumChans; ++i)
              {
                  if (!((i >= iEndChan) && (i < iStartChan)))
                  {
                      stYUM.pcIsChanGood[i] = YAPP_FALSE;
                  }
              }
          }
          else
          {
              iStartChan = iSubBand * iChansPerSubBand;
              iEndChan = (iSubBand + 1) * iChansPerSubBand;
              for (i = 0; i < iNumChans; ++i)
              {
                  if (!((i >= iStartChan) && (i < iEndChan)))
                  {
                      stYUM.pcIsChanGood[i] = YAPP_FALSE;
                  }
              }
          }
          stYUMOut.iNumGoodChans = iChansPerSubBand;
      }
  
      fStatBW = stYUMOut.iNumGoodChans * fChanBW;  /* in MHz */
      (void) printf("Usable bandwidth                  : %g MHz\n", fStatBW);
      fNoiseRMS = 1.0 / sqrt(fStatBW * dTSamp * 1e3);
      (void) printf("Expected noise RMS                : %g\n", fNoiseRMS);
  
      /* open the data file for reading */
      g_pFData = fopen(pcFileSpec, "r");
      if (NULL == g_pFData)
      {
          (void) fprintf(stderr,
                         "ERROR: Opening file %s failed! %s.\n",
                         pcFileSpec,
                         strerror(errno));
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
  
      /* allocate memory for the primary and secondary buffers, based on the
         number of channels and time samples */
      g_pfBuf0 = (float *) YAPP_Malloc((size_t) iNumChans * iBlockSize,
                                       sizeof(float),
                                       YAPP_FALSE);
      if (NULL == g_pfBuf0)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s!\n",
                         strerror(errno));
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
      g_pfBuf1 = (float *) YAPP_Malloc((size_t) iNumChans * iBlockSize,
                                       sizeof(float),
                                       YAPP_FALSE);
      if (NULL == g_pfBuf1)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s!\n",
                         strerror(errno));
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
  
      if (1 == iNumReads)
      {
          cIsLastBlock = YAPP_TRUE;
      }
  
      /* allocate memory for storing the dedispersed data */
      g_pfDedispData = (float *) YAPP_Malloc((size_t) iBlockSize,
                                             sizeof(float),
                                             YAPP_FALSE);
      if (NULL == g_pfDedispData)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s!\n",
                         strerror(errno));
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
  
      if (iDecFac != 0)
      {
          /* allocate memory for the primary and secondary buffers, based on the
             number of channels and time samples */
          g_pfDecDedispData = (float *) YAPP_Malloc((size_t) iOutNumChans
                                                    * iBlockSize,
                                                    sizeof(float),
                                                    YAPP_FALSE);
          if (NULL == g_pfDecDedispData)
          {
              (void) fprintf(stderr,
                             "ERROR: Memory allocation failed! %s!\n",
                             strerror(errno));
              YAPP_CleanUp();
              return YAPP_RET_ERROR;
          }
      }
  
      if (YAPP_FORMAT_FIL == iFormat)
      {
          /* TODO: Need to do this only if the file contains the header */
          /* skip the header */
          (void) fseek(g_pFData, (long) stYUM.iHeaderLen, SEEK_SET);
          /* skip data, if any are to be skipped */
          (void) fseek(g_pFData, lBytesToSkip, SEEK_CUR);
      }
      else
      {
          /* skip data, if any are to be skipped */
          (void) fseek(g_pFData, lBytesToSkip, SEEK_SET);
      }
  
      /* read the first block of data */
      (void) printf("Reading data block %d.", iReadBlockCount);
      (void) fflush(stdout);
      iReadItems = YAPP_ReadData(g_pFData,
                                 g_pfBuf0,
                                 fSampSize,
                                 iTotSampsPerBlock);
      if (YAPP_RET_ERROR == iReadItems)
      {
          (void) fprintf(stderr, "ERROR: Reading data failed!\n");
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
      pfPriBuf = g_pfBuf0;
      pfSpectrum = g_pfBuf0;
      iPrimaryBuf = BUF_0;
      --iNumReads;
      ++iReadBlockCount;
  
      /* calculate the number of time samples in the block - this may not
         be iBlockSize for the last block, and should be iBlockSize for
         all other blocks */
      iNumSamps = iReadItems / iNumChans;
  
      if (YAPP_FORMAT_SPEC == iFormat)
      {
          /* flag bad time sections, and if required, normalise within the beam
             flip time section and perform gain correction */
          for (i = 0; i < iNumSamps; ++i)
          {
              if ((dTNow >= (*stYUM.padBadTimes)[iBadTimeSect][BADTIME_BEG])
                  && (dTNow <= (*stYUM.padBadTimes)[iBadTimeSect][BADTIME_END]))
              {
                  cIsInBadTimeRange = YAPP_TRUE;
                  g_pcIsTimeGood[((iReadBlockCount-1)*iBlockSize)+i]
                      = YAPP_FALSE;
              }
  
              if ((YAPP_TRUE == cIsInBadTimeRange)
                  && (dTNow > (*stYUM.padBadTimes)[iBadTimeSect][BADTIME_END]))
              {
                  cIsInBadTimeRange = YAPP_FALSE;
                  ++iBadTimeSect;
              }
  
              /* get the beam flip time section corresponding to this
                 sample */
              if (dTNow > dTNextBF)
              {
                  dTNextBF += stYUM.dTBFInt;
                  ++iTimeSect;
                  if (iTimeSect >= stYUM.iBFTimeSects)
                  {
                      (void) fprintf(stderr,
                                     "ERROR: Beam flip time section anomaly "
                                     "detected!\n");
                      YAPP_CleanUp();
                      return YAPP_RET_ERROR;
                  }
              }
  
              pfTimeSectGain = stYUM.pfBFGain + (iTimeSect * iNumChans);
              pfSpectrum = pfPriBuf + i * iNumChans;
              for (j = 0; j < iNumChans; ++j)
              {
                  if (stYUM.pcIsChanGood[j])
                  {
                      pfSpectrum[j] = (pfSpectrum[j]
                                       / stYUM.pfBFTimeSectMean[iTimeSect])
                                      - pfTimeSectGain[j];
                  }
                  else    /* remove bad channels (just for plotting; bad channels
                             are not considered in dedispersion) */
                  {
                      pfSpectrum[j] = 0.0;
                  }
              }
  
              dTNow += dTSampInSec;   /* in s */
          }
      }
  
      if (cHasGraphics)
      {
          /* open the PGPLOT graphics device */
          if (cIsPlotToFile)
          {
              /* build the name of the PGPLOT device */
              pcFilename = YAPP_GetFilenameFromPath(pcFileSpec);
              (void) strcpy(acDev, pcFilename);
              (void) strcat(acDev, PLOT_DDS_SUFFIX);
              (void) strcat(acDev, EXT_PS);
              (void) strcat(acDev, PG_DEV_PS);
  
              iRet = cpgopen(acDev);
              if (iRet <= 0)
              {
                  (void) fprintf(stderr,
                                 "ERROR: Opening graphics device %s failed!\n",
                                 acDev);
                  YAPP_CleanUp();
                  return YAPP_RET_ERROR;
              }
              cpgask(YAPP_FALSE);
          }
          else
          {
              iRet = cpgopen(PG_DEV);
              if (iRet <= 0)
              {
                  (void) fprintf(stderr,
                                 "ERROR: Opening graphics device %s failed!\n",
                                 PG_DEV);
                  YAPP_CleanUp();
                  return YAPP_RET_ERROR;
              }
              cpgask(YAPP_TRUE);
          }
  
         /* set the background colour to white and the foreground colour to
            black, if user requires so */
         if (YAPP_TRUE == iInvCols)
         {
             cpgscr(0, 1.0, 1.0, 1.0);
             cpgscr(1, 0.0, 0.0, 0.0);
         }
  
          /* set up the plot's X-axis */
          g_pfXAxis = (float *) YAPP_Malloc((size_t) iBlockSize,
                                            sizeof(float),
                                            YAPP_FALSE);
          if (NULL == g_pfXAxis)
          {
              (void) fprintf(stderr,
                             "ERROR: Memory allocation failed! %s!\n",
                             strerror(errno));
              cpgclos();
              YAPP_CleanUp();
              return YAPP_RET_ERROR;
          }
  
          /* set up the image plot's Y-axis (frequency) */
          g_pfYAxis = (float *) YAPP_Malloc((size_t) iNumChans,
                                            sizeof(float),
                                            YAPP_FALSE);
          if (NULL == g_pfYAxis)
          {
              (void) fprintf(stderr,
                             "ERROR: Memory allocation failed! %s!\n",
                             strerror(errno));
              cpgclos();
              YAPP_CleanUp();
              return YAPP_RET_ERROR;
          }
          if (stYUM.cIsBandFlipped)
          {
              for (i = 0; i < iNumChans; ++i)
              {
                  g_pfYAxis[i] = stYUM.fFMax - i * fChanBW;
              }
          }
          else
          {
              for (i = 0; i < iNumChans; ++i)
              {
                  g_pfYAxis[i] = stYUM.fFMin + i * fChanBW;
              }
          }
  
          /* allocate memory for the cpgimag() plotting buffer */
          g_pfPlotBuf = (float *) YAPP_Malloc((size_t) iNumChans * iBlockSize,
                                              sizeof(float),
                                              YAPP_FALSE);
          if (NULL == g_pfPlotBuf)
          {
              (void) fprintf(stderr,
                             "ERROR: Memory allocation failed! %s!\n",
                             strerror(errno));
              cpgclos();
              YAPP_CleanUp();
              return YAPP_RET_ERROR;
          }
      }
  
      /* generate dedispersed data file name and config file name from the input
         file name */
      pcFilename = YAPP_GetFilenameFromPath(pcFileSpec);
      if (0 == iNumSubBands)  /* not doing sub-band dedispersion */
      {
          (void) sprintf(acFileDedisp,
                         "%s.%s%g",
                         pcFilename,
                         INFIX_DEDISPERSE,
                         dDM);
      }
      else
      {
          (void) sprintf(acFileDedisp,
                         "%s.%s%g.%s%d",
                         pcFilename,
                         INFIX_DEDISPERSE,
                         dDM,
                         INFIX_SUBBAND,
                         iSubBand);
      }
      if (YAPP_FORMAT_DTS_TIM == iOutputFormat)
      {
          (void) strcat(acFileDedisp, EXT_TIM);
      }
      else if (YAPP_FORMAT_DTS_DDS == iOutputFormat)
      {
          (void) strcat(acFileDedisp, EXT_DEDISPSPEC);
      }
      else    /* if filterbank */
      {
          (void) strcat(acFileDedisp, EXT_FIL);
      }
  
      /* update header for .fil or .tim file format */
      if ((YAPP_FORMAT_FIL == iOutputFormat)
          || (YAPP_FORMAT_DTS_TIM == iOutputFormat))
      {
          stYUMOut.iNumBits = YAPP_SAMPSIZE_32;
          /* enter the start time corrected for dispersion */
          stYUMOut.dTStart = stYUM.dTStart - (fStartOffset / 86400);
          stYUMOut.dDM = dDM;
          if (YAPP_FORMAT_DTS_TIM == iOutputFormat)
          {
              stYUMOut.cIsBandFlipped = YAPP_FALSE;
          }
  
          /* if doing sub-band dedispersion, update the number of channels,
             bandwidth, etc. */
          if (iNumSubBands > 0)
          {
              stYUMOut.iNumChans = iChansPerSubBand;
              stYUMOut.fBW = stYUMOut.iNumChans * stYUMOut.fChanBW;
              if (stYUM.cIsBandFlipped)
              {
                  stYUMOut.fFMin = stYUM.fFMax
                                   - (iStartChan * stYUMOut.fChanBW);
              }
              else
              {
                  stYUMOut.fFMin = stYUM.fFMin
                                   + (iStartChan * stYUMOut.fChanBW);
              }
              stYUMOut.fFMax = stYUMOut.fFMin
                              + ((stYUMOut.iNumChans - 1) * stYUMOut.fChanBW);
              if (0 == (stYUMOut.iNumChans % 2))   /* even number of channels */
              {
                  stYUMOut.fFCentre = (stYUMOut.fFMin - (stYUMOut.fChanBW / 2))
                                     + ((stYUMOut.iNumChans / 2)
                                        * stYUMOut.fChanBW);
              }
              else                                /* odd number of channels */
              {
                  stYUMOut.fFCentre = stYUMOut.fFMin
                                     + (((float) stYUMOut.iNumChans / 2)
                                        * stYUMOut.fChanBW);
              }
          }
  
          if (iDecFac != 0)
          {
              stYUMOut.iNumChans = iOutNumChans;
              stYUMOut.fChanBW = stYUM.fChanBW * iDecFac;
          }
  
          /* write metadata to disk */
          iRet = YAPP_WriteMetadata(acFileDedisp, iOutputFormat, stYUMOut);
          if (iRet != YAPP_RET_SUCCESS)
          {
              (void) fprintf(stderr,
                             "ERROR: Writing metadata failed for file %s!\n",
                             acFileDedisp);
              YAPP_CleanUp();
              return YAPP_RET_ERROR;
          }
      }
  
      if (iDecFac != 0)
      {
          /* allocate memory for the output delay table */
          g_piOutOffsetTab = (int *) YAPP_Malloc((size_t) iOutNumChans,
                                                 sizeof(int),
                                                 YAPP_FALSE);
          if (NULL == g_piOutOffsetTab)
          {
              (void) fprintf(stderr,
                             "ERROR: Memory allocation failed! %s!\n",
                             strerror(errno));
              YAPP_CleanUp();
              return YAPP_RET_ERROR;
          }
          iRet = YAPP_CalcDelays(dDM, stYUMOut, fLaw, &g_piOutOffsetTab, &iOutMaxOffset);
          if (iRet != YAPP_RET_SUCCESS)
          {
              (void) fprintf(stderr,
                             "ERROR: Calculating delays failed!\n");
              YAPP_CleanUp();
              return YAPP_RET_ERROR;
          }
      }
  
      /* open the output file for appending data */
      pFDedispData = fopen(acFileDedisp, "a");
      if (NULL == pFDedispData)
      {
          fprintf(stderr,
                  "ERROR: Opening file %s failed! %s.\n",
                  acFileDedisp,
                  strerror(errno));
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
  
      /* set up the plots */
      if (cHasGraphics)
      {
          cpgsubp(1, 3);
          cpgsch(PG_CH_3P);
      }
  
      /* dedisperse the data */
      while (iNumReads >= 0)
      {
          /* for optimisation - calculate ((iReadBlockCount - 1) * iBlockSize) */
          iReadSmpCount = (iReadBlockCount - 1) * iBlockSize;
  
          if (cHasGraphics)
          {
              for (i = 0; i < iBlockSize; ++i)
              {
                  g_pfXAxis[i] = (iReadSmpCount * dTSampInSec)
                                 + (i * dTSampInSec);
              }
  
              cpgpanl(1, 1);
              /* erase just before plotting, to reduce flicker */
              cpgeras();
  
              pfSpectrum = pfPriBuf;
              fDataMin = pfSpectrum[0];
              fDataMax = pfSpectrum[0];
              for (j = 0; j < iBlockSize; ++j)
              {
                  pfSpectrum = pfPriBuf + j * iNumChans;
                  for (k = 0; k < iNumChans; ++k)
                  {
                      if (pfSpectrum[k] < fDataMin)
                      {
                          fDataMin = pfSpectrum[k];
                      }
                      if (pfSpectrum[k] > fDataMax)
                      {
                          fDataMax = pfSpectrum[k];
                      }
                  }
              }
  
              #ifdef DEBUG
              (void) printf("Minimum value of data             : %g\n",
                            fDataMin);
              (void) printf("Maximum value of data             : %g\n",
                            fDataMax);
              #endif
  
              /* get the transpose of the two-dimensional array */
              k = 0;
              l = 0;
              m = 0;
              for (i = 0; i < iBlockSize; ++i)
              {
                  pfSpectrum = pfPriBuf + i * iNumChans;
                  for (j = 0; j < iNumChans; ++j)
                  {
                      l = m + k * iBlockSize;
                      g_pfPlotBuf[l] = pfSpectrum[j];
                      ++k;
                  }
                  k = 0;
                  ++m;
              }
  
              Plot2D(g_pfPlotBuf, fDataMin, fDataMax,
                     g_pfXAxis, iBlockSize, dTSampInSec,
                     g_pfYAxis, iNumChans, fChanBW,
                     "Time - Start Time (s)", "Frequency (MHz)",
                     "Before Dedispersion",
                     iColourMap);
          }
  
          /* read the next data block, for reading samples that would be
             shifted in */
          if (!(cIsLastBlock))
          {
              (void) printf("\rReading data block %d.", iReadBlockCount);
              (void) fflush(stdout);
              if (BUF_0 == iPrimaryBuf)
              {
                  iReadItems = YAPP_ReadData(g_pFData,
                                             g_pfBuf1,
                                             fSampSize,
                                             iTotSampsPerBlock);
                  pfSecBuf = g_pfBuf1;
              }
              else
              {
                  iReadItems = YAPP_ReadData(g_pFData,
                                             g_pfBuf0,
                                             fSampSize,
                                             iTotSampsPerBlock);
                  pfSecBuf = g_pfBuf0;
              }
              if (ferror(g_pFData))
              {
                  (void) fprintf(stderr, "ERROR: File read failed!\n");
                  if (cHasGraphics)
                  {
                      cpgclos();
                  }
                  (void) fclose(pFDedispData);
                  YAPP_CleanUp();
                  return YAPP_RET_ERROR;
              }
              if (iReadItems < iTotSampsPerBlock)
              {
                  iDiff = iTotSampsPerBlock - iReadItems;
  
                  /* reset remaining elements to '\0' */
                  if (BUF_0 == iPrimaryBuf)
                  {
                      (void) memset((g_pfBuf1 + iReadItems),
                                    '\0',
                                    (sizeof(float) * iDiff));
                  }
                  else
                  {
                      (void) memset((g_pfBuf0 + iReadItems),
                                    '\0',
                                    (sizeof(float) * iDiff));
                  }
              }
  
              /* calculate the number of time samples in the block - this may not
                 be iBlockSize for the last block, and should be iBlockSize for
                 all other blocks */
              iNumSamps = iReadItems / iNumChans;
  
              if (YAPP_FORMAT_SPEC == iFormat)
              {
                  /* flag bad time sections, and if required, normalise within
                     the beam flip time section and perform gain correction */
                  for (i = 0; i < iNumSamps; ++i)
                  {
                      if ((dTNow >= (*stYUM.padBadTimes)[iBadTimeSect][BADTIME_BEG])
                          && (dTNow
                              <= (*stYUM.padBadTimes)[iBadTimeSect][BADTIME_END]))
                      {
                          cIsInBadTimeRange = YAPP_TRUE;
                          g_pcIsTimeGood[((iReadBlockCount-1)*iBlockSize)+i]
                              = YAPP_FALSE;
                      }
  
                      if ((YAPP_TRUE == cIsInBadTimeRange)
                          && (dTNow > (*stYUM.padBadTimes)[iBadTimeSect][BADTIME_END]))
                      {
                          cIsInBadTimeRange = YAPP_FALSE;
                          ++iBadTimeSect;
                      }
  
                      /* get the beam flip time section corresponding to this
                         sample */
                      if (dTNow > dTNextBF)
                      {
                          dTNextBF += stYUM.dTBFInt;
                          ++iTimeSect;
                          if (iTimeSect >= stYUM.iBFTimeSects)
                          {
                              (void) fprintf(stderr,
                                             "ERROR: Beam flip time section "
                                             "anomaly detected!\n");
                              if (cHasGraphics)
                              {
                                  cpgclos();
                              }
                              (void) fclose(pFDedispData);
                              YAPP_CleanUp();
                              return YAPP_RET_ERROR;
                          }
                      }
  
                      pfTimeSectGain = stYUM.pfBFGain + (iTimeSect * iNumChans);
                      pfSpectrum = pfSecBuf + i * iNumChans;
                      for (j = 0; j < iNumChans; ++j)
                      {
                          if (stYUM.pcIsChanGood[j])
                          {
                              pfSpectrum[j] = (pfSpectrum[j]
                                               / stYUM.pfBFTimeSectMean[iTimeSect])
                                              - pfTimeSectGain[j];
                          }
                          else    /* remove bad channels (just for plotting; bad
                                     channels are not considered in
                                     dedispersion) */
                          {
                              pfSpectrum[j] = 0.0;
                          }
                      }
                  }
  
                  dTNow += dTSampInSec;   /* in s */
              }
          }
  
          /* clear the g_pfDedispData array */
          (void) memset(g_pfDedispData,
                        '\0',
                        (sizeof(float) * iBlockSize));
  
          if (0 == iDecFac)
          {
              for (k = 0; k < iBlockSize; ++k)
              {
                  pfSpectrum = pfPriBuf + k * iNumChans;
                  for (l = 0; l < iNumChans; ++l)
                  {
                      if (stYUM.pcIsChanGood[l])
                      {
                          /* get the offset for the corresponding DM and frequency
                             channel from the offset table */
                          iOffset = g_piOffsetTab[l];
                          /* apply the delay - shift all time samples up */
                          if ((k + iOffset) >= iBlockSize)
                          {
                              if (!(cIsLastBlock))
                              {
                                  m = k + iOffset - iBlockSize;
                                  pfOffsetSpec = pfSecBuf + m * iNumChans;
                                  pfSpectrum[l] = pfOffsetSpec[l];
                                  if (g_pcIsTimeGood[iReadSmpCount+k+iOffset])
                                  {
                                      g_pfDedispData[k] += pfSpectrum[l];
                                      ++iEffcNumGoodChans;
                                  }
                              }
                          }
                          else
                          {
                              pfOffsetSpec = pfPriBuf
                                             + (k + iOffset) * iNumChans;
                              pfSpectrum[l] = pfOffsetSpec[l];
                              if (g_pcIsTimeGood[iReadSmpCount+k+iOffset])
                              {
                                  g_pfDedispData[k] += pfSpectrum[l];
                                  ++iEffcNumGoodChans;
                              }
                          }
                      }
                  }
  
                  /* get the average over all the good channels */
                  /*if (iEffcNumGoodChans != 0)
                  {
                      g_pfDedispData[k] /= iEffcNumGoodChans;
                  }
                  else
                  {
                      g_pfDedispData[k] = 0.0;
                  }
  
                  g_pfDedispData[k] /= fNoiseRMS;*/
  
                  /* reset the effective number of good channels */
                  iEffcNumGoodChans = 0;
              }
          }
          else
          {
              for (k = 0; k < iBlockSize; ++k)
              {
                  pfSpectrum = pfPriBuf + k * iNumChans;
                  pfOutSpec = g_pfDecDedispData + k * iOutNumChans;
                  for (l = 0; l < iNumChans; ++l)
                  {
                      n = l / iDecFac;
                      if (((float) l / iDecFac) - (l /iDecFac) == 0)
                      {
                          g_pfDedispData[k] = 0.0;
                      }
                      if (stYUM.pcIsChanGood[l])
                      {
                          /* get the offset for the corresponding DM and frequency
                             channel from the offset table */
                          iOffset = g_piOffsetTab[l] - g_piOutOffsetTab[n];
                          /* apply the delay - shift all time samples up */
                          if ((k + iOffset) >= iBlockSize)
                          {
                              if (!(cIsLastBlock))
                              {
                                  m = k + iOffset - iBlockSize;
                                  pfOffsetSpec = pfSecBuf + m * iNumChans;
                                  pfSpectrum[l] = pfOffsetSpec[l];
                                  if (g_pcIsTimeGood[iReadSmpCount+k+iOffset])
                                  {
                                      g_pfDedispData[k] += pfSpectrum[l];
                                      pfOutSpec[n] = g_pfDedispData[k];
                                      ++iEffcNumGoodChans;
                                  }
                              }
                          }
                          else
                          {
                              pfOffsetSpec = pfPriBuf
                                             + (k + iOffset) * iNumChans;
                              pfSpectrum[l] = pfOffsetSpec[l];
                              if (g_pcIsTimeGood[iReadSmpCount+k+iOffset])
                              {
                                  g_pfDedispData[k] += pfSpectrum[l];
                                  pfOutSpec[n] = g_pfDedispData[k];
                                  ++iEffcNumGoodChans;
                              }
                          }
                      }
                  }
  
                  /* get the average over all the good channels */
                  /*if (iEffcNumGoodChans != 0)
                  {
                      g_pfDedispData[k] /= iEffcNumGoodChans;
                  }
                  else
                  {
                      g_pfDedispData[k] = 0.0;
                  }
  
                  g_pfDedispData[k] /= fNoiseRMS;*/
  
                  /* reset the effective number of good channels */
                  iEffcNumGoodChans = 0;
              }
          }
  
          if (cHasGraphics)
          {
              for (i = 0; i < iBlockSize; ++i)
              {
                  g_pfXAxis[i] = (iReadSmpCount * dTSampInSec)
                                 + (i * dTSampInSec) - fStartOffset;
              }
  
              cpgpanl(1, 2);
              /* erase just before plotting, to reduce flicker */
              cpgeras();
  
              pfSpectrum = pfPriBuf;
              fDataMin = pfSpectrum[0];
              fDataMax = pfSpectrum[0];
              for (j = 0; j < iBlockSize; ++j)
              {
                  pfSpectrum = pfPriBuf + j * iNumChans;
                  for (k = 0; k < iNumChans; ++k)
                  {
                      if (pfSpectrum[k] < fDataMin)
                      {
                          fDataMin = pfSpectrum[k];
                      }
                      if (pfSpectrum[k] > fDataMax)
                      {
                          fDataMax = pfSpectrum[k];
                      }
                  }
              }
  
              #ifdef DEBUG
              (void) printf("Minimum value of data             : %g\n",
                            fDataMin);
              (void) printf("Maximum value of data             : %g\n",
                            fDataMax);
              #endif
  
              /* get the transpose of the two-dimensional array */
              k = 0;
              l = 0;
              m = 0;
              for (i = 0; i < iBlockSize; ++i)
              {
                  pfSpectrum = pfPriBuf + i * iNumChans;
                  for (j = 0; j < iNumChans; ++j)
                  {
                      l = m + k * iBlockSize;
                      g_pfPlotBuf[l] = pfSpectrum[j];
                      ++k;
                  }
                  k = 0;
                  ++m;
              }
  
              Plot2D(g_pfPlotBuf, fDataMin, fDataMax,
                     g_pfXAxis, iBlockSize, dTSampInSec,
                     g_pfYAxis, iNumChans, fChanBW,
                     "Time - Start Time (s)",
                     "Frequency (MHz)",
                     "After Delay Correction",
                     iColourMap);
          }
  
          /* pfPriBuf now contains dedispersed, non-collapsed data. write this to
             disk if required */
          if (YAPP_FORMAT_FIL == iOutputFormat)
          {
              if (0 == iDecFac)
              {
                  /* kludgy fix to the last block problem */
                  if (!cIsLastBlock)
                  {
                      (void) fwrite(pfPriBuf,
                                    sizeof(float),
                                    iNumChans * iBlockSize,
                                    pFDedispData);
                  }
                  else
                  {
                      (void) fwrite(pfPriBuf,
                                    sizeof(float),
                                    iNumChans * (iBlockSize - iMaxOffset),
                                    pFDedispData);
                  }
              }
              else
              {
                  /* kludgy fix to the last block problem */
                  if (!cIsLastBlock)
                  {
                      (void) fwrite(g_pfDecDedispData,
                                    sizeof(float),
                                    iOutNumChans * iBlockSize,
                                    pFDedispData);
                  }
                  else
                  {
                      (void) fwrite(g_pfDecDedispData,
                                    sizeof(float),
                                    iOutNumChans * (iBlockSize - iMaxOffset),
                                    pFDedispData);
                  }
  
              }
          }
          else
          {
              /* kludgy fix to the last block problem */
              if (!cIsLastBlock)
              {
                  (void) fwrite(g_pfDedispData,
                                sizeof(float),
                                iBlockSize,
                                pFDedispData);
              }
              else
              {
                  (void) fwrite(g_pfDedispData,
                                sizeof(float),
                                iBlockSize - iMaxOffset,
                                pFDedispData);
              }
          }
  
          if (cHasGraphics)
          {
              cpgpanl(1, 3);
              /* erase just before plotting, to reduce flicker */
              cpgeras();
  
              fDataMin = g_pfDedispData[0];
              fDataMax = g_pfDedispData[0];
              for (l = 0; l < iBlockSize; ++l)
              {
                  if (g_pfDedispData[l] < fDataMin)
                  {
                      fDataMin = g_pfDedispData[l];
                  }
                  if (g_pfDedispData[l] > fDataMax)
                  {
                      fDataMax = g_pfDedispData[l];
                  }
              }
  
              #ifdef DEBUG
              (void) printf("Minimum value of data             : %g\n",
                            fDataMin);
              (void) printf("Maximum value of data             : %g\n",
                            fDataMax);
              #endif
  
              cpgsvp(PG_2D_VP_ML, PG_2D_VP_MR, PG_2D_VP_MB, PG_2D_VP_MT);
  
              cpgswin(g_pfXAxis[0], g_pfXAxis[iBlockSize-1], fDataMin, fDataMax);
              cpgbox("BCNST", 0.0, 0, "BCNST", 0.0, 0);
              cpglab("Time - Start Time (s)",
                     "Total Power (arbitrary units)",
                     "Dedispersed Time Series");
              cpgsci(PG_CI_PLOT);
              cpgline(iBlockSize, g_pfXAxis, g_pfDedispData);
              cpgsci(PG_CI_DEF);
  
              if (!(cIsLastBlock))
              {
                  if (!(cIsNonInteractive))
                  {
                      /* draw the 'next' and 'exit' buttons */
                      cpgsvp(PG_VP_BUT_ML,
                             PG_VP_BUT_MR,
                             PG_VP_BUT_MB,
                             PG_VP_BUT_MT);
                      cpgswin(PG_BUT_L, PG_BUT_R, PG_BUT_B, PG_BUT_T);
                      cpgsci(PG_BUT_FILLCOL); /* set the fill colour */
                      cpgrect(PG_BUTNEXT_L,
                              PG_BUTNEXT_R,
                              PG_BUTNEXT_B,
                              PG_BUTNEXT_T);
                      cpgrect(PG_BUTEXIT_L,
                              PG_BUTEXIT_R,
                              PG_BUTEXIT_B,
                              PG_BUTEXIT_T);
                      cpgsci(0);  /* set colour index to white */
                      cpgtext(PG_BUTNEXT_TEXT_L, PG_BUTNEXT_TEXT_B, "Next");
                      cpgtext(PG_BUTEXIT_TEXT_L, PG_BUTEXIT_TEXT_B, "Exit");
  
                      fButX = (PG_BUTNEXT_R - PG_BUTNEXT_L) / 2;
                      fButY = (PG_BUTNEXT_T - PG_BUTNEXT_B) / 2;
  
                      while (YAPP_TRUE)
                      {
                          iRet = cpgcurs(&fButX, &fButY, &cCurChar);
                          if (0 == iRet)
                          {
                              (void) fprintf(stderr,
                                             "WARNING: Reading cursor "
                                             "parameters failed!\n");
                              break;
                          }
  
                          if (((fButX >= PG_BUTNEXT_L)
                               && (fButX <= PG_BUTNEXT_R))
                              && ((fButY >= PG_BUTNEXT_B)
                                  && (fButY <= PG_BUTNEXT_T)))
                          {
                              /* animate button click */
                              cpgsci(PG_BUT_FILLCOL);
                              cpgtext(PG_BUTNEXT_TEXT_L,
                                      PG_BUTNEXT_TEXT_B,
                                      "Next");
                              cpgsci(0);  /* set colour index to white */
                              cpgtext(PG_BUTNEXT_CL_TEXT_L,
                                      PG_BUTNEXT_CL_TEXT_B,
                                      "Next");
                              (void) usleep(PG_BUT_CL_SLEEP);
                              cpgsci(PG_BUT_FILLCOL); /* set colour index to fill
                                                         colour */
                              cpgtext(PG_BUTNEXT_CL_TEXT_L,
                                      PG_BUTNEXT_CL_TEXT_B,
                                      "Next");
                              cpgsci(0);  /* set colour index to white */
                              cpgtext(PG_BUTNEXT_TEXT_L,
                                      PG_BUTNEXT_TEXT_B,
                                      "Next");
                              cpgsci(1);  /* reset colour index to black */
                              (void) usleep(PG_BUT_CL_SLEEP);
  
                              break;
                          }
                          else if (((fButX >= PG_BUTEXIT_L)
                                    && (fButX <= PG_BUTEXIT_R))
                              && ((fButY >= PG_BUTEXIT_B)
                                  && (fButY <= PG_BUTEXIT_T)))
                          {
                              /* animate button click */
                              cpgsci(PG_BUT_FILLCOL);
                              cpgtext(PG_BUTEXIT_TEXT_L,
                                      PG_BUTEXIT_TEXT_B,
                                      "Exit");
                              cpgsci(0);  /* set colour index to white */
                              cpgtext(PG_BUTEXIT_CL_TEXT_L,
                                      PG_BUTEXIT_CL_TEXT_B,
                                      "Exit");
                              (void) usleep(PG_BUT_CL_SLEEP);
                              cpgsci(PG_BUT_FILLCOL); /* set colour index to fill
                                                         colour */
                              cpgtext(PG_BUTEXIT_CL_TEXT_L,
                                      PG_BUTEXIT_CL_TEXT_B,
                                      "Exit");
                              cpgsci(0);  /* set colour index to white */
                              cpgtext(PG_BUTEXIT_TEXT_L,
                                      PG_BUTEXIT_TEXT_B,
                                      "Exit");
                              cpgsci(1);  /* reset colour index to black */
                              (void) usleep(PG_BUT_CL_SLEEP);
  
                              cpgclos();
                              (void) fclose(pFDedispData);
                              YAPP_CleanUp();
                              return YAPP_RET_SUCCESS;
                          }
                      }
                  }
                  else
                  {
                      /* pause before erasing */
                      (void) usleep(PG_PLOT_SLEEP);
                  }
              }
          }
  
          if (BUF_0 == iPrimaryBuf)
          {
              iPrimaryBuf = BUF_1;
              pfPriBuf = g_pfBuf1;
          }
          else
          {
              iPrimaryBuf = BUF_0;
              pfPriBuf = g_pfBuf0;
          }
  
          --iNumReads;
          ++iReadBlockCount;
          if (0 == iNumReads)
          {
              cIsLastBlock = YAPP_TRUE;
          }
      }
  
      (void) printf("DONE!\n");
  
      if (cHasGraphics)
      {
          cpgclos();
      }
  
      (void) fclose(pFDedispData);
      YAPP_CleanUp();
  
      return YAPP_RET_SUCCESS;
  }
  
  /*
   * Prints usage information
   */
  void PrintUsage(const char *pcProgName)
  {
      (void) printf("Usage: %s [options] <data-file>\n",
                    pcProgName);
      (void) printf("    -h  --help                          ");
      (void) printf("Display this usage information\n");
      (void) printf("    -s  --skip <time>                   ");
      (void) printf("The length of data in seconds, to be\n");
      (void) printf("                                        ");
      (void) printf("skipped\n");
      (void) printf("    -p  --proc <time>                   ");
      (void) printf("The length of data in seconds, to be\n");
      (void) printf("                                        ");
      (void) printf("processed\n");
      (void) printf("                                        ");
      (void) printf("(default is all)\n");
      (void) printf("    -n  --nsamp <samples>               ");
      (void) printf("Number of samples read in one block\n");
      (void) printf("                                        ");
      (void) printf("(default is 4096 samples)\n");
      (void) printf("    -d  --dm <dm>                       ");
      (void) printf("DM at which to de-disperse\n");
      (void) printf("                                        ");
      (void) printf("(default is 10.0)\n");
      (void) printf("    -l  --law <law>                     ");
      (void) printf("Dispersion law\n");
      (void) printf("                                        ");
      (void) printf("(default is 2.0)\n");
      (void) printf("    -b  --nsubband <nsubband>           ");
      (void) printf("Number of sub-bands\n");
      (void) printf("                                        ");
      (void) printf("(must be < number of channels)\n");
      (void) printf("    -u  --subband <subband>             ");
      (void) printf("Sub-band to dedisperse\n");
      (void) printf("                                        ");
      (void) printf("(within [0, nsubband))\n");
      (void) printf("    -o  --out-format <format>           ");
      (void) printf("Output format - 'dds', 'tim', or 'fil'\n");
      (void) printf("                                        ");
      (void) printf("(default is 'tim')\n");
      (void) printf("    -f  --f-decimate <factor>           ");
      (void) printf("Decimation-in-frequency factor,\n");
      (void) printf("                                        ");
      (void) printf("for 'fil' output\n");
      (void) printf("    -g  --graphics                      ");
      (void) printf("Turn on plotting\n");
      (void) printf("    -m  --colour-map <name>             ");
      (void) printf("Colour map for plotting\n");
      (void) printf("                                        ");
      (void) printf("(default is 'jet')\n");
      (void) printf("    -i  --invert                        ");
      (void) printf("Invert background and foreground\n");
      (void) printf("                                        ");
      (void) printf("colours in plots\n");
      (void) printf("    -e  --non-interactive               ");
      (void) printf("Run in non-interactive mode\n");
      (void) printf("    -v  --version                       ");
      (void) printf("Display the version\n");
  
      return;
  }
  
  int YAPP_CalcDelays(double dDM,
                      YUM_t stYUM,
                      float fLaw,
                      int** ppiOffsetTab,
                      int* piMaxOffset)
  {
      int i = 0;
      float fF1 = 0.0;
      float fF2 = 0.0;
      double dDelay = 0.0;
      float fFMaxCalc = stYUM.fFMax;     /* reference frequency */
  
      /* calculate quadratic delays */
      /* NOTE: delay may not be 0 for the highest frequency channel,
         but the offset samples may be (depending on the sampling rate) */
  #ifdef DEBUG
      {
          FILE *pFFileDelaysQuad = NULL;
  
          pFFileDelaysQuad = fopen(YAPP_FILE_DELAYS_QUAD, "w");
          if (NULL == pFFileDelaysQuad)
          {
              fprintf(stderr,
                      "ERROR: Opening file %s failed! %s.\n",
                      YAPP_FILE_DELAYS_QUAD,
                      strerror(errno));
              YAPP_CleanUp();
              return YAPP_RET_ERROR;
          }
  #endif
      if (dDM < 0)
      {
          if (stYUM.cIsBandFlipped)
          {
              fF1 = fFMaxCalc;
              fF2 = stYUM.fFMax;
              for (i = stYUM.iNumChans - 1; i >= 0; --i)
              {
                  dDelay = (double) -4.148741601e6
                           * (((double) 1.0 / pow(fF1, fLaw))
                              - ((double) 1.0 / pow(fF2, fLaw)))
                           * dDM;    /* in ms */
                  (*ppiOffsetTab)[i] = (int) (dDelay / stYUM.dTSamp);
  #ifdef DEBUG
                  (void) fprintf(pFFileDelaysQuad,
                                 "%d %g %d\n",
                                 i,
                                 dDelay,
                                 (*ppiOffsetTab)[i]);
  #endif
                  fF2 -= stYUM.fChanBW;
              }
              *piMaxOffset = (*ppiOffsetTab)[0];
          }
          else
          {
              fF1 = fFMaxCalc;
              fF2 = stYUM.fFMax;
              for (i = 0; i < stYUM.iNumChans; ++i)
              {
                  dDelay = (double) -4.148741601e6
                           * (((double) 1.0 / pow(fF1, fLaw))
                              - ((double) 1.0 / pow(fF2, fLaw)))
                           * dDM;    /* in ms */
                  (*ppiOffsetTab)[i] = (int) (dDelay / stYUM.dTSamp);
  #ifdef DEBUG
                  (void) fprintf(pFFileDelaysQuad,
                                 "%d %g %d\n",
                                 i,
                                 dDelay,
                                 (*ppiOffsetTab)[i]);
  #endif
                  fF2 -= stYUM.fChanBW;
              }
              *piMaxOffset = (*ppiOffsetTab)[stYUM.iNumChans-1];
          }
      }
      else
      {
          if (stYUM.cIsBandFlipped)
          {
              fF1 = fFMaxCalc;
              fF2 = stYUM.fFMax;
              for (i = 0; i < stYUM.iNumChans; ++i)
              {
                  dDelay = (double) -4.148741601e6
                           * (((double) 1.0 / pow(fF1, fLaw))
                              - ((double) 1.0 / pow(fF2, fLaw)))
                           * dDM;    /* in ms */
                  (*ppiOffsetTab)[i] = (int) (dDelay / stYUM.dTSamp);
  #ifdef DEBUG
                  (void) fprintf(pFFileDelaysQuad,
                                 "%d %g %d\n",
                                 i,
                                 dDelay,
                                 (*ppiOffsetTab)[i]);
  #endif
                  fF2 -= stYUM.fChanBW;
              }
              *piMaxOffset = (*ppiOffsetTab)[stYUM.iNumChans-1];
          }
          else
          {
              fF1 = fFMaxCalc;
              fF2 = stYUM.fFMax;
              for (i = stYUM.iNumChans - 1; i >= 0; --i)
              {
                  dDelay = (double) -4.148741601e6
                           * (((double) 1.0 / pow(fF1, fLaw))
                              - ((double) 1.0 / pow(fF2, fLaw)))
                           * dDM;    /* in ms */
                  (*ppiOffsetTab)[i] = (int) (dDelay / stYUM.dTSamp);
  #ifdef DEBUG
                  (void) fprintf(pFFileDelaysQuad,
                                 "%d %g %d\n",
                                 i,
                                 dDelay,
                                 (*ppiOffsetTab)[i]);
  #endif
                  fF2 -= stYUM.fChanBW;
              }
              *piMaxOffset = (*ppiOffsetTab)[0];
          }
      }
  #ifdef DEBUG
          (void) fclose(pFFileDelaysQuad);
      }
  #endif
  
      return YAPP_RET_SUCCESS;
  }
  
  