line.c

Go to the documentation of this file.

/* 
Copyright 2021 Arizona Board of Regents on behalf of the University of Arizona.

For commercial uses, to obtain a license to sell and or sublicense copies of
the Software please contact the University of Arizona at Tech Launch Arizona:
lewish@tla.arizona.edu.

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction for non-commercial use educational and
research purposes, including without limitation the rights to use, copy,
modify, merge, publish, distribute, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:

The entire above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

1. Additional Required Provisions

1.1.  Arbitration. The parties agree that if a dispute arises between them
concerning this Agreement, the parties may be required to submit the matter to
arbitration pursuant to Arizona law.

1.2.  Applicable Law and Venue. This Agreement shall be interpreted pursuant
to the laws of the State of Arizona. Any arbitration or litigation between the
Parties shall be conducted in Pima County, ARIZONA, and LICENSEE hereby
submits to venue and jurisdiction in Pima County, ARIZONA.

1.3.  Non-Discrimination. The Parties agree to be bound by state and federal
laws and regulations governing equal opportunity and non-discrimination and
immigration.

1.4.   Appropriation of Funds. The Parties recognize that performance by
ARIZONA may depend upon appropriation of funds by the State Legislature of
ARIZONA. If the Legislature fails to appropriate the necessary funds, or if
ARIZONA’S appropriation is reduced during the fiscal year, ARIZONA may cancel
this Agreement without further duty or obligation. ARIZONA will notify
LICENSEE as soon as reasonably possible after it knows of the loss of funds.

1.5.   Conflict of Interest. This Agreement is subject to the provisions of
A.R.S. 38-511 and other conflict of interest regulations. Within three years
of the EFFECTIVE DATE, ARIZONA may cancel this Agreement if any person
significantly involved in initiating, negotiating, drafting, securing, or
creating this Agreement for or on behalf of ARIZONA becomes an employee or
consultant in any capacity of LICENSEE with respect to the subject matter of
this Agreement.

 */


#include "extern.h"

#include <gsl/gsl_statistics_double.h>

/*  Changes:
 *
 *  Added nbox function: twf; March 01, 2016
 *
 */

double workingTsys    = 0.0;
double workingTcal    = 0.0;
double workingTime    = 0.0;
float workingX[MAX_CHANS];
float workingY[MAX_CHANS];

char rmsResultsStr[256];

#define USE_GSL 1

int computeRms(y, n, rms)
float *y, *rms;
int n;
{
  int i, verbose=0;
#ifdef USE_GSL
//#warning Using the GSL version of computeRMS()
  double meangsl, rmsgsl, tss, gsv;
  double d[MAX_CHANS];
#else
  double mean, sums;
#endif

  if(isSet(PF_V_FLAG))
  {
    st_report("computeRms(): Verbose Set");
    verbose = 1;
  }

  clearRMS();

  if(n==0)
  {
    return(1);
  }

#ifdef USE_GSL

  for(i=0;i<n;i++)
  {
    d[i] = y[i];
  }

  meangsl = gsl_stats_mean(d, 1, n);

  if(verbose)
  {
    st_print("GSL Mean      = %12.6f mK\n", meangsl);
  }

  tss = gsl_stats_tss_m(d, 1, n, meangsl);

  gsv = gsl_stats_variance_m(d, 1, n, meangsl);

  if(verbose)
  {
    st_print("GSL Variance  = %12.6f mK\n", gsv       * 1000.0);
    st_print("GSL TSS       = %12.6f mK\n", tss       * 1000.0);
    st_print("GSL sqrt(TTS) = %12.6f mK\n", sqrt(tss) * 1000.0);
  }

  rmsgsl  = gsl_stats_absdev_m(&d[0], 1, n, meangsl);
  if(verbose)
  {
    st_print("GSL Abs Dev   = %12.6f mK\n", rmsgsl * 1000.0);
  }

  rmsgsl  = gsl_stats_sd_m(d, 1, n, meangsl);

  if(verbose)
  {
    st_print("GSL rms       = %12.6f mK\n", rmsgsl * 1000.0);
  }

  if(verbose)
  {
    st_print("GSL Skew      = %12.6f mK\n", gsl_stats_skew(d, 1, n) * 1000.0);
    st_print("GSL Kurtosis  = %12.6f mK\n", gsl_stats_kurtosis(d, 1, n) * 1000.0);
    st_print("GSL AutoCorr  = %12.6f mK\n", gsl_stats_lag1_autocorrelation(d, 1, n) * 1000.0);
  }

  /* Use the gsl version */
  *rms = rmsgsl;

#else

  mean = 0.0;
  for(i=0;i<n;i++)
  {
    mean += y[i];
  }

  mean /= (float)n;

//  st_print("My Mean       = %f\n", mean);

  sums = 0.0;
  for(i=0;i<n;i++)
  {
    sums += ((y[i] - mean) * (y[i] - mean));
  }

  sums /= (float)n;

  if(sums > 0.0)
  {
    *rms = sqrt(sums);
  }
  else
  {
    *rms = 100.0;
  }

//  st_print("My rms        = %f\n", *rms * 1000.0);

#endif

  return(0);
}



int do_specCombClear()
{
//  st_print("Resetting Accumulators\n");

  bzero((char *)&accum_image, sizeof(struct LP_IMAGE_INFO));
  strcpy(accum_image.name, "ACCUM");

  bzero((char *)&workingX, sizeof(workingX));
  bzero((char *)&workingY, sizeof(workingY));

  workingTsys = 0.0;
  workingTcal = 0.0;
  workingTime = 0.0;

  return(0);
}


int line_accum(image)
struct LP_IMAGE_INFO *image;
{
  int k;
  float *i, iwt, *o, owt;
  double raoff, decoff;
  static double xoff, yoff, roff, doff, noint;

  if(env->accumreset)
  {
    do_specCombClear();

    clearFBOffets();

    strcpy(accum_image.magic, "LP_IMAGE_INFO");
    
    bcopy((char *)&image->head, (char *)&accum_image.head, sizeof(accum_image.head));

    accum_image.ainfo.bscan      = (float)image->head.scan;

    xoff  = image->head.az_offset;
    yoff  = image->head.el_offset;

    roff  = image->head.epocra  - image->head.xsource;
    doff  = image->head.epocdec - image->head.ysource;

    noint = image->head.noint;

    env->avgEl = 0.0;
  }

  if(image->head.noint != noint)
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: NOINT Differs: %.1f, %.1f\n", 
                        image->head.scan, image->head.noint, noint);
    }

    incrementRejNoint();

    return(1);
  }

    /* Has user supplied an upper limit on acceptable tsys? */
  if((env->tsyslimit > 0.0 && image->head.stsys > env->tsyslimit) || gsl_isnan(image->head.stsys))
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: TSYS %.1f >  TSLimit: %.1f\n", 
                        image->head.scan, image->head.stsys, env->tsyslimit);
    }

    incrementBogusTsys();

    return(1);
  }

  if(!env->ignoremapoff && image->head.az_offset != xoff)
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: AZ_OFFSET Differs: %.1f, %.1f\n", 
                        image->head.scan, image->head.az_offset*3600.0, xoff*3600.0);
    }

    incrementRejOffsets();

    return(1);
  }

  if(!env->ignoremapoff && image->head.el_offset != yoff)
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: EL_OFFSET Differs: %.1f, %.1f\n", 
                        image->head.scan, image->head.el_offset*3600.0, yoff*3600.0);
    }

    incrementRejOffsets();

    return(1);
  }

  raoff  = image->head.epocra  - image->head.xsource;
  decoff = image->head.epocdec - image->head.ysource;

  if(!env->ignoremapoff && raoff != roff)
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: RA_OFFSET Differs: %.1f, %.1f\n", 
                        image->head.scan, raoff*3600.0, roff*3600.0);
    }

    incrementRejOffsets();

    return(1);
  }

  if(!env->ignoremapoff && decoff != doff)
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: DEC_OFFSET Differs: %.1f, %.1f\n", 
                        image->head.scan, decoff*3600.0, doff*3600.0);
    }

    incrementRejOffsets();

    return(1);
  }

  if(lp_fcmp(accum_image.head.freqres, image->head.freqres, 0.01 ))
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: Acuum Freq Res %.3f vs. Scan Freq Res %.3f, differ\n", 
                        image->head.scan, accum_image.head.freqres, image->head.freqres);
    }

    incrementRejFreqres();

    return(1);
  }

  if(fabs(accum_image.head.epocra - image->head.epocra ) > 0.001 && !env->ignorepos)
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: Ra %.3f & %.3f, differ\n", 
                        image->head.scan, accum_image.head.epocra, image->head.epocra);
    }

    incrementRejRa();


    return(1);
  }

  if(fabs(accum_image.head.epocdec - image->head.epocdec) > 0.001 && !env->ignorepos)
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: Dec %.3f & %.3f, differ\n", 
                        image->head.scan, accum_image.head.epocdec, image->head.epocdec);
    }

    incrementRejDec();

    return(1);
  }

  if((int)accum_image.head.sideband != (int)image->head.sideband )
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: Sb %d & %d, differ\n", 
                        image->head.scan, (int)accum_image.head.sideband, (int)image->head.sideband);
    }

    incrementRejSb();

    return(1);
  }

  if(fabs(accum_image.head.restfreq / 1000.0 - image->head.restfreq / 1000.0) > env->dfreq )
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: Rest Freq %.5f & %.5f, differ\n", image->head.scan,
                                        accum_image.head.restfreq / 1000.0, 
                                        image->head.restfreq / 1000.0);
    }

    incrementRejFreq();

    return(1);
  }

  if(env->dvel == 0.0) /* Check for uninitialize value */
  {
    env->dvel = 0.001;
  }

  if(fabs(accum_image.head.refpt_vel - image->head.refpt_vel) > env->dvel  && !env->ignorepos)
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("line_accum(%.2f), Rejected: Vel %.3f & %.3f, differ\n", 
                        image->head.scan, accum_image.head.refpt_vel, image->head.refpt_vel);
    }

    incrementRejVel();

    return(1);
  }

  /* Do a weighted average base on integration_time / tsys^2 */

  accum_image.ainfo.escan   = (float)image->head.scan;

  workingTime += image->head.inttime;
  workingTsys += (image->head.stsys * image->head.inttime);
  workingTcal += (image->head.tcal  * image->head.inttime);

  iwt = image->head.inttime / ( image->head.stsys * image->head.stsys );
  owt = accum_image.ainfo.wt;

  i = image->yplot1;
  o = workingY;

  for(k=0;k<accum_image.head.noint;k++, i++, o++)
  {
    *o = (*o * owt + *i * iwt) / (owt + iwt);
  }

  accum_image.ainfo.wt   += iwt;
  accum_image.ainfo.escan = (float)image->head.scan;
  accum_image.ainfo.scanct++;

  env->avgEl += image->head.el;

  printScan(&image->head, 0);

  env->accumreset = 0; /* Don't reset again until told to */

  return(0);
}


int line_accumSelected()
{
  int i, j = 0, indx, cols, rows;
  struct HEADER *p;
  double ret;
  char tmp[80], tbuf[256];

  zero_bl_averages();

  env->accumreset = 1; /* Reset accumalators */

  getTerminalSize(&cols, &rows);

  rows -= 2;

  resetRejections();

  freeX();
  freeY();

  /* Clear the titleBlock which may contain out-dated info */
  sprintf(tbuf, "titleblock /c");
  parseCmd(tbuf);

//  if(env->markers == 0)
//  {
//    env->markers = 16;
//    st_report("Setting Markers to Default: %d", env->markers);
//  }

  st_report("Spectral Line Accum");

  if(env->autobad)
  {
    st_report("Flagged Bad Channels will be removed");
  }

  if(!env->numofscans)
  {
    st_warn("No Valid Data Files Loaded");

    env->accumreset = 1;

    return(1);
  }

  if(env->accumreset)
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      if(env->freq > 0.0)
      {
        st_report("Looking for scans using Freq = %.3f +/- %3.3f MHz", env->freq, env->dfreq*1000.0);
      }
      else
      {
        st_report("Looking for scans using ALL Freqs");
      }

      if(env->check_date)
      {
        st_report("Looking for scans taken on %.4f", env->date);
      }
      else
      {
        st_report("Looking for scans taken on ALL dates");
      }

      if(env->check_source)
      {
        st_report("Looking for scans on %s", env->source);
      }
      else
      {
        st_report("Looking for scans on ALL Sources");
      }

      if(env->check_scans)
      {
        st_report("Looking for scans starting at: %.2f up to %.2f", env->bscan, env->escan);
      }
      else
      {
        st_report("Looking for Entire Range of scan numbers");
      }

      switch(env->bkmatch)
      {
        case BKMATCH_FEED:
            st_report("Looking for scans ending in %.2f", env->if_feed);
          break;

        case BKMATCH_RESOLUTION:
            st_report("Looking for scans with backend resolution: %.3f", env->bkres);
          break;

        case BKMATCH_NAME:
            st_report("Looking for scans with backend name: %s", env->backend);
          break;

        case BKMATCH_NAMERES:
            st_report("Looking for scans with backend name: %s & backend resolution: %.3f", env->backend, env->bkres);
          break;
      }
    }

    env->dobaseline = 0;
  }

  p    = &sa[0];
  indx = 0;

  for(i=0;i<env->numofscans;i++,p++)
  {
    if(verboseLevel & VERBOSE_CRITIRIA)
    {
      printf("Passing Scan %.2f to meetsCritiria(): ", p->scan);
    }

    /* Check the most obvious critiria */
    if(!meetsCritiria(p, CHECK_ALL))
    {
      continue;
    }

    /* Check for proper LINE* Mode */
    if(strncmp(p->obsmode, "LINEBSP",  7) && 
       strncmp(p->obsmode, "LINEPS",   6) &&
       strncmp(p->obsmode, "LINEFS",   6) &&
       strncmp(p->obsmode, "LINETPON", 8) &&
       strncmp(p->obsmode, "LINEAPS",  7)) /* We only process BSP, APS & PS so far; treat FS the same */
    {
      incrementRejMode();
      continue; /* Not Spectral */
    }

    if(p->openpar[5] == 5555.5) /* a BSP focal */
    {
      continue; /* Ignore */
    }

    /* Is this scan ignored??? */
    if(isScanIgnored(p->scan))
    {
      incrementRejIgnored();
      continue;
    }

    if(!env->allow_five_pt && isSpecFivePoint(p))
    {
      incrementRejFive();
      continue;
    }

    /* Made it this far; read it in */
    ret = read_scan(p->scan, env->filename, h0->yplot1, &h0->head, env->autobad);
    env->curScan = p->scan;

    if(ret > 0.0)
    {
      if(!(j%rows))
      {
        st_print("%s   Res    Xoff   Yoff\n", listHeader);
        j = 1;
      }

      /* Go accum this scan */
      if(!line_accum(h0))
      {
        indx++; /* Accum went well */
        j++;
      }
    }
  }

  if(indx == 0)
  {
    st_print("No Scans Found Matching Search Criteria\n");

    env->last_scan_count = 0;

    printRejections();

    return(1);
  }

  env->avgEl /= (double)indx;

  /* Normalize the accum's */

  accum_image.head.stsys   = workingTsys / workingTime;
  accum_image.head.tcal    = workingTcal / workingTime;
  accum_image.head.inttime = workingTime;
  accum_image.head.effint  = accum_image.head.inttime;

  accum_image.binfo.doBaseLine = env->dobaseline;
  bcopy((char *)&env->baselines, (char *)&accum_image.binfo.baselines, sizeof(accum_image.binfo.baselines));

  st_report(" ");

  st_print("%s Avg Tsys: %.1f for %d scans; ", prompt,
                accum_image.head.stsys,
                accum_image.ainfo.scanct);

  env->last_scan_count = accum_image.ainfo.scanct;

  if(accum_image.head.inttime / 60.0 > 120)
  {
    st_print("%.2f Hours Int Time\n", accum_image.head.inttime / 3600.0);
  }
  else
  {
    st_print("%.2f Minutes Int Time\n", accum_image.head.inttime / 60.0);
  }

  /* now copy the workingY into the image struct */

  for(i=0;i<accum_image.head.noint;i++)
  {
    accum_image.yplot1[i] = workingY[i];
  }

  line_removeBadChans(&accum_image, 1);

  strcpy(tmp, h0->name);  /* Save the name first */
  bcopy((char *)&accum_image, (char *)h0, sizeof(struct LP_IMAGE_INFO));
  strcpy(h0->name, tmp);   /* Restore the name */

  h0->narray = 1;
  h0->nsize  = h0->head.noint * sizeof(float);

  h0->head.nostac = env->last_scan_count;
  h0->head.fscan  = h0->ainfo.bscan;
  h0->head.lscan  = h0->ainfo.escan;
  h0->head.rms    = env->rms;

  linePlot(h0);

  env->accumreset = 1;

  return(0);
}




int line_determine_min_max(image)
struct LP_IMAGE_INFO *image;
{
  int i;
  float *f;
                                                   /* redetermine min and max */
  image->datamax = -9999999999.9;
  image->datamin = 9999999999.0;

  f = &image->yplot1[env->bdrop];
  for(i=env->bdrop;i<image->head.noint-env->edrop;i++,f++)
  {
    if(*f > image->datamax)
    {
      image->datamax  = *f;
      image->peakChan = env->peakChan = i;
    }

    if(*f < image->datamin)
    {
      image->datamin = *f;
    }
  }

  if(fabs(image->datamax - image->datamin) < DATASPREAD)
  {
    st_print("Scan #%7.2f is all %.4f\n", image->head.scan, image->datamin);

    image->datamin -= 1.0;
    image->datamax += 1.0;

    return(1);
  }

  env->datamin = image->datamin;
  env->datamax = image->datamax;

  return(0);
}


int processSpectra(image)
struct LP_IMAGE_INFO *image;
{

  return(0);
}


int line_display_spectra(image, plot)
struct LP_IMAGE_INFO *image;
int plot;
{
  if(env->procspec)
  {
    if(processSpectra(image))
    {
      return(1);
    }
  }

  line_removeBadChans(image, 1);

  image->narray = 1;

  if(isSet(PF_F_FLAG)) /* Spec focus contains several plots */
  {
    image->head.noint = image->head.datalen / sizeof(float);
  }

  image->nsize  = image->head.noint * sizeof(float);

  linePlot(image);

  return(0);
}



int removeFittedData(image)
struct LP_IMAGE_INFO *image;
{
  int i;
  float *f, *g;

  if(env->rm_fit_data && image->pinfo.bplot)
  {
    f = &image->yplot1[0];
    g = &image->bplot[0];

    for(i=0;i<image->head.noint;i++,f++,g++)
    {
      *f = *f - *g;
    }

    env->rm_fit_data = image->pinfo.bplot = 0;

    image->yplot1[0] = 0.0; /* Most times the 0'th entry cant be removed */
  }

  return(0);
}


  /* Draw the plot */
int linePlot(image)
struct LP_IMAGE_INFO *image;
{
  if(image->head.noint == 0.0)
  {
    st_fail("Nothing in Array %s to Display", image->name);
    return(0);
  }

  st_report("Displaying Plot: %s", image->name);

  copyEnv2Image(image);

  /* Go fill in x-axis values */
  computeXaxis(image, 0);

  line_determine_min_max(image);

/* Check if we are running external viewer */
  if(xgraphic)
  {
    sendImageStruct(image, env->mode);

    validate(BSP);

    return(0);
  }
  else
  {
    st_report("No connection to xgraphic to send plot");
  }

  return(0);
}


int do_F_Command()
{
  float s;

  if(strlen(args->secondArg))
  {
    s = atof(args->secondArg);

    s = floor(s);

    subArgs(2, "%.2f", s+0.01);

    getScan(1);

    showScan();
  }

  return(0);
}


int do_S_Command()
{
  float s;

  if(strlen(args->secondArg))
  {
    s = atof(args->secondArg);

    s = floor(s);

    if(env->issmt)
    {
      subArgs(2, "%.2f", s+0.02);
    }
    else
    {
      subArgs(2, "%.2f", s+0.03);
    }

    getScan(1);

    showScan();
  }

  return(0);
}


int setHWDC(image)
struct LP_IMAGE_INFO *image;
{

  if(!strncmp(image->head.backend, "ARO", 3))
  {
    env->hwdc = 1.0;
  }
  else
  if(!strncmp(image->head.backend, "MAC", 3)  || !strncmp(image->head.backend, "HS_", 3))
  {
    env->hwdc = 1.0;
  }
  else
  if(!strncmp(image->head.backend, "FFB", 3))
  {
    env->hwdc = 1.0;
  }
  else
  {
    env->hwdc = 0.9;
  }

  st_report("Auto Setting HWDC to: %.3f", env->hwdc);

  return(0);

}


/* Unipops: 2.0*stsys/(abs(freqres)*1.0e6*inttime)^0.5
 */
int printRMS()
{
  int ret;
  double sq;

  struct LP_IMAGE_INFO *image = h0;

  /* Go compute the baseline which computes 
     the rms of the baseline region */
  ret = computeBaseLine(image);

  if(env->hwdc < 0.1)
  {
    setHWDC(image);
  }

  sq = sqrt(fabs(image->head.freqres) * 1e6 * image->head.inttime * env->hwdc);

  if(sq > 0.0)
  {
    env->trms = image->head.stsys * 2.0 / sq;

    env->p2p  = env->rms / sin(45.0 * CTR);
  }
  else
  {
    st_fail("sqrt(freqres * inttime * hwdc) = %f", sq);
    env->trms = 99999.0;
  }

  st_report("Using HWDC: %.3f", env->hwdc);


  sprintf(rmsResultsStr, "%.2f / %.2f mK over %d chans; Factor: %f", 
                                env->rms * 1000.0, env->trms * 1e3, ret,
                                env->rms / env->trms);

  st_report("Computed rms = %s", rmsResultsStr);

  if(isSet(PF_P_FLAG))
  {
    st_report("Computed P2P = %.2f mK", env->p2p * 1e3);
  }

  strcmprs(rmsResultsStr);

  return(0);
}


int selfRej()
{
  double ratio;

  if(strlen(args->secondArg))
  {
    ratio = atof(args->secondArg);

    if(env->rms / env->trms > ratio)
    {
      st_fail("Scan %.2f BAD", h0->head.scan);

      addIgnore(h0->head.scan);
    }
  }

  return(0);
}


int clearRMS()
{
  env->rms       = 0.0;
  env->slope     = 0.0;
  env->intercept = 0.0;

  return(0);
}



int line_showTable()
{
  int i, bc, ec, eFlag = 0, first=1;
  struct LP_IMAGE_INFO *image = h0;

  eFlag = isSet(PF_EXP_FLAG);;

  bc = env->bdrop;
  ec = image->head.noint - env->edrop;

  if(image->head.noint > 0)
  {
    for(i=bc;i<ec;i++)
    {
      if(first || !(i%8))
      {
        first = 0;
        st_print("\n%5d: ", i);
      }

      if(eFlag)
      {
        st_print("%14.6e ", image->yplot1[i]);
      }
      else
      {
        st_print("%10.6f ", image->yplot1[i]);
      }
    }

    st_print("\n");
  }
  else
  {
    st_fail("No data loaded");
  }

  return(0);
}


int line_MinMax(image)
struct LP_IMAGE_INFO *image;
{
  int i;
  float *f, min, max;

  max = -1.0e127; /* there are probably standard macros, but this works */
  min = 9.99e127;

  f = &image->yplot1[0];
  for(i=0;i<image->head.noint;i++)                                /* compute min and max */
  {
    if(f[i] > max)
    {
      max = f[i];
    }
    else
    if(f[i] < min)
    {
      min = f[i];
    }
  }

  image->datamin = min;
  image->datamax = max;

  st_report("Max = %f, Min = %f", image->datamax, image->datamin);

  return(0);
}



int doMeasure()
{
  int n;
  double from, to, diff, freeD, wireD;

  n = getNumericResponse("Enter Number of Cycles to Measure: ");

  _doTCur("Click the Left Mouse Button on First Peak");
  from = env->fcur;

  _doTCur("Click the Left Mouse Button on Last  Peak");
  to   = env->fcur;

  diff = fabs(from - to);

  diff *= 1e6;  /* Convert to Hz */

  diff /= (double)n; /* Reduce to one cycle */

  freeD = 1.0 / diff * SPEED_OF_LIGHT;
  wireD = 1.0 / diff * SPEED_OF_LIGHT_COAX;

  freeD /= 2.0;
  wireD /= 2.0;

  st_report("Spacing:               %7.3f MHz", diff / 1e6);
  st_report("Free Space Distance:   %7.3f Meters", freeD);
  st_report("Coaxial Wire Distance: %7.3f Meters", wireD);

  return(0);
}


int cmdShowAverage()
{
  if(args->ntokens > 1)
  {
    showAverage(args->secondArg);
  }
  else
  {
    showAverage(NULL);
  }

  return(0);
}

int showAverage(buf)
char *buf;
{
  int i, n, j=0;
  double avg = 0.0, ptmb = 0.0, rejs[4];
  struct LP_IMAGE_INFO *image = h0;

  for(i=0;i<4;i++)
  {
    rejs[i] = 1.0 / (1.0 + pow(10.0, -1.0 * env->rejections[i] / 10.0));
  }

  if(isPlanet(env->source))
  {
    st_report("Source is Planet");

    getPlanet();

    if(!strcasecmp(env->source, "Mercury"))
    {
      ptmb = env->mercury_tmb;
    }
    else
    if(!strcasecmp(env->source, "Venus"))
    {
      ptmb = env->venus_tmb;
    }
    else
    if(!strcasecmp(env->source, "Mars"))
    {
      ptmb = env->mars_tmb;
    }
    else
    if(!strcasecmp(env->source, "Jupiter"))
    {
      ptmb = env->jupiter_tmb;
    }
    else
    if(!strcasecmp(env->source, "Saturn"))
    {
      ptmb = env->saturn_tmb;
    }
    else
    if(!strcasecmp(env->source, "Uranus"))
    {
      ptmb = env->uranus_tmb;
    }
    else
    if(!strcasecmp(env->source, "Neptune"))
    {
      ptmb = env->neptune_tmb;
    }
  }

  if(env->mode == MODE_SPEC || env->mode == MODE_OTF) /* Examining Calibrate works the same */
  {
    if(isSet(PF_B_FLAG))  /* Use only baseline regions */
    {
      if(isSet(PF_S_FLAG))
      {
        j   = computeEachBaseLine(h0);
        avg = env->rst;
      }
      else
      {
        setFlag(PF_N_FLAG, 1);

        j = computeBaseLine(h0);

        setFlag(PF_N_FLAG, 0);

        avg = env->rst;
      }
    }

    if(j == 0) /* Because computeBaseLine() can return 0 if no baselines defined */
    {
      n = image->head.noint;

      for(i=env->bdrop;i<(n-env->edrop);i++)
      {
        avg += image->yplot1[i];
        j++;
      }

      avg /= (double)j;
    }

    if(buf)
    {
      st_report("Average of %s, using %d channels: %.6fK", buf, j, avg);
    }
    else
    {
      st_report("Average Bandpass Temperature, using %d channels: %.6fK", j, avg);
    }

    env->bp_avg = avg;

    if(ptmb > 0.0) /* It's a planet? */
    {
      st_report("Using Rejections Coeff: %f %f %f %f", rejs[0], rejs[1], rejs[2], rejs[3]);

      st_report("Tmb of %s:       %.2fK", env->source, ptmb);

      env->beam_eff = env->bp_avg / ptmb * 100.0 * rejs[0];

      st_report("Beam Eff for %s: %.2f%%", env->source, env->beam_eff);
    }
  }
  else
  if(env->mode == MODE_CONT)
  {
    if(ptmb > 0.0)
    {
      st_print("      Channel:      1          2          3          4\n");

      st_print("   Average Tr: %9.4f  %9.4f  %9.4f  %9.4f\n", 
                        local_ars.avgTr[0],
                        local_ars.avgTr[1],
                        local_ars.avgTr[2],
                        local_ars.avgTr[3]);

      st_print("  Planet Temp: %9.4f  %9.4f  %9.4f  %9.4f\n", 
                        ptmb, ptmb, ptmb, ptmb);

      st_print("Main Beam Eff: %9.2f  %9.2f  %9.2f  %9.2f\n", 
                        local_ars.avgTr[0] / ptmb * 100.0 * rejs[0],
                        local_ars.avgTr[1] / ptmb * 100.0 * rejs[1],
                        local_ars.avgTr[2] / ptmb * 100.0 * rejs[2],
                        local_ars.avgTr[3] / ptmb * 100.0 * rejs[3]);

      st_print("   Rejections: %9d  %9d  %9d  %9d\n", 
                                (int)env->rejections[0],
                                (int)env->rejections[1],
                                (int)env->rejections[2],
                                (int)env->rejections[3]);
                        
    }
  }

  return(0);
}



int doPshow()
{
  int nsize;
  double s, min=999999.9, max=-999999.9, min1, min2, max1, max2;
  double diff = 0.01;

  if(strlen(args->secondArg))
  {
    doPageCmd();

    clearAllImages();

    s = atof(args->secondArg);

    subArgs(2, "%.2f", s);
    _h0();                              /* Point to H0 */

    if(getScan(1))
    {
      return(1);
    }

    line_removeBadChans(h0, 1);

    if(isSet(PF_B_FLAG))
    {
      st_report("Removing Baselines");

      env->dobaseline = 1;

      line_doBaseLine(h0);
    }

    if(isSet(PF_I_FLAG))
    {
      doSpecInvert();
    }

    line_MinMax(h0);

    min1 = h0->datamin;
    max1 = h0->datamax;
    if(min1 < min)
    {
      min = min1;
    }
    if(max1 > max)
    {
      max = max1;
    }

    if(!strncmp(h0->head.backend, "FB_2", 4))
    {
      diff = 0.02; /* Filters in series inciments by 0.02 */
    }

    subArgs(2, "%.2f", s + diff);
    _h1();                              /* Point to H1 */

    if(getScan(1))
    {
      return(1);
    }

    line_removeBadChans(h0, 1);

    if(isSet(PF_B_FLAG))
    {
      st_report("Removing Baselines");

      env->dobaseline = 1;

      line_doBaseLine(h0);
    }

    if(isSet(PF_I_FLAG))
    {
      doSpecInvert();
    }

    line_MinMax(h0);

    min2 = h0->datamin;
    max2 = h0->datamax;
    if(min2 < min)
    {
      min = min2;
    }
    if(max2 > max)
    {
      max = max2;
    }

    nsize  = h0->head.noint * sizeof(float);

    env->ylower = min;
    env->yupper = max;
    env->yscale = Y_FIXED_SCALE;
    H0.nsize    = nsize;
    H0.narray   = 2;                    /* Let xgraphic know there are 2 arrays */

    bcopy((char *)&H1.yplot1, (char *)&H0.yplot2, nsize); /* Copy H1 array to H0's yplot2 */

    _h0();                              /* Point back to H0 */

    linePlot(h0);                       /* Plot it */
  }
  else
  {
    st_fail("Missing arg to pshow");
    return(1);
  }

  return(0);
}



int removeBadChans()
{
  if(env->mode == MODE_SPEC)
  {
    line_removeBadChans(h0, 1);
  }

  return(0);
}


int doQshow()
{
  int i, nsize;
  double s, minT=999999.9, maxT=-999999.9, min[4], max[4];
  double frac[] = { 0.00, 0.01, 0.02, 0.03 };

  if(strlen(args->secondArg))
  {
    doPageCmd();

    clearAllImages();

    s = atof(args->secondArg);

/* First */
    for(i=0;i<4;i++)
    {
      switch(i)
      {
        case 0: _h0();  /* Point to H0 */ break;
        case 1: _h1();  /* Point to H1 */ break;
        case 2: _h2();  /* Point to H2 */ break;
        case 3: _h3();  /* Point to H3 */ break;
      }

      subArgs(2, "%.2f", s + frac[i]);

      if(getScan(1))
      {
        return(1);
      }

      line_removeBadChans(h0, 1);

      if(isSet(PF_B_FLAG))
      {
        st_report("Removing Baselines");

        env->dobaseline = 1;

        line_doBaseLine(h0);
      }

      if(isSet(PF_I_FLAG))
      {
        doSpecInvert();
      }

      line_MinMax(h0);

      min[i] = h0->datamin;
      max[i] = h0->datamax;
      if(min[i] < minT)
      {
        minT = min[i];
      }
      if(max[i] > maxT)
      {
        maxT = max[i];
      }

      strncpy(H0.mhp[i].backend, h0->head.backend, 8);
      H0.mhp[i].backend[8] = '\0';

      H0.mhp[i].scan      = h0->head.scan;
      H0.mhp[i].restfreq  = h0->head.restfreq;
      H0.mhp[i].stsys     = h0->head.stsys;
      H0.mhp[i].tcal      = h0->head.tcal;
      H0.mhp[i].freqres   = h0->head.freqres;
      H0.mhp[i].sideband  = h0->head.sideband;
      H0.mhp[i].refpt     = h0->head.refpt;
      H0.mhp[i].refpt_vel = h0->head.refpt_vel;
      H0.mhp[i].deltax    = h0->head.deltax;
    }

    _h0();                              /* Point back to H0 */

    nsize  = h0->head.noint * sizeof(float);

    env->ylower = minT;
    env->yupper = maxT;
    env->yscale = Y_FIXED_SCALE;
    H0.nsize    = nsize;
    H0.narray   = 4;                    /* Let xgraphic know there are 2 arrays */

    bcopy((char *)&H1.yplot1, (char *)&H0.yplot2, nsize); /* Copy H1 array to H0's yplot2 */
    bcopy((char *)&H2.yplot1, (char *)&H0.yplot3, nsize); /* Copy H2 array to H0's yplot3 */
    bcopy((char *)&H3.yplot1, (char *)&H0.yplot4, nsize); /* Copy H3 array to H0's yplot4 */

    linePlot(h0);                       /* Plot it */
  }
  else
  {
    st_fail("Missing arg to pshow");

    return(1);
  }

  return(0);
}


int do_psc()
{
  int    i, save_bkmatch, save_ckfeed, req[2];
  int    nsize, first=1, rb = 0;
  double minT=999999.9, maxT=-999999.9, min[4], max[4];

  /* Check if user wants baselines removed */
  if(isSet(PF_B_FLAG))
  {
    rb = 1;
    st_report("Removing Baselines from Plots");
  }

  /* save environment */
  save_bkmatch    = env->bkmatch;
  save_ckfeed     = env->check_feed;

  env->bkmatch    = BKMATCH_FEED;
  env->check_feed = 1;

  if(!strncmp(args->currentCmd, "psc12", 5))
  {
    req[0] = 1;
    req[1] = 2;
  }
  else
  if(!strncmp(args->currentCmd, "psc34", 5))
  {
    req[0] = 3;
    req[1] = 4;
  }
  else
  if(!strncmp(args->currentCmd, "psc13", 5))
  {
    req[0] = 1;
    req[1] = 3;
  }
  else
  if(!strncmp(args->currentCmd, "psc24", 5))
  {
    req[0] = 2;
    req[1] = 4;
  }
  else
  {
    req[0] = 0;
    req[1] = 0;
  }

  _h1();  /* Point to H1 */
  for(i=0;i<4;i++)
  {
    switch(i)
    {
      case 0:
        if(i+1 == req[0] || i+1 == req[1])
        {
          showC1();

          if(rb)
          {
            doBaseline();
          }

          line_MinMax(h0);

          min[i] = h0->datamin;
          max[i] = h0->datamax;
          if(min[i] < minT)
          {
            minT = min[i];
          }
          if(max[i] > maxT)
          {
            maxT = max[i];
          }

          nsize = h0->head.noint * sizeof(float);

          if(first)
          {
            _copyImage(h0, &H0);

            _h2();  /* Point to H2 */

            first = 0;
          }
          else
          {
            bcopy((char *)&h0->yplot1, (char *)&H0.yplot2, nsize); /* Copy H1 array to H0's yplot2 */
          }
        }
        break;

      case 1:
        if(i+1 == req[0] || i+1 == req[1])
        {
          showC2();

          if(rb)
          {
            doBaseline();
          }

          line_MinMax(h0);

          min[i] = h0->datamin;
          max[i] = h0->datamax;
          if(min[i] < minT)
          {
            minT = min[i];
          }
          if(max[i] > maxT)
          {
            maxT = max[i];
          }

          nsize = h0->head.noint * sizeof(float);

          if(first)
          {
            _copyImage(h0, &H0);

            _h2();  /* Point to H2 */

            first = 0;
          }
          else
          {
            bcopy((char *)&h0->yplot1, (char *)&H0.yplot2, nsize); /* Copy H1 array to H0's yplot2 */
          }
        }
        break;

      case 2:
        if(i+1 == req[0] || i+1 == req[1])
        {
          showC3();

          if(rb)
          {
            doBaseline();
          }

          line_MinMax(h0);

          min[i] = h0->datamin;
          max[i] = h0->datamax;
          if(min[i] < minT)
          {
            minT = min[i];
          }
          if(max[i] > maxT)
          {
            maxT = max[i];
          }

          nsize = h0->head.noint * sizeof(float);

          if(first)
          {
            _copyImage(h0, &H0);

            _h2();  /* Point to H2 */

            first = 0;
          }
          else
          {
            bcopy((char *)&h0->yplot1, (char *)&H0.yplot2, nsize); /* Copy H1 array to H0's yplot2 */
          }
        }
        break;

      case 3:
        if(i+1 == req[0] || i+1 == req[1])
        {
          showC4();

          if(rb)
          {
            doBaseline();
          }

          line_MinMax(h0);

          min[i] = h0->datamin;
          max[i] = h0->datamax;
          if(min[i] < minT)
          {
            minT = min[i];
          }
          if(max[i] > maxT)
          {
            maxT = max[i];
          }

          nsize = h0->head.noint * sizeof(float);

          if(first)
          {
            _copyImage(h0, &H0);

            _h2();  /* Point to H2 */

            first = 0;
          }
          else
          {
            bcopy((char *)&h0->yplot1, (char *)&H0.yplot2, nsize); /* Copy H1 array to H0's yplot2 */
          }
        }
        break;
    }
  }

  env->ylower = minT;
  env->yupper = maxT;
  env->yscale = Y_FIXED_SCALE;
  H0.nsize    = nsize;
  H0.narray   = 2;                      /* Let xgraphic know there are 2 arrays */

  _h0();                                /* Point back to H0 */

  h0->ainfo.escan = H2.ainfo.escan; /* Tag the last scan */

  linePlot(h0);                 /* Plot it */

  /* restore environment */
  env->bkmatch    = save_bkmatch;
  env->check_feed = save_ckfeed;

  return(0);
}


int doBzero()
{
  int i, n;

  if(strlen(args->secondArg))
  {
    n = atoi(args->secondArg);

    if(n > 0 && n < h0->head.noint)
    {
      st_report("Zeroing out %d channels at beginning of data", n);

      for(i=0;i<n;i++)
      {
        h0->yplot1[i] = 0.0;
      }
    }
    else
    if(n < 0 && abs(n) < h0->head.noint)
    {
      n = abs(n);

      st_report("Zeroing out %d channels at end of data", n);

      for(i=h0->head.noint-n;i<h0->head.noint;i++)
      {
        h0->yplot1[i] = 0.0;
      }
    }

  }

  return(0);
}


/* Invert the spectra in h0;
   Copied from unipops/source/conline/line/au2.F
 */
int doSpecInvert()
{
  int i, noint;
  float *s, *d, hold[MAX_CHANS];

  noint = h0->head.noint;

  st_report("Inverting %d Spectal Data", noint);

  /* Save the data */
  bcopy((char *)&h0->yplot1[0], (char *)&hold[0], sizeof(float)*noint);
  
  /* Invert the data */
  s = &hold[noint-1];
  d = &h0->yplot1[0];
  for(i=0;i<noint;i++,d++,s--)
  {
    *d = *s;
  }

  /* fix header variables */

  h0->head.deltax  *= -1.0;
  h0->head.refpt    =  h0->head.noint - h0->head.refpt + h0->head.spn;
  h0->head.freqres *= -1.0;

  return(0);
}

/* OSB -- Other Side Band

   Copied from unipops/source/conline/line/au8.F

 */

int doOSB()
{
  struct HEADER *h;
  double fshift;

  h = &h0->head;

  /* flip sign of c12fr, and c12dx */
  h->freqres *= -1.0;
  h->deltax  *= -1.0;

  /* recalculate the frequencies, first the shift */
  if(h->sideband == 0.0 || h->sideband == 2.0)  /* USB -> LSB */
  {
    fshift = -2.0 * h->firstif;
  }
  else /* LSB -> USB */
  {
    fshift =  2.0 * h->firstif;
  }

  if(h->obsfreq != 0.0)
  {
    h->restfreq = h->restfreq * (1.0 + fshift / h->obsfreq);
  }
  else
  {
    h->restfreq = h->restfreq + fshift;
  }

  /* then the center freq, which is just a straight shift */
  h->obsfreq = h->obsfreq + fshift;

  /* reset the sideband value */
  if(h->sideband == 0.0)
  {
    h->sideband = 1.0;
  }
  else
  if(h->sideband == 1.0)
  {
    h->sideband = 0.0;
  }
  else
  if(h->sideband == 2.0)
  {
    h->sideband = 3.0;
  }
  else
  {
    h->sideband = 2.0;
  }

  /* and set/unset a flag in class 10 */
  /* I don't care to do this as; well don't see any need to */
  /* OK, talked me into it; set a flag in openpar[0] */

  if(h->openpar[0] == 54321.12345)
  {
    h->openpar[0] = 0.0;
  }
  else
  {
    h->openpar[0] = 54321.12345;
  }

  return(0);
}

#define SIN_PAT 1
#define COS_PAT 2
#define TAN_PAT 3
#define SQR_PAT 4

int fillPattern()
{
  int i, n, pat = SIN_PAT;
  float *f;
  double freq, fr, f2, f1, nchans, delta;

  if(args->ntokens < 2)
  {
    st_fail("Missing args for fillpattern: Usage: fillpattern freq");

    return(1); 
  }

  if(isSet(PF_S_FLAG))
  {
    pat = SIN_PAT;
  }

  if(isSet(PF_T_FLAG))
  {
    pat = TAN_PAT;
  }

  if(isSet(PF_C_FLAG))
  {
    pat = COS_PAT;
  }

  if(isSet(PF_Q_FLAG))
  {
    pat = SQR_PAT;
  }

  freq   = atof(args->secondArg);
  fr     = fabs(h0->head.freqres);
  n      = h0->head.noint;
  nchans = freq / fr;
  f1     = nchans / 360.0;
  delta  = 1.0 / f1;

  st_report("freq = %.2f, fr = %.2f, nchans = %.2f, f1 = %.2f, delta = %.2f", 
                freq, fr, nchans, f1, delta);

  f = &h0->yplot1[0];
  f2 = 0.0;
  for(i=0;i<n;i++,f++)
  {
    switch(pat)
    {
      case SIN_PAT: *f = (float)sin(f2 * DEG_TO_RAD); break;
      case COS_PAT: *f = (float)cos(f2 * DEG_TO_RAD); break;
      case SQR_PAT: 
        *f = (float)sin(f2 * DEG_TO_RAD); 
        if(*f < -0.0)
        {
          *f = -1.0;
        }
        if(*f > 0.0)
        {
          *f = 1.0;
        }
        break;

      case TAN_PAT: 
        *f = (float)tan(f2 * DEG_TO_RAD); 
        if(*f < -10.0)
        {
          *f = -10.0;
        }
        if(*f > 10.0)
        {
          *f = 10.0;
        }
        break;
    }

    f2 += delta;

    if(f2 >= 360.0)
    {
      f2 = 0.0;
    }
  }

  return(0);
}


/* Set y-azis to 10 * log10(y/peak) */
int lineConvert2LogPlot()
{
  int i, n;
  float *p, max=-1e12;

  n = h0->head.noint;

  st_report("Converting %d floats to log-plot", n);

  
  p = &h0->yplot1[0];
  for(i=0;i<n;i++,p++)
  {
    if(*p > max)
    {
      max = *p;
    }
  }

  p = &h0->yplot1[0];
  for(i=0;i<n;i++,p++)
  {
    *p = 10.0 * log10(*p / max);
  }

  h0->pinfo.yplotmode = LOG_PLOT;

  return(0);
}



int line_save_data(h)
struct LP_IMAGE_INFO *h;
{
  int    i, dlen;
  char fname[256];
  FILE *fp;

  if(args->ntokens < 2)
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  dlen = (int)h->head.noint;

  sprintf(fname, "%s", args->secondArg);

  if((fp = fopen(fname, "w")) == NULL)
  {
    st_fail("Unable to open %s for writting", fname);

    return(0);
  }

  st_report("Saving %d records to: %s", dlen, fname);

  fprintf(fp, "# Scan Number    %8.2f\n", h->head.scan);
  fprintf(fp, "# Bscan Number   %8.2f\n", h->ainfo.bscan);
  fprintf(fp, "# Escan Number   %8.2f\n", h->ainfo.escan);
  fprintf(fp, "# Scan Count     %8d\n",   h->ainfo.scanct);
  fprintf(fp, "# \n");

  for(i=0;i<dlen;i++)
  {
    fprintf(fp, "%5d %12.6f\n", i, h->yplot1[i]);
  }

  fclose(fp);

  return(0);
}


int deleteReport()
{
  system("rm -f report.log");

  return(0);
}


int writeReport(buf)
char *buf;
{
  FILE *fp;
  char *fname = "report.log";

  if((fp = fopen(fname, "a")) == NULL)
  {
    st_fail("Unable to open %s for writting", fname);

    return(0);
  }

  fprintf(fp, "%s\n", buf);

  fclose(fp);

  return(0);
}


/* Print various scan varibles */
int printReport()
{
  int tmp;
  char outBuf[256], tbuf[512], src[256];

  if(args->ntokens > 1)
  {
    if(!strncmp(args->secondArg, "delete", strlen(args->secondArg)))
    {
      deleteReport();

      return(0);
    }
  }

  strcpy(outBuf, "REPORT: ");

  if(isSet(PF_O_FLAG))
  {
    sprintf(src, "%16.16s ", h0->head.object);
    strcmprs(src);

    sprintf(tbuf, "Source: %s ", src);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_S_FLAG))
  {
    sprintf(tbuf, "Scan: %.2f ", h0->head.scan);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_D_FLAG))
  {
    sprintf(tbuf, "Date: %.4f ", h0->head.utdate);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_U_FLAG))
  {
    sprintf(tbuf, "UT: %.2f ", h0->head.ut);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_L_FLAG))
  {
    sprintf(tbuf, "LST: %.2f ", h0->head.lst);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_A_FLAG))
  {
    sprintf(tbuf, "Az: %.1f ", h0->head.az);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_E_FLAG))  /* El flag */
  {
    sprintf(tbuf, "El: %.1f ", h0->head.el);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_I_FLAG))
  {
    sprintf(tbuf, "Time: %.0f ", h0->head.inttime);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_P_FLAG))
  {
    sprintf(tbuf, "Pk: %.2f ", env->datamax);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_R_FLAG))
  {
    strcpy(tbuf, "rms"); /* Go compute rms */

    tmp = stquiet;
    stquiet = 1;

    parseCmd(tbuf);

    stquiet = tmp;

    sprintf(tbuf, "rms: %.2f mK ", env->rms * 1000.0);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_T_FLAG))
  {
    sprintf(tbuf, "Tsys: %.1f ", h0->head.stsys);
    strcat(outBuf, tbuf);
  }

  if(isSet(PF_W_FLAG))
  {
    writeReport(outBuf);
  }
  else
  {
    strcat(outBuf, "\n");
    st_print(outBuf);
  }

  return(0);
}


int isSpecFivePoint(h)
struct HEADER *h;
{
  if(h->openpar[8] == 9998.50 || h->openpar[8] == 9999.50)
  {
    return(1);
  }

  return(0);
}


int line_setScale()
{
  int i;

  if(args->ntokens > 1)
  {
    if(strlen(args->secondArg))
    {
      env->cal_scale = atof(args->secondArg);

      st_report("Scaling Data by %.3f", env->cal_scale);

      h0->head.tcal  /= env->cal_scale;
      h0->head.stsys /= env->cal_scale;

      for(i=0;i<h0->head.noint;i++)
      {
        h0->yplot1[i] /= env->cal_scale;
      }

      xx();
    }
    else
    {
      printH_FlagHelp(args->firstArg);
    }
  }

  return(0);
}


int doZero()
{
  int bd, ed, i, n;
  double avg;
  struct LP_IMAGE_INFO *image = h0;

  _doCCur("Click Left Mouse Button on Left Edge of Desired Region");

  bd = (int)round(env->ccur);

  if(bd < 1)
  {
    bd = 1;
  }

  _doCCur("Click Left Mouse Button on Right Edge of Desired Region");

  ed = (int)round(env->ccur);

  n = image->head.noint;

  if(ed >= n)
  {
    ed = n -1;
  }

  if(isSet(PF_Z_FLAG)) /* Force them to zero ?? */
  {
    avg = 0.0;
  }
  else /* Default is to average the end points */
  {
    avg = (image->yplot1[bd] + image->yplot1[ed]) / 2.0;
  }

  st_report("Setting Channels %d - %d to %.4f", bd, ed, avg);

  for(i=bd;i<ed;i++)
  {
    image->yplot1[i] = avg;
  }

  return(0);
}


int doTrec()
{
  int i, n, ret;
  char tbuf[256];
  double bscan, havg, savg, cavg;
  double yfac, syfac, avyfac, trec;
  double xval, yval;

  if(args->ntokens > 1)
  {
    if(strlen(args->secondArg))
    {
      bscan = atof(args->secondArg);

      st_report("T-Receiver: %.2f", bscan);

      /* Get hot load */
      _h1();
      bscan += 0.20;
      subArgs(2, "%.2f", bscan);
      ret = ggetScan("Vane");

      if(ret)
      {
        st_fail("doTrec(): Hot-Load Scan %.2f not found", bscan);

        return(1);
      }

      showAverage("Vane");
      havg = env->bp_avg;

      /* Get sky load */
      _h2();
      bscan += 0.20;
      subArgs(2, "%.2f", bscan);
      ret = ggetScan("Sky");


      if(ret)
      {
        st_fail("doTrec(): Sky Scan %.2f not found", bscan);

        return(1);
      }

      showAverage("Sky");
      savg = env->bp_avg;

      syfac = havg / savg;

      st_report("Sky Yfac = %.2f", syfac);

      /* Get cold load */
      _h3();
      bscan += 0.20;
      subArgs(2, "%.2f", bscan);
      ret = ggetScan("Cold");

      if(ret)
      {
        st_fail("doTrec(): Cold-Load Scan %.2f not found", bscan);

        return(1);
      }

      showAverage("Cold");
      cavg = env->bp_avg;

      avyfac = havg / cavg;
      trec = (273.0+H1.head.tchop-avyfac*H1.head.tcold) / (avyfac - 1.0);

      st_report("Yfac = %.2f, Trec = %.2f", avyfac, trec);

      _h0();

      hzero(&H0);

        /* Copy the header */
      bcopy((char *)&H1.head, (char *)&H0.head, sizeof(H0.head));

      n = H1.head.noint;

      for(i=env->bdrop;i<(n-env->edrop);i++)
      {
        if(H3.yplot1[i] > 0.0)
        {
          yfac = H1.yplot1[i] / H3.yplot1[i];
        }
        else
        {
          yfac = H1.yplot1[i] / cavg;
        }

        H0.yplot1[i] = (273.0+H1.head.tchop-yfac*H1.head.tcold) / (yfac - 1.0);
      }

      H0.head.scan = atof(args->secondArg); /* Set the scan number back */

      showScan();

      getResponse(WAIT_FOR_PLOT); /* Wait for xgraphics to plot and respond */

      sprintf(tbuf, "hline /Cmg /c 0 %f", trec);
      parseCmd(tbuf);

      strcpy(h0->title, "T-receiver");

      n = (int)(H1.head.noint * 0.10);

      xval = H0.x1data[n];
      yval = (env->yplotmax - env->yplotmin) * 0.90 + env->yplotmin;

      sprintf(tbuf, "label set 0 %f %f 5 0 T_{rec}: %.1f, Y_{fac}: %.2f /3", xval, yval, trec, avyfac);
      parseCmd(tbuf);
    }
    else
    {
      printH_FlagHelp(args->firstArg);
    }
  }

  return(0);
}


int doYFactor()
{
  int    ret;
  char   tbuf[256];
  double bscan, havg, savg;
  double yfac;

  if(args->ntokens > 1)
  {
    if(strlen(args->secondArg))
    {
      bscan = atof(args->secondArg);

      st_report("doYFactor: %.2f", bscan);

      /* Get hot load */
      _h1();
      bscan += 0.20;
      subArgs(2, "%.2f", bscan);
      ret = ggetScan("Vane");

      if(ret)
      {
        st_fail("doYFactor(): Hot-Load Scan %.2f not found", bscan);

        return(1);
      }

      showAverage("Vane");
      havg = env->bp_avg;

      /* Get sky load */
      _h2();
      bscan += 0.20;
      subArgs(2, "%.2f", bscan);
      ret = ggetScan("Sky");


      if(ret)
      {
        st_fail("doYFactor(): Sky Scan %.2f not found", bscan);

        return(1);
      }

      showAverage("Sky");
      savg = env->bp_avg;

      yfac = havg / savg;

      st_report("Sky Yfac = %.2f", yfac);

      sprintf(tbuf, "image /y h1 h2");
      parseCmd(tbuf);

      _h0();

      strcpy(h0->title, "Y-Factors");

      showScan();

      getResponse(WAIT_FOR_PLOT); /* Wait for xgraphics to plot and respond */

      sprintf(tbuf, "hline /Cmg /c 0 %f", yfac);
      parseCmd(tbuf);

    }
    else
    {
      printH_FlagHelp(args->firstArg);
    }
  }

  return(0);
}


/* Duplicate the Unipops 'nsac' command.
 *
 * Recall any number of nsaves; shift them if necessary; sum them;
 */
int doNsac()
{
  int i, n, ret, nsave, nscans=0;
  int argsPassed = 0;
  char tbuf[256];
  double freq, newfreq, fscan = 1e9, lscan = 0.0;
  char source[80], newsource[80];

  if(!strlen(env->savefile))
  {
    setColor(RED);
    st_report("Save file not set");
    setColor(BLACK);

    return(1);
  }

  saveEnv();

  argsPassed = args->ntokens - 1;
  
  if(!argsPassed)
  {
    n = getNumericResponse("Enter Number of 'nsaves' to average < 9: ");
  }
  else
  {
    n = atoi(args->tok[1]);
  }
  
  if(n == 0)
  {
    setColor(RED);
    st_report("nsac(): Cancelled");
    setColor(BLACK);
    restoreEnv();

    return(0);
  }

  if(n > 9)
  {
    setColor(RED);
    st_report("nsac(): Number of nsaves, %d, too large; 9 Max", n);
    setColor(BLACK);
    restoreEnv();

    return(1);
  }

//  st_report("ntokens = %d, n+2 = %d\n", args->ntokens, n+2);

  if(argsPassed && args->ntokens != n+2)
  {
    setColor(RED);
    st_report("nsac(): Wrong Number of arguments for command");
    setColor(BLACK);
    restoreEnv();

    return(1);
  }

  if(isSet(PF_C_FLAG))
  {
    clearAllImages();
  }

  /* Fetch all the nsaves first */
  for(i=0;i<n;i++)
  {
    if(!argsPassed)
    {
      nsave = getNumericResponse("Enter nsave: %d of %d: ", i+1, n);
    }
    else
    {
      nsave = atoi(args->tok[i+2]);
    }

    if(nsave == 0)
    {
      setColor(RED);
      st_report("nsac(): Cancelled");
      setColor(BLACK);
      restoreEnv();

      return(0);
    }

    env->nsave = nsave;

    sprintf(tbuf, "H%d", i+1);
    subArgs(2, tbuf);
    setHeader();
                                /* The first one is the refereance */
    ret = recallNSave();

    if(ret)
    {
      setColor(RED);
      st_report("Error retreiving nasave %d", nsave);
      setColor(BLACK);
      restoreEnv();

      return(1);
    }

    sprintf(tbuf, "NSAVE_%d", nsave);
    subArgs(2, tbuf);
    strcpy(h0->title, tbuf);

    if(i == 0)
    {
      strncpy(source, h0->head.object, 16);
      source[16] = '\0';
      strcmprs(source);

      freq = h0->head.restfreq;
    }
    else
    {
      strncpy(newsource, h0->head.object, 16);
      newsource[16] = '\0';
      strcmprs(newsource);

      if(strncmp(source, newsource, 16))
      {
        setColor(RED);
        st_report("Sources are not the same: %s, %s", source, newsource);
        setColor(BLACK);
        restoreEnv();

        return(1);
      }

      newfreq = h0->head.restfreq;
      if(fabs(newfreq-freq) > 50.0) /* 50 MHz should cover all shifts */
      {
        setColor(RED);
        st_report("RestFreq's are not the same: %f, %f", newfreq, freq);
        setColor(BLACK);
        restoreEnv();

        return(1);
      }

      subArgs(2, "H1");
      ret = compAutoShift();

      if(ret)
      {
        setColor(RED);
        st_report("nsac(): Shift Failed");
        setColor(BLACK);
        restoreEnv();

        return(1);
      }
    }

    h0->ainfo.bscan  = h0->head.fscan;
    h0->ainfo.escan  = h0->head.lscan;
    h0->ainfo.scanct = h0->head.nostac;

    nscans += h0->head.nostac;

    if(h0->head.fscan < fscan)
    {
      fscan = h0->head.fscan;
    }

    if(h0->head.lscan > lscan)
    {
      lscan = h0->head.lscan;
    }
  }

  /* Now sum them all up */

  for(i=0;i<n;i++)
  {
    if(i == 0) /* Copy the first one into holder 0 */
    {
      subArgs(2, "H1");
      subArgs(3, "H0");

      copyImage();
    }
    else
    {
      sprintf(tbuf, "H%d", i+1);
      subArgs(2, "H0");
      subArgs(3, tbuf);

      doSum();
    }
  }

  subArgs(2, "H0");
  setHeader();
  strcpy(h0->title, "Sum");

  h0->ainfo.bscan  = h0->head.fscan  = fscan;
  h0->ainfo.escan  = h0->head.lscan  = lscan;
  h0->ainfo.scanct = h0->head.nostac = nscans;

  showScan();

  imageSummary();

  restoreEnv();

  return(0);
}


int doDedit()
{
  int i, n, ret, which = 0;
  char fname[256], sysBuf[512], *string, line[256];
  float *f, ff;
  FILE *fp;

  if(isSet(PF_0_FLAG))
  {
    which = 0;
  }
  else
  if(isSet(PF_1_FLAG))
  {
    which = 1;
  }
  else
  if(isSet(PF_2_FLAG))
  {
    which = 2;
  }
  else
  if(isSet(PF_3_FLAG))
  {
    which = 3;
  }

  if(isSet(PF_S_FLAG))  /* Save */
  {
    sprintf(fname, "/tmp/%s/dedit.txt", session);

    st_report("Saving Data[%d] array to: %s", which, fname);

    if((fp = fopen(fname, "w")) == NULL)
    {
      st_fail("Unable to open %s for writting", fname);

      return(0);
    }

    n =  h0->head.noint;

    switch(which)
    {
      case 0: f = &h0->yplot1[0]; break;
      case 1: f = &h0->yplot2[0]; break;
      case 2: f = &h0->yplot3[0]; break;
      case 3: f = &h0->yplot4[0]; break;
    }

    for(i=0;i<n;i++, f++)
    {
      fprintf(fp, "%5d %.10f\n", i, *f);
    }

    fclose(fp);

    return(0);
  }

  if(isSet(PF_E_FLAG))  /* Edit */
  {
    sprintf(fname, "/tmp/%s/dedit.txt", session);

    st_report("Editing %s", fname);

    sprintf(sysBuf, "/usr/bin/vim %s", fname);

    system(sysBuf);

    return(0);
  }

  if(isSet(PF_R_FLAG))  /* Read */
  {
    sprintf(fname, "/tmp/%s/dedit.txt", session);

    st_report("Reading into Data[%d], %s", which, fname);

    if((fp = fopen(fname, "r")) == NULL)
    {
      st_fail("Unable to open %s for reaading", fname);

      return(0);
    }

    n =  h0->head.noint;

    while((string = fgets(line, sizeof(line), fp)) != NULL)
    {
      ret = sscanf(line, "%d %f", &i, &ff);

      if(ret == 2)
      {
                        /* I can't assume the array is still in channel 
                           order, so use index */
        if(i>= 0 && i < n)
        {
          switch(which)
          {
            case 0: h0->yplot1[i] = ff; break;
            case 1: h0->yplot2[i] = ff; break;
            case 2: h0->yplot3[i] = ff; break;
            case 3: h0->yplot4[i] = ff; break;
            case 4: h0->bplot [i] = ff; break;
            case 5: h0->gplot [i] = ff; break;
          }
        }
        else
        {
          st_print("doDedit(): %d, greater then noint: %d", i, n);
        }
      }
      else
      {
        st_print("doDedit(): sacanf() returned %d, and not 2", ret);
      }
    }

    xx();

    fclose(fp);

    return(0);
  }

  if(isSet(PF_U_FLAG))  /* User action flag */
  {
    if(args->ntokens > 1)
    {
      sprintf(fname, "/tmp/%s/dedit.txt", session);

      sprintf(sysBuf, "%s %s", args->secondArg, fname);

      st_report("Executing: %s", fname);

      system(sysBuf);

      return(0);
    }

    return(0);
  }

  printH_FlagHelp(args->firstArg);

  return(0);
}