gridmap.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"

struct BEAM_MAP beam_map;

char   middleRowFname[256];

double sqrarg;

float global_gdata[MAX_CONT_SAMPLES];
int   use_global_gdata;

int    doingDive = 0;

#define SQR(a) (sqrarg=(a),sqrarg*sqrarg)

#define MAX_NUM_BEAM_SEP 2048

struct BEAM_MAPS_SEP {
        int used;
        int pad;
        double scan;
        double seperation;
        double bw1;
        double bw2;
        double peak1;
        double peak2;
        double ratio;
};

struct BEAM_MAPS_SEP beam_map_sep[MAX_NUM_BEAM_SEP];

struct PARSE_FLAGS gmap_pflags[MAX_PARSE_FLAGS];

struct SCAN_HEADERS {
        int   ver;
        float scan;

        struct HEADER head;

        struct SCAN_HEADERS *next;
};

struct SCAN_HEADERS *scan_headers=NULL;


struct HEADER *do_addHeader(ver, scan, head)
int ver;
float scan;
struct HEADER *head;
{
  struct HEADER *h;

  if((h = checkScanHeader(ver, scan)))
  {
    return(h);
  }
  else
  {
    if((h = addScanHeader(ver, scan, head)))
    {
      return(h);
    }
    else
    {
      log_msg("do_addHeader(): Error adding scan %.2f, in ver %d", scan, ver);
    }
  }

  return((struct HEADER *)NULL);
}


struct HEADER *checkScanHeader(v, s)
int v;
float s;
{
  struct SCAN_HEADERS *m;
  float ss;

  if(scan_headers == NULL)
  {
    return((struct HEADER *)NULL);
  }

  /* users can pass a scan number like 1234.0, which will ignore all
     scans with that prefix */
  ss = (float)(int)s;

  for(m=scan_headers;m->next; m=m->next)
  {
    if(m)
    {
      if((m->scan == s || m->scan == ss) && m->ver == v)
      {
        return(&m->head);
      }
    }
  }

  /* Check the last one, if valid */
  if(m)
  {
    if((m->scan == s || m->scan == ss) && m->ver == v)
    {
      return(&m->head);
    }
  }

  return((struct HEADER *)NULL);
}


struct HEADER *addScanHeader(v, s, head)
int v;
float s;
struct HEADER *head;
{
  struct SCAN_HEADERS *m;

  /* Create a fresh list */
  if(scan_headers == NULL)
  {
    if((scan_headers=(struct SCAN_HEADERS *)malloc((unsigned)(sizeof(struct SCAN_HEADERS))))==NULL)
    {
      log_msg("Error in addIgnoreScans malloc(%d)", sizeof(struct SCAN_HEADERS));

      return((struct HEADER *)NULL);
    }

    scan_headers->ver  = v;
    scan_headers->scan = s;

    bcopy((char *)head, (char *)&scan_headers->head, sizeof(struct HEADER));

    scan_headers->next = NULL;

    return(&scan_headers->head);
  }

  /* Find next available slot */
  for(m=scan_headers;m->next; m=m->next)
  {
    ;
  }

  /* Append new scan */
  if((m->next=(struct SCAN_HEADERS *)malloc((unsigned)(sizeof(struct SCAN_HEADERS))))==NULL)
  {
    log_msg("Error in addIgnoreScans malloc(%d)", sizeof(struct SCAN_HEADERS));

    return((struct HEADER *)NULL);
  }

  m->next->ver  = v;
  m->next->scan = s;

  bcopy((char *)head, (char *)&m->next->head, sizeof(struct HEADER));

  m->next->next = NULL;

  return(&m->next->head);
}



int remScanHeaders()
{
  struct SCAN_HEADERS *p;
  struct SCAN_HEADERS *q;

  for(p=scan_headers; p; p=q)
  {
    q = p->next;

    if(p)
    {
      free((char *)p);
    }
  }

  scan_headers = NULL;

  return(0);
}



int clearScreen()
{
  fprintf(stdout, "%c[0m" , 0x1b);
  fflush(stdout);
  fprintf(stdout, "%c[2J", 0x1b);
  fflush(stdout);

  return(0);
}


int gotoxy(x, y)
int x, y;
{
  fprintf(stdout, "%c[%d;%dH" , 0x1b, y, x);
  fflush(stdout);

  return(0);
}


int gridmap_goWriteCenterRow(xdata, ydata, aa1, aa2, n, scanno)
double *xdata, *ydata, *aa1, *aa2, scanno;
int n;
{
  int i;
  FILE *fp;
  double g1, g2, y, dyda[4];

  sprintf(middleRowFname, "/tmp/%s/middle_row.%04.2f.dat", session, scanno);

  if((fp = fopen(middleRowFname, "w")) == NULL)
  {
    fprintf(stderr, "Unable to open output file %s\n", middleRowFname);

    exit(1);
  }

  for(i=0;i<n;i++)
  {
    fgauss(xdata[i], aa1, &y, dyda);
    g1 = y;

    fgauss(xdata[i], aa2, &y, dyda);
    g2 = y;

    fprintf(fp, "%f %f %f %f\n", xdata[i], ydata[i], g1, g2);

    global_gdata[i] = g1 + g2;
  }

  fclose(fp);


  return(0);
}


int gridmap_showMiddleRow()
{
  char tbuf[512];

  if(strlen(middleRowFname) > 0)
  {
    sprintf(tbuf, "show_mid.csh %s %s &", middleRowFname, session);
    system(tbuf);
  }
  else
  {
    st_print("No Middle Row to Plot\n");
    st_print("Have you done a 'gridmap' and a 'mapgauss'?\n");
  }

  return(0);
}


int gridmap_write_datafile()
{
  int x, y;
  FILE *fp;
  char outfile[256];
  struct BEAM_MAP_ELEM *p;

  if(beam_map.map_h == 0 || beam_map.map_w == 0)
  {
    st_fail("No map loaded");
    return(0);
  }

  sprintf(outfile, "/tmp/%s/beam_map.dat", session);

  if((fp = fopen(outfile, "w")) == NULL)
  {
    st_print("Unable to open output file %s\n", outfile);

    return(1);
  }

  for(y=beam_map.plot_yb;y<beam_map.plot_ye;y++)
  {
    p = &beam_map.bmes[y][beam_map.plot_xb];

    if(!beam_map.bmes[y][0].n) /* Row not populated yet? */
    {
      continue;
    }

    for(x=beam_map.plot_xb;x<beam_map.plot_xe;x++,p++)
    {
      fprintf(fp, "%5.2f %5.2f %10.6f\n", p->eloff, p->azoff, p->pval);
    }

    fprintf(fp, "\n");
  }

  fclose(fp);

  return(0);
}


int gridmap_write_csvfile()
{
  int x, y;
  FILE *fp;
  char outfile[256], source[80];
  double pvalue;
  struct BEAM_MAP_ELEM *p;
  struct HEADER *h;

  if(beam_map.map_h == 0 || beam_map.map_w == 0)
  {
    st_fail("No map loaded");
    return(0);
  }

  p = &beam_map.bmes[beam_map.map_h /2][0];
  h = p->headers[0];                                    /* Grab the center row header */

  if(isSet(PF_A_FLAG))                                  /* Look for Auto-Name flag */
  {
    strncpy(source, h->object, 16);
    source[16] = '\0';
    strcmprs(source);
    sprintf(outfile, "%s-BeamMap-%.1f-IF%1d.csv", source, h->restfreq / 1e3, (int)round(env->if_feed * 100.0));
    strcmprs(outfile);
  }
  else
  if(args->ntokens > 1)                                 /* Look for a passed filename */
  {
    strcpy(outfile, args->secondArg);

    if(!strstr(outfile, ".csv"))                        /* check for file-extension */
    {
      strcat(outfile, ".csv");
    }
  }
  else                                  /* If not, Use a default */
  {
    sprintf(outfile, "beam_map.csv");
  }

  if((fp = fopen(outfile, "w")) == NULL)
  {
    st_print("Unable to open output file %s\n", outfile);

    return(1);
  }


  fprintf(fp, "# MAP Source  (Name): %16.16s\n", h->object);
  fprintf(fp, "# MAP X-Axis  (Cols): %.0f\n",    h->noxpts);
  fprintf(fp, "# MAP Y-Axis  (Rows): %.0f\n",    h->noypts);
  fprintf(fp, "# MAP Delta-X (Secs): %.0f\n",    h->deltaxr);
  fprintf(fp, "# MAP Delta-Y (Secs): %.0f\n",    h->deltayr);
  fprintf(fp, "# MAP Freq    (GHz ): %.6f\n",    h->restfreq / 1e3);
  fprintf(fp, "# MAP Bk-End        : %8.8s\n",   h->backend);

  for(y=0;y<beam_map.map_h;y++)
  {
    fprintf(fp, "%.1f, ", (float)y - (float)beam_map.map_h / 2.0 + 0.5);

    for(x=0;x<beam_map.map_w;x++)
    {
      pvalue = beam_map.bmes[y][x].pval;

      fprintf(fp, "%f", pvalue);

      if(x < (beam_map.map_w-1))
      {
        fprintf(fp, ", ");
      }
    }

    fprintf(fp, "\n");
  }

  if(verboseLevel & VERBOSE_GRIDMAP)
  {
    st_report("Wrote %d x %d map to %s", beam_map.map_w, beam_map.map_h, outfile);
  }

  fclose(fp);

  return(0);
}


int gridMapStats()
{
  int    i, x, y, j, n, mid;
  int    plotval = 0, scanct  = 0, scans   = 0, tsys    = 0, 
         rval    = 0, itime   = 0, offsets = 0, writeit = 0;
  int    xstart, xend, ystart, yend;
  char   fname[256], blanks[80];
  struct BEAM_MAP_ELEM *p;
  struct BEAM_MAP *bm;
  FILE   *fp, *saveFp;

  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  if(beam_map.map_h == 0 || beam_map.map_w == 0)
  {
    st_fail("No map loaded");
    return(0);
  }

  bm = &beam_map;

  if(isSet(PF_C_FLAG))                                  /* Display Scan Counts; default */
  {
    scanct = 1;
  }
  else
  if(isSet(PF_S_FLAG))                                  /* Display Scan Counts; default */
  {
    scans = 1;
  }
  else
  if(isSet(PF_R_FLAG))                                  /* Display Raw Value */
  {
    rval = 1;
  }
  else
  if(isSet(PF_P_FLAG))                                  /* Display Plot Value */
  {
    plotval = 1;
  }
  else
  if(isSet(PF_T_FLAG))                                  /* Display Tsys Value */
  {
    tsys = 1;
  }
  else
  if(isSet(PF_I_FLAG))                                  /* Display Intigration Time */
  {
    itime = 1;
  }
  else
  if(isSet(PF_O_FLAG))                                  /* Display Map Offsets */
  {
    offsets = 1;
  }

  /* The WRITE flag indicates if the user wants
     to write out the results too;

     Hijack the library printing function to log the
     output for the duration of this function */
  if(isSet(PF_W_FLAG))
  {
    if(args->ntokens > 1 && strlen(args->secondArg) > 1)
    {
      strcpy(fname, args->secondArg);

      st_report("gridMapStats(): Writing outputfile: %s", fname);

      if((fp = fopen(fname, "w") ) <= (FILE *)0)
      {
        st_fail("Unable to open %s", fname);
        writeit = 0;
      }
      else
      {
        saveFp = errPrint_fp;
        errPrint_fp = fp;
        writeit = 1;
      }
    }
    else
    {
      st_fail("Missing filename arg for '/w' flag");

      printH_FlagHelp(args->firstArg);

      return(1);
    }
  }

  strcpy(blanks, "                                     ");

  n = 8 - strlen(bm->source);
  if(n > 0)
  {
    blanks[n] = '\0';
    st_print("# MAP Source  (Name):  %s%.16s\n", blanks, bm->source);
  }
  else
  {
    st_print("# MAP Source  (Name):  %.16s\n", bm->source);
  }

  st_print("# MAP X-Axis  (Cols):  %8d\n",     bm->map_w);
  st_print("# MAP Y-Axis  (Rows):  %8d\n",     bm->map_h);
  st_print("# MAP Delta-X (Secs):  %8.0f\n",   bm->map_step_x);
  st_print("# MAP Delta-Y (Secs):  %8.0f\n",   bm->map_step_y);
  st_print("# MAP Size-X  (Secs):  %8.0f\n",   bm->map_size_x);
  st_print("# MAP Size-Y  (Secs):  %8.0f\n",   bm->map_size_y);

  /* Slice through the map horizonatally in the middle and look for
     values to determine the beginning and ending rows; only works if
     map is somewhat populated */
  xstart = -1;
  xend   =  0;

  mid = bm->map_h / 2;
  for(i=0;i<bm->plot_xe;i++)
  {
    if((bm->bmes[mid][i].rval != 0.0) && (xstart == -1))
    {
      xstart = i;
    }

    if(bm->bmes[mid][i].rval != 0.0)
    {
      xend = i;
    }
  }

  st_print("# MAP Start-X:         %8d\n",     xstart +1);
  st_print("# MAP End-X:           %8d\n",     xend   +1);

  /* Now do the same vertically through center */
  ystart = -1;
  yend   =  0;

  mid = bm->map_w / 2;
  for(i=0;i<bm->plot_ye;i++)
  {
    if((bm->bmes[i][mid].rval != 0.0) && (ystart == -1))
    {
      ystart = i;
    }

    if(bm->bmes[i][mid].rval != 0.0)
    {
      yend = i;
    }
  }

  st_print("# MAP Start-Y:         %8d\n",   ystart +1);
  st_print("# MAP End-Y:           %8d\n",   yend   +1);

  st_print("# MAP Freq    (GHz ):  %8.1f\n", bm->freq / 1e3);
  st_print("# MAP Num of Scans:    %8d\n",   bm->nscans);
  st_print("# MAP Blines Samples:  %8d\n",   bm->nbase);
  st_print("# MAP Blines Removed:  %8d\n",   bm->base_removed);
  st_print("# MAP First Scan:      %8.2f\n", bm->fscan);
  st_print("# MAP Last  Scan:      %8.2f\n", bm->lscan);
  st_print("# MAP Y-Label:         %8s\n",   bm->yaxis_label);

  st_print("#\n");

  /* No use going any farther */
  if((plotval + scanct + scans + tsys + rval + itime + offsets) == 0)
  {
    if(writeit)
    {
      fclose(fp);
      errPrint_fp = saveFp;
    }

    return(0);
  }

  if(scans)
  {
    for(y=0;y<bm->map_h;y++)
    {
      st_print("%2d ", y);

      p = &bm->bmes[y][0];

      st_print("%7.1f ", p->eloff);

      for(j=0;j<p->n;j++)
      {
        if(p && p->headers[j]->scan)
        {
          st_print("%4.2f ", p->headers[j]->scan); /* Stored as (int)(scan * 100.0) */
        }
        else
        {
          st_print("%4.2f ", 0.0);
        }
      }

      st_print("\n");
    }

    if(writeit)
    {
      fclose(fp);
      errPrint_fp = saveFp;
    }

    return(0);
  }


  /* don't place the 'Col' header in the output file */
  if(!writeit)
  {
    st_print("#\n# Col");

    for(x=0;x<bm->map_w;x++)
    {
      if(x == 0)
      {
        st_print("%4d,", x);
      }
      else
      {
        if(offsets)
        {
          st_print("%9d,", x);
        }
        else
        {
          st_print("%8d,", x);
        }
      }
    }

    st_print("\n#\n");
  }

  if(offsets)
  {
    for(y=0;y<bm->map_h;y++)
    {
        /* Don't write leading 'Row' values to output file */
      if(!writeit)
      {
        st_print("%2d ", y);
        st_print("%8.2f", bm->bmes[y][0].eloff);
      }

      for(x=1;x<bm->map_w;x++)
      {
        st_print("%8.2f", bm->bmes[y][x].azoff);

        if(x < (bm->map_w-1))
        {
          st_print(", ");
        }
      }

      st_print("\n");
    }

    if(writeit)
    {
      fclose(fp);
      errPrint_fp = saveFp;
    }

    return(0);
  }

  for(y=0;y<bm->map_h;y++)
  {
    if(!writeit)
    {
      st_print("%2d: ", y);
    }

    for(x=0;x<bm->map_w;x++)
    {

      if(scanct)
      {
        st_print("%2d", bm->bmes[y][x].n);
      }
      else
      if(rval)
      {
        st_print("%7.3f", bm->bmes[y][x].rval);
      }
      else
      if(plotval)
      {
        st_print("%7.3f", bm->bmes[y][x].pval);
      }
      else
      if(tsys)
      {
        st_print("%3.0f", bm->bmes[y][x].tsys);
      }
      else
      if(itime)
      {
        st_print("%2.0f", bm->bmes[y][x].itime);
      }

      if(x < (bm->map_w-1))
      {
        st_print(", ");
      }
    }

    st_print("\n");
  }

  if(writeit)
  {
    fclose(fp);
    errPrint_fp = saveFp;
  }

  return(0);
}



int gridmap_writeDemFile(off, how)
int off, how;
{
  int i, interp=0, max;
  FILE *fp;
  char source[80], fname[256], psfilename[256];
  double clevel, cval;
//  struct BEAM_MAP_ELEM *p;

  sprintf(fname, "/tmp/%s/bmap.dem", session);
  if((fp = fopen(fname, "w")) == NULL)
  {
    fprintf(stderr, "Unable to open output file \"%s\"\n", fname);

    return(1);
  }

  interp = isSet(PF_I_FLAG);            /* Interpolate level ?? */

//  p = &beam_map.bmes[off][0];

//  strncpy(source, p->headers[0]->object, 16);
//  source[16] = '\0';

  strcpy(source, env->source);

  fprintf(fp, "reset\n");
  fprintf(fp, "set title \"Beam Map of %s\"\n", source);
  fprintf(fp, "set xlabel \"Map El Offset\"\n");
  fprintf(fp, "set ylabel \"Map Az Offset\"\n");
  fprintf(fp, "set zlabel \"%s\"\n", beam_map.yaxis_label);

  fprintf(fp, "#\n");
  fprintf(fp, "set border 4095\n");
  fprintf(fp, "set hidden3d\n");
  fprintf(fp, "#\n");

  if(isSet(PF_A_FLAG))
  {
    max = beam_map.map_size_x > beam_map.map_size_y ? beam_map.map_size_x : beam_map.map_size_y;

    fprintf(fp, "set xrange [%.2f:%.2f]\n", max / 2.0 - max, max / 2.0);
    fprintf(fp, "set yrange [%.2f:%.2f]\n", max / 2.0 - max, max / 2.0);
  }
  else
  {
    fprintf(fp, "set xrange [%.2f:%.2f]\n", beam_map.map_size_y / 2.0 - beam_map.map_size_y, beam_map.map_size_y / 2.0);
    fprintf(fp, "set yrange [%.2f:%.2f]\n", beam_map.map_size_x / 2.0 - beam_map.map_size_x, beam_map.map_size_x / 2.0);
  }

  if(env->gmap_zmin != 0.0 || env->gmap_zmax != 0.0)
  {
    if(isSet(PF_V_FLAG))
    {
      st_report("Using User supplied Z-Range of: [%.2f %.2f]", env->gmap_zmin, env->gmap_zmax);
    }

    fprintf(fp, "set zrange  [%.2f:%.2f]\n", env->gmap_zmin, env->gmap_zmax);
  }
  else
  {
    fprintf(fp, "set zrange  [%.2f:%.2f]\n", beam_map.map_min, beam_map.map_max);
  }

  fprintf(fp, "set cbrange [%.2f:%.2f]\n", beam_map.map_min, beam_map.map_max);
  fprintf(fp, "#\n");
  fprintf(fp, "set samples %d*%d\n", beam_map.map_w, beam_map.map_h);
  fprintf(fp, "set isosamples 60\n");

  fprintf(fp, "set ticslevel 0.1\n");

  fprintf(fp, "#\n");

  if(isSet(PF_B_FLAG))
  {
    if(isSet(PF_B_FLAG) > 1)            /* Reverse plot order */
    {
      fprintf(fp, "set pm3d at bs ");
    }
    else
    {
      fprintf(fp, "set pm3d at sb ");
    }

    if(interp)
    {
      fprintf(fp, "interpolate %d,%d\n", interp+1, interp+1);
    }
    else
    {
      fprintf(fp, "\n");
    }
  }
  else
  if(isSet(PF_W_FLAG))
  {
    fprintf(fp, "unset pm3d\n");
  }
  else
  {
    fprintf(fp, "set pm3d at s ");

    if(interp)
    {
      fprintf(fp, "interpolate %d,%d\n", interp+1, interp+1);
    }
    else
    {
      fprintf(fp, "\n");
    }
  }

  fprintf(fp, "#\n");

  fprintf(fp, "set palette color\n");

  if(isSet(PF_R_FLAG) == 1)
  {
    fprintf(fp, "set palette defined\n");
  }
  else
  if(isSet(PF_R_FLAG) == 2)
  {
    fprintf(fp, "set palette model RGB\n");
    fprintf(fp, "set palette defined ( 0 0 0 0, 1 0 0 1, 3 0 1 0, 4 1 0 0, 6 1 1 1 )\n");
  }
  else
  if(isSet(PF_R_FLAG) >= 3)
  {
    fprintf(fp, "set palette model HSV\n");
    fprintf(fp, "set palette defined ( 0 0 1 1, 1 1 1 1 )\n");
  }
  else
  if(isSet(PF_G_FLAG))
  {
    fprintf(fp, "set palette model RGB\n");
    fprintf(fp, "set palette defined ( 0 0 0 0, 1 1 1 1 )\n");
  }
  else
  {
    fprintf(fp, "set palette model RGB\n");
    fprintf(fp, "set palette defined ( 0 \"green\", 1 \"blue\", 2 \"red\", 3 \"orange\" )\n");
  }

  if(env->gmap_gamma != 0.0)
  {
    fprintf(fp, "set palette gamma %.2f\n", env->gmap_gamma);
  }
  else
  {
    fprintf(fp, "set palette gamma 1.0\n");
  }

  fprintf(fp, "#\n");

  fprintf(fp, "set view %.1f, %.1f, 1.0, 1.0\n", env->map_angle[0], env->map_angle[1]);

  if(how == GNUPLOT)
  {
    if(doingDive)
    {
      fprintf(fp, "set terminal x11 enhanced %d size %d, %d position 500, 0\n", 21, 1024, 1080 );
    }
    else
    {
      fprintf(fp, "set terminal x11 enhanced %d size %d, %d position 930, 0\n", 20, 1100, 1100 );
    }
  }
  else
  {
    fprintf(fp, "set terminal postscript enhanced size 17, 9.5625\n");
  }

  if(isSet(PF_L_FLAG))                  /* want contours */
  {
    if(isSet(PF_S_FLAG) == 1)
    {
      fprintf(fp, "set contour surface\n");
    }
    else
    if(isSet(PF_S_FLAG) == 2)
    {
      fprintf(fp, "set contour both\n");
    }
    else
    {
      fprintf(fp, "set contour base\n");
    }

    clevel = fabs(beam_map.map_min-beam_map.map_max) / 21.0;
  
    if(isSet(PF_V_FLAG))
    {
      st_print("Contour levels set to:\n");

      cval = beam_map.map_min;

      for(i=0;i<21;i++)
      {
        st_print("%6.1f ", cval);

        if(i == 6 || i == 13)
        {
          st_print("\n");
        }

        cval += clevel;
      }

      st_print("\n");
    }

    fprintf(fp, "set cntrparam levels incr %.2f,%.2f,%.2f\n", beam_map.map_min, clevel, beam_map.map_max);
  }

  if(how == GNUPLOT)
  {
    fprintf(fp, "splot \"/tmp/%s/beam_map.dat\" every 1::0::1024 with lines\n", session);
//    fprintf(fp, "pause -1 \"Hit Return When Finished\"\n");
  }
  else
  {
    sprintf(psfilename, "/tmp/%s/ccal_plot.ps", session);

    fprintf(fp, "set out \"%s\"\n", psfilename);
    fprintf(fp, "splot \"/tmp/%s/beam_map.dat\" every 1::0::1024 with lines\n", session);
    fprintf(fp, "\n");

    st_print("\nMade PS file: %s\n\n", psfilename);
  }

  fclose(fp);

  return(0);
}



int gridMapLimits()
{
  int xb, xe, yb, ye;

  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  if(isSet(PF_R_FLAG)) /* Reset */
  {
    beam_map.plot_xb = 0;
    beam_map.plot_xe = beam_map.map_w;
    beam_map.plot_yb = 0;
    beam_map.plot_ye = beam_map.map_h;

    return(0);
  }

  if(isSet(PF_S_FLAG)) /* Show */
  {
    st_report("Gridmap Plot Limits: Az: %d - %d; El: %d - %d",
        beam_map.plot_xb,
        beam_map.plot_xe,
        beam_map.plot_yb,
        beam_map.plot_ye);

    return(0);
  }

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

    return(1);
  }

  xb = atoi(args->secondArg);
  xe = atoi(args->thirdArg);
  yb = atoi(args->fourthArg);
  ye = atoi(args->fifthArg);

  if(xb < 0 || xb > beam_map.map_w)
  {
    st_fail("xb val: %d, out of range of %d %d", xb, 0, beam_map.map_w);

    return(0);
  }

  if(xe == xb)
  {
    st_fail("xe val: %d, can not be same as xb: %d", xe, xb);

    return(0);
  }

  if(xe > beam_map.map_w || xe < xb)
  {
    st_fail("xe val: %d, out of range of %d %d", xe, xb, beam_map.map_w);

    return(0);
  }

  if(yb < 0 || yb > beam_map.map_h)
  {
    st_fail("yb val: %d, out of range of %d %d", yb, 0, beam_map.map_h);

    return(0);
  }

  if(ye == yb)
  {
    st_fail("ye val: %d, can not be same as yb: %d", ye, yb);

    return(0);
  }

  if(ye > beam_map.map_h || ye <= yb)
  {
    st_fail("ye val: %d, out of range of %d %d", ye, yb, beam_map.map_h);

    return(0);
  }

  beam_map.plot_xb = xb;
  beam_map.plot_xe = xe;
  beam_map.plot_yb = yb;
  beam_map.plot_ye = ye;

  return(0);
}


int gridmap_fitGauss()
{
  int i, j, x, y, indx1, indx2, ret, hm, glen, request;
  double xdata[512], ydata[512];
  double arcsec, max1=0.0, max2=0.0;
  double xd[512], yd[512], sig[512];
  double aa1[4], aa2[4], peak1, peak2, beam1, beam2;
  char tbuf[256], ttbuf[256];
  struct BEAM_MAP_ELEM *p;

  if(beam_map.map_w == 0 || beam_map.map_h == 0)
  {
    st_print("\n No Map to fit a gaussian to\n");
    st_print("Have you done a 'gridmap' yet?\n\n");

    return(1);
  }
 
  glen = 5;

  bzero((char *)&xd[0], sizeof(xd));
  bzero((char *)&yd[0], sizeof(yd));

   y = beam_map.map_h / 2;
  hm = beam_map.map_w / 2;  /* Half map */

  if(args->ntokens > 1)
  {
    request = atoi(args->secondArg);

    request--; /* change it to zero based */

    if(request >= 0 && request < beam_map.map_h)
    {
      y = request;
    }
  }

  p = &beam_map.bmes[y][0];

  if(!p->itime)
  {
    st_print("\n No Map in row %d to fit a gaussian to\n", y);

    return(1);
  }

  arcsec = p->headers[0]->deltaxr; /* Grab the first header in this cell */

  tbuf[0] = '\0';
  for(j=0;j<p->n;j++)
  {
    sprintf(ttbuf, "%.2f ", p->headers[j]->scan);

    strcat(tbuf, ttbuf);
  }

  st_report( "Fitting Gaussian Using Scan(s) %s; Row %d; Arcsec = %.1f", tbuf, (int)p->headers[0]->ycell0-1, arcsec);

  for(x=0;x<beam_map.map_w;x++) /* Create the x,y-axis data; Find Min/Max too */
  {
    xdata[x] = (float)x * arcsec;

    if(isSet(PF_P_FLAG))
    {
      ydata[x] = beam_map.bmes[y][x].pval;
    }
    else
    {
      ydata[x] = beam_map.bmes[y][x].rval;
    }

    sig[x] = 1.0;

    if(ydata[x] > max1)
    {
      max1  = ydata[x];
      indx1 = x;
    }
    if(ydata[x] < max2)
    {
      max2  = ydata[x];
      indx2 = x;
    }
  }

  // Set up for gausian fitting

  if(max1 == 0.0 || max2 == 0.0)
  {
    st_fail("No data to do gausfit");
    return(1);
  }

  st_report("Guess, Plus  Peak: %7.2f @ Sample %d or %.2f\"",
                                max1, indx1, xdata[indx1]);
  st_report("Guess, Minus Peak: %7.2f @ Sample %d or %.2f\"",
                                max2, indx2, xdata[indx2]);
  st_report("Guess, Beam Sep:    %7.2f\"",
                                fabs(indx1-indx2) * arcsec);

  /* Compute ideal gaussian 1/2 power width in samples */
  glen = (int)(p->headers[0]->bfwhm / arcsec) + 1;

// Load first peak
  for(i=indx1-glen,j=0;i<indx1+glen;i++,j++)
  {
    xd[j] = xdata[i];
    yd[j] = ydata[i];
  }

  st_report("Fitting a %d Point Gauss to Plus Peak", j+1);
  aa1[1] = max1;
  aa1[2] = xdata[indx1];
  aa1[3] = (double)j * arcsec; /* Beam width */

  ret = gauss_fit(&xd[0], &yd[0], sig, j, aa1, isSet(PF_V_FLAG));

  if(!ret)
  {
    st_fail("1st Half Gauss Fit Failed");

    for(i=0;i<j;i++)
    {
      st_print("%2d %4.1f %4.1f\n", i, xd[i], yd[i]);
    }

    return(1);
  }

  st_report("Found Plus  Beam @ %7.2f\", Peak = %7.2f, Beam Width = %7.2f",
                aa1[2], aa1[1], fabs(aa1[3]));

  peak1 = aa1[2];
  beam1 = fabs((aa1[2] - xdata[hm]));

// Load second peak
  for(i=indx2-glen,j=0;i<indx2+glen;i++,j++)
  {
    xd[j] = xdata[i];
    yd[j] = ydata[i];
  }

  st_report("Fitting a %d Point Gauss to Minus Peak", j+1);
  aa2[1] = max2;
  aa2[2] = xdata[indx2];
  aa2[3] = (double)j * arcsec;

  ret = gauss_fit(&xd[0], &yd[0], sig, j, aa2, isSet(PF_V_FLAG));

  if(!ret)
  {
    st_fail("2nd Half Gauss Fit Failed");

    for(i=0;i<j;i++)
    {
      st_print("%2d %4.1f %4.1f\n", i, xd[i], yd[i]);
    }

    return(1);
  }

  st_report("Found Minus Beam @ %7.2f\", Peak = %7.2f, Beam Width = %7.2f",
                aa2[2], aa2[1], fabs(aa2[3]));

  peak2 = aa2[2];
  beam2 = fabs((aa2[2] - xdata[hm]));

  st_report("Plus  Beam  +%7.2f\" From Center", beam2);
  st_report("Minus Beam  -%7.2f\" From Center", beam1);
  st_report("Beam Sep:    %7.2f\"", fabs(peak2-peak1));

  beam_map.gmap_results = fabs(peak2-peak1);
  beam_map.gmap_bw1     = fabs(aa1[3]);
  beam_map.gmap_bw2     = fabs(aa2[3]);
  beam_map.gmap_peak1   = fabs(aa1[1]);
  beam_map.gmap_peak2   = fabs(aa2[1]);

  st_print("Beam Ratio:  %7.2f  ", beam_map.gmap_peak1 / beam_map.gmap_peak2);

  if(beam_map.gmap_peak1 / beam_map.gmap_peak2 > 1.0)
  {
    st_print(" <- Neg Beam Stronger ");
  }
  else
  {
    st_print(" Pos Beam Stronger -> ");
  }

  st_print("\n\n");

  p = &beam_map.bmes[y][0];

  use_global_gdata = 1;
  gridmap_goWriteCenterRow(xdata, ydata, aa1, aa2, beam_map.map_w, p->headers[0]->scan);

  return(0);
}



int gridmap_lsqfit(x, y, bdata, ndata, sig, mwt, a, b, siga, sigb, chi2, q)
double x[], y[], sig[], *a, *b, *siga, *sigb, *chi2, *q;
int bdata, ndata, mwt;
{
  int    i;
  double wt, t, sxoss, sx=0.0, sy=0.0, st2=0.0, ss, sigdat;

  *b = 0.0;

  if(mwt)
  {
    ss = 0.0;

    for(i=bdata;i<ndata;i++)
    {
      wt  = 1.0 / SQR(sig[i]);
      ss += wt;
      sx += x[i] * wt;
      sy += y[i] * wt;
    }
  }
  else
  {
    for(i=bdata;i<ndata;i++)
    {
      sx += x[i];
      sy += y[i];
    }

    ss = ndata - bdata;
  }

  sxoss = sx / ss;

  if(mwt)
  {
    for(i=bdata;i<ndata;i++)
    {
      t    = (x[i] - sxoss) / sig[i];
      st2 += t * t;
      *b  += t * y[i] / sig[i];
    }
  }
  else
  {
    for(i=bdata;i<ndata;i++)
    {
      t    = x[i] - sxoss;
      st2 += t * t;
      *b  += t * y[i];
    }
  }

  *b   /= st2;
  *a    = (sy - sx * (*b)) / ss;
  *siga = sqrt((1.0 + sx*sx / (ss * st2)) / ss);
  *sigb = sqrt(1.0 / st2);
  *chi2 = 0.0;

  if(mwt == 0)
  {
    for(i=bdata;i<ndata;i++)
    {
      *chi2 += SQR(y[i] - (*a) - (*b) * x[i]);
    }

    *q=1.0;
    sigdat = sqrt((*chi2) / (ndata - bdata -2));

    *siga *= sigdat;
    *sigb *= sigdat;
  }

  return(0);
}


int gridmap_findMinMax()
{
  int x, y;
  double val;

  beam_map.map_max = -1e9;
  beam_map.map_min =  1e9;

  for(y=0;y<beam_map.map_h;y++)
  {
    for(x=0;x<beam_map.map_w;x++)
    {
      val = beam_map.bmes[y][x].pval;

      if(val > beam_map.map_max)
      {
        beam_map.map_max = val;
      }
      if(val < beam_map.map_min)
      {
        beam_map.map_min = val;
      }
    }
  }

  beam_map.map_min -= 0.5;
  beam_map.map_max += 0.5;

  if(beam_map.map_min < env->ylower)
  {
    st_report("Truncating Z-Scale upper to %f", env->yupper);
    beam_map.map_min = env->ylower;
  }

  if(beam_map.map_max > env->yupper)
  {
    st_report("Truncating Z-Scale lower to %f", env->ylower);
    beam_map.map_max = env->yupper;
  }

  if(verboseLevel & VERBOSE_GRIDMAP)
  {
    st_report("Plot: Min/Max Map Values: %f/%f", beam_map.map_min, beam_map.map_max);
  }

  return(0);
}


/* Remove a baseline from the Raw Array */
int gridmap_removeBaselines(nb)
int nb;
{
  int    x, y, j, start, end;
  double xd[32], yd[32], sig[512];
  double a, b, siga, sigb, chi2, q, xdata[512];

  if(beam_map.base_removed)
  {
    st_report("Baselines already removed from Map");

    return(0);
  }

  if(nb == 0)
  {
    st_report("mapbaselines set to 0; skipping removal");

    return(0);
  }

  if(nb < 0)
  {
    st_report("gridmap_removeBaselines(): Removing Linear Baseline Using Center %d samples", abs(nb));
  }
  else
  {
    st_report("gridmap_removeBaselines(): Removing Linear Baseline Using %d samples on both ends", nb);
  }

  bzero((char *)&xd[0], sizeof(xd));
  bzero((char *)&yd[0], sizeof(yd));

  for(x=0;x<beam_map.map_w;x++) /* Create the x-axis data */
  {
    xdata[x] = (float)x - (float)beam_map.map_w / 2.0;
  }

  for(y=0;y<beam_map.map_h;y++) /* For each row */
  {
    j = 0;

    if(nb < 0)
    {
      start = (beam_map.map_w / 2) - (abs(nb) / 2);
      end   = (beam_map.map_w / 2) + (abs(nb) / 2);

      for(x=start;x<end;x++)
      {
        xd[j] = xdata[x];
        yd[j] = beam_map.bmes[y][x].rval;

        j++;
      }
    }
    else
    {
      for(x=0;x<nb;x++)
      {
        xd[j] = xdata[x];
        yd[j] = beam_map.bmes[y][x].rval;

        j++;
      }

      for(x=beam_map.map_w-nb;x<beam_map.map_w;x++)
      {
        xd[j] = xdata[x];
        yd[j] = beam_map.bmes[y][x].rval;

        j++;
      }
    }

    gridmap_lsqfit(xd, yd, 0, j-1, sig, 0, &a, &b, &siga, &sigb, &chi2, &q);

    if(a != 0.0 && b != 0.0) /* Remove a baseline where there is data */
    {
      for(x=0;x<beam_map.map_w;x++)
      {
        beam_map.bmes[y][x].rval = beam_map.bmes[y][x].rval -(a+b*xdata[x]);
      }
    }
  }

  beam_map.base_removed = 1;
  beam_map.nbase = abs(nb);

  return(0);
}


int scanAlreadyLoaded(p, scan)
struct BEAM_MAP_ELEM *p;
double scan;
{
  int i;

  for(i=0;i<p->n;i++)
  {
    if(fabs(p->headers[i]->scan-scan) < 0.005)
    {
      return(1);
    }
  }

//  p->scans[p->n] = scan;

  return(0);
}



int gmapFillOffsets(h)
struct HEADER *h;
{
  int    x, y;
  struct BEAM_MAP_ELEM *p;

  for(y=0;y<beam_map.map_h;y++)
  {
    p = &beam_map.bmes[y][0];

    for(x=0;x<beam_map.map_w;x++,p++)
    {
      p->azoff = ((float)x - (float)beam_map.map_w / 2.0)       * h->deltaxr;
      p->eloff = ((float)y - (float)beam_map.map_h / 2.0 + 0.5) * h->deltayr;
    }
  }

  return(0);
}


/* Grid the data into the Raw Array */
int gridmap_grid_data(data, h)
float *data;
struct HEADER *h;
{
  int    x, y, dlen, skip = 0;
  float  ydata[MAX_CONT_SAMPLES];
  double iwt, owt, beamthrow;
  struct BEAM_MAP_ELEM *p;
  struct BEAM_MAP *bm;

  if(processContinuum(h, data, &ydata[0], &dlen, "Az-El Grid Map"))
  {
    return(1);
  }

  bm = &beam_map;

  if(!bm->offsets_filled) /* the first time through ?? */
  {
    bzero(bm->source, sizeof(beam_map.source));
    strncpy(bm->source, h->object, 16);
    strcmprs(bm->source);

    bm->freq = h->restfreq;
    bm->beamthrow = fabs(h->mbeam[0]) + fabs(h->pbeam[0]);

    gmapFillOffsets(h);

    bm->offsets_filled = 1;

    st_report("Setting Beam Map to: Source: %s, Freq %f, Beam Throw: %f", bm->source, bm->freq, bm->beamthrow);
  }

  /* Check to make sure this scan matches the map cube */

  if(strncmp(h->object, bm->source, strlen(bm->source)))
  {
    st_fail("Scan %.2f, not same Source: %.16s as map: %f", h->scan, h->object, bm->freq);
    return(1);
  }

  if(fabs(h->restfreq-bm->freq) > 1.0)
  {
    st_fail("Scan %.2f, not same freq: %f as map: %f", h->scan, h->restfreq, bm->freq);
    return(1);
  }

  beamthrow = fabs(h->mbeam[0]) + fabs(h->pbeam[0]);

  if(fabs(beamthrow-bm->beamthrow) > 2.0)
  {
    st_fail("Scan %.2f, not same BeamThrow: %f as map: %f", h->scan, beamthrow, bm->beamthrow);
    return(1);
  }

  y = (int)h->ycell0 - 1;
  st_report("Gridding %.2f; Grid-Row %2d", h->scan, y);

  if(y < 0 || y >= bm->map_h)
  {
    st_fail("Offset %d out of range of map %d", y, bm->map_h);

    return(1);
  }

  iwt = h->samprat / (h->stsys * h->stsys);     /* New data */

  p = &bm->bmes[y][0];


  st_report("gridmap_grid_data(): Gridding %d samples", dlen);

  for(x=0;x<dlen;x++,p++)
  {
    if(p->n < MAX_GMAP_SCANS)
    {
      p->headers[p->n] = do_addHeader(0, h->scan, h);

      if(!scanAlreadyLoaded(p, h->scan))
      {
        p->x      = x;
        p->y      = y;
        p->tsys  *= p->itime; /* Un-roll previous tys */

        p->itime += h->samprat;
        p->tsys  += (h->stsys * h->samprat);
        p->tsys  /= p->itime;

//        p->azoff = ((float)x - (float)bm->map_w / 2.0)       * h->deltaxr;
        p->eloff = ((float)y - (float)bm->map_h / 2.0 + 0.5) * h->deltayr;

    /* Compute weights and weight-average data */
        owt = p->wt;                            /* Data already in array */

        p->rval  = (p->rval * owt + ydata[x] * iwt) / (owt + iwt);      /* Average them using a weighted scale */

        p->wt += iwt;

        p->n++;
      }
      else
      {
        skip++;
      }
    }
    else
    {
      st_print("Too many scans for array: limit: %d", MAX_GMAP_SCANS);
    }
  }

  if(skip)
  {
    st_report("Skipping scan %.2f; already loaded", h->scan);
    return(1);
  }

  return(0);
}


/* Limit value to upper and lower limits */
double limitValue(l, v, u)
double l, v, u;
{
  double val;

  if(v > u)
  {
    val = u;
  }
  else
  if(v < l)
  {
    val = l;
  }
  else
  {
    val = v;
  }

  return(val);
}



/* Copy the Raw Array to the Plot Array appling any limits, checking for NAN's and
 * conversion to dB's if requested 
 */
int prepareMap()
{
  int x, y, dBmap = 0, isnanCount=0;
  double pvalue, max_value = -1e9, min_value = 1e9, ylower, yupper;

  if(isSet(PF_D_FLAG))
  {
    freeY();                    /* Can't use upper and lower limits and decibel scale */
    dBmap = 1;
  }

  yupper = env->yupper;
  ylower = env->ylower;

        /* Find map Max, Min, check for NAN values */
  for(y=0;y<beam_map.map_h;y++)
  {
    for(x=0;x<beam_map.map_w;x++)
    {
      if(gsl_isnan(beam_map.bmes[y][x].rval))
      {
        beam_map.bmes[y][x].rval = 0.0; /* Correct it */
        isnanCount++;
      }

      pvalue = beam_map.bmes[y][x].rval;

      max_value = GSL_MAX(pvalue, max_value);
      min_value = GSL_MIN(pvalue, min_value);
    }

    if(isnanCount)
    {
      st_report("Corrected %d bad cells in Map", isnanCount);
    }
  }

  if(isSet(PF_V_FLAG))
  {
    st_report("Raw:  Min/Max Map Values: %f/%f", min_value, max_value);
  }

  if(dBmap)
  {
    strcpy(beam_map.yaxis_label, "dB's");
  }
  else
  {
    strcpy(beam_map.yaxis_label, "Intensity (K)");
  }

  for(y=0;y<beam_map.map_h;y++)
  {
    for(x=0;x<beam_map.map_w;x++)
    {
      pvalue = limitValue(ylower, beam_map.bmes[y][x].rval, yupper);

      if(dBmap)
      {
        if(env->lowerdb > 0.0 && fabs(pvalue) < env->lowerdb)
        {
          pvalue = env->lowerdb;
        }

        if(pvalue != 0.0)
        {
          pvalue = log10(fabs(pvalue) / max_value) * 10.0;
        }
        else
        {
          pvalue = env->lowerdb;
        }

        if(env->lowerdb < 0.0 && pvalue < env->lowerdb)
        {
          pvalue = env->lowerdb;
        }

      }

      beam_map.bmes[y][x].pval = pvalue;
    }
  }

  return(0);
}


int findMapColIndx(x1, y1)
double x1;
int    y1;
{
  int x, n, indx = 0;
  double diff, closest = 1e9;

  if(verboseLevel & VERBOSE_GRIDMAP)
  {
    st_report("findMapColIndx(%f, %d)", x1, y1);
  }

  n = beam_map.map_w;

  for(x=0;x<n;x++)
  {
    diff = fabs(beam_map.bmes[y1][x].azoff - x1);

    if(verboseLevel & VERBOSE_GRIDMAP)
    {
      st_report("Checking fabs(%f - %f): %f", beam_map.bmes[y1][x].azoff, x1, diff);
    }

    if(diff < closest)
    {
      closest = diff;
      indx = x;
    }
  }

  if(verboseLevel & VERBOSE_GRIDMAP)
  {
    st_report("findMapColIndx(): Closest = %f, indx = %d", closest, indx);
  }

  return(indx);
}


int findMapRowIndx(x1, y1)
int    x1;
double y1;
{
  int y, n, indx = 0;
  double diff, closest = 1e9;

  if(verboseLevel & VERBOSE_GRIDMAP)
  {
    st_report("findMapRowIndx(%d, %f)", x1, y1);
  }

  n = beam_map.map_h;

  for(y=0;y<n;y++)
  {
    diff = fabs(beam_map.bmes[y][x1].eloff - y1);

  if(verboseLevel & VERBOSE_GRIDMAP)
  {
    st_report("Checking fabs(%f - %f): %f", beam_map.bmes[y][x1].eloff, y1, diff);
  }

    if(diff < closest)
    {
      closest = diff;
      indx = y;
    }
  }

  if(verboseLevel & VERBOSE_GRIDMAP)
  {
    st_report("findMapRowIndx(): Closest = %f, indx = %d", closest, indx);
  }

  return(indx);
}


int gridMapSlice()
{
  int x, y, cr, start, end, dlen, old_plotstyle;
  double xb, xe, yb, ye;
  struct BEAM_MAP_ELEM *p;

  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

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

    return(1);
  }

  if(beam_map.map_w == 0)
  {
    st_fail("No Grid Map Loaded");

    return(1);
  }

  cr = atoi(args->secondArg); /* Col or row number */

  cr--;   /* Make it 0 based */

  if(isSet(PF_Y_FLAG)) /* DO columns ?? */
  {
    if(cr < 0 || cr > beam_map.map_w)
    {
      st_fail("Col number %d, out of range: %d to %d", cr, 0, beam_map.map_w);
      return(1);
    }
  }
  else
  {
    if(cr < 0 || cr > beam_map.map_h)
    {
      st_fail("Row number %d, out of range: %d to %d", cr, 0, beam_map.map_h);
      return(1);
    }
  }

  if(args->ntokens > 3) /* User selected a sub range */
  {
    if(isSet(PF_Y_FLAG)) /* DO columns ?? */
    {
      yb = atof(args->thirdArg);
      ye = atof(args->fourthArg);

      start = findMapRowIndx(cr, yb);
      end   = findMapRowIndx(cr, ye);

      dlen  = (end - start);
    }
    else
    {
      xb = atof(args->thirdArg);
      xe = atof(args->fourthArg);

      start = findMapColIndx(xb, cr);
      end   = findMapColIndx(xe, cr);

      dlen  = (end - start);
    }
  }
  else
  {
    start = 0;

    if(isSet(PF_Y_FLAG)) /* DO columns ?? */
    {
      end   = beam_map.map_h;
    }
    else
    {
      end   = beam_map.map_w;
    }

    dlen  = end - start;
  }

  if(verboseLevel & VERBOSE_GRIDMAP)
  {
    st_report("Start = %d, End = %d, Dlen = %d", start, end, dlen);
  }

  if(isSet(PF_Y_FLAG)) /* DO columns ?? */
  {
    /* Pick center row for header */

    p = &beam_map.bmes[beam_map.map_h/2][0];

    if(p->headers[0])
    {
      clearFillContStruct(p->headers[0], dlen, "Map El Offset", beam_map.yaxis_label);
    }
    else
    {
      st_fail("gmapSlice(): No header found");
      return(1);
    }

    for(y=start;y<end;y++)
    {
      cont_plain_send.xdata[y-start] = beam_map.bmes[y][cr].eloff;
      cont_plain_send.ydata[y-start] = beam_map.bmes[y][cr].pval;
    }

    if(xgraphic)
    {
      sendCont(SEND_PLAIN);
    }
  }
  else
  {
    p = &beam_map.bmes[cr][0];

    if(p->headers[0])
    {
      clearFillContStruct(p->headers[0], dlen, "Map Az Offset", beam_map.yaxis_label);
    }
    else
    {
      st_fail("gmapSlice(): No header found");
      return(1);
    }

    for(x=start;x<end;x++)
    {
      cont_plain_send.xdata[x-start] = beam_map.bmes[cr][x].azoff;
      cont_plain_send.ydata[x-start] = beam_map.bmes[cr][x].pval;
    }

    if(use_global_gdata)
    {
      for(x=start;x<end;x++)
      {
        cont_plain_send.gdata[x] = global_gdata[x];
      }

      cont_plain_send.ngauss = beam_map.map_w;
    }

    if(xgraphic)
    {
      old_plotstyle = env->plotstyle;

      env->plotstyle = PLOT_STYLE_LINES;

      sendCont(SEND_PLAIN);

      env->plotstyle = old_plotstyle;
    }
  }

  return(0);
}



int gridMapPlot()
{
  int win;
  char name[256], sysBuf[256];

  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  if(notSet(PF_F_FLAG)) /* Auto set map angle or leave it fixed */
  {
    doGetMapAngle();
  }

  if(isSet(PF_0_FLAG))
  {
    env->map_angle[0] =  60.0;
    env->map_angle[1] =  30.0;
  }
  else
  if(isSet(PF_1_FLAG))
  {
    env->map_angle[0] = 180.0;
    env->map_angle[1] =  90.0;
  }
  else
  if(isSet(PF_2_FLAG))
  {
    env->map_angle[0] =  90.0;
    env->map_angle[1] =   0.0;
  }
  else
  if(isSet(PF_3_FLAG))
  {
    env->map_angle[0] =   0.0;
    env->map_angle[1] =   0.0;
  }
  else
  if(isSet(PF_4_FLAG))
  {
    env->map_angle[0] =   0.0;
    env->map_angle[1] =  90.0;
  }
  else
  if(isSet(PF_5_FLAG))
  {
    env->map_angle[0] =  90.0;
    env->map_angle[1] = 270.0;
  }
  else
  if(isSet(PF_6_FLAG))
  {
    env->map_angle[0] =  90.0;
    env->map_angle[1] = 180.0;
  }
  else
  if(isSet(PF_7_FLAG))
  {
    env->map_angle[0] =  60.0;
    env->map_angle[1] = 220.0;
  }
  else
  if(isSet(PF_8_FLAG))
  {
    env->map_angle[0] =  90.0;
    env->map_angle[1] =  90.0;
  }
  else
  if(isSet(PF_9_FLAG))
  {
    env->map_angle[0] = 270.0;
    env->map_angle[1] = 270.0;
  }


  /* Remember the flags for any move that's played later */
  bcopy(args->pflags, (char *)&gmap_pflags, sizeof(gmap_pflags));

  prepareMap();         /* Go prepare the map for plotting */
  gridmap_findMinMax(); /* Determine it's bounds */

  if(isSet(PF_C_FLAG))
  {
    gridmap_write_csvfile();
    validate(CONTMAP);

    return(0);
  }

  if(isSet(PF_P_FLAG) == 1) /* Make a PS file */
  {
    doMapPlot(PSPLOT);
  }
  else
  {
    doMapPlot(GNUPLOT); /* Display a GUI plot */
  }

  if(isSet(PF_P_FLAG) == 2) /* Make a PNG file */
  {
    if(doingDive)
    {
      win = getPlotWinId(21);
    }
    else
    {
      win = getPlotWinId(20);
    }

    if(!win)
    {
      st_fail("gridMapPlot(): Unable to determine GNUPLOT windo ID");

      validate(CONTMAP);

      return(1);
    }

    sprintf(name, "%s.%d.png", env->source, (int)(env->if_feed*100.0));

    sprintf(sysBuf, "import -window 0x%x -silent %s", win, name);
    system(sysBuf);

    st_report("Created file: %s", name);
  }

  validate(CONTMAP);

  return(0);
}

int doMapPlot(how)
int how;
{
  char sysBuf[256];

  if(beam_map.map_h > 0 && beam_map.map_w > 0)
  {
    gridmap_write_datafile();
    gridmap_writeDemFile(beam_map.map_h / 2.0, how);

    sprintf(sysBuf, "xgraphic load '/tmp/%s/bmap.dem'", session);

    if(verboseLevel & VERBOSE_GRIDMAP)
    {
      st_report("Parsing Command: %s", sysBuf);
    }

    parseCmd(sysBuf);
  }
  else
  {
    st_warn("No Map to Plot");
    st_warn("Have you done a 'gridmap'?");
  }

  return(0);
}


int mapGaussPlot()
{
  gridmap_showMiddleRow();

  return(0);
}


int mapGaussFit()
{
  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  gridmap_fitGauss();

  return(0);
}


int mapBaseline()
{
  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  gridmap_removeBaselines(env->mapbaselines);

  gridmap_findMinMax();

  return(0);
}


int gridMapGrid()
{
  int i, n=0, ret, delay=0, first = 1, w, h, append=0;
  struct HEADER *p;

  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  if(!env->numofscans)
  {
    st_warn("No data file loaded");

    return(1);
  }

  if(notSet(PF_A_FLAG))                 /* Do not zero out beam_map struct if we are appending */
  {
    bzero((char *)&beam_map, sizeof(beam_map));

    remScanHeaders();                   /* Zero out scan Headers */
  }
  else
  {
    beam_map.base_removed = 0;          /* Allow baselines to be removed again */
    append = 1;
  }

  if(isSet(PF_D_FLAG))
  {
    delay = 1;
  }

  p = &sa[0];

  for(i=0;i<env->numofscans;i++,p++)
  {
    if(!meetsCritiria(p, CHECK_ALL))
    {
      continue;
    }

    if(strncmp(p->obsmode, "CONTMAP", 7))
    {
      continue;
    }

    if(isScanIgnored(p->scan))
    {
      continue;
    }

    if((ret = read_scan( p->scan, env->filename, h0->yplot1, &h0->head, 0)) > 0.0)
    {
      env->curScan = p->scan;

      w = (int)h0->head.noxpts;
      h = (int)h0->head.noypts;

      if(w == 0 || h == 0 || w > MAX_MAP_W || h > MAX_MAP_H)
      {
        fprintf(stderr, "Error in grid map size: %d %d\n", w, h);

        return(1);
      }

      if(first)
      {
        beam_map.map_w = w;
        beam_map.map_h = h;
        beam_map.map_step_x = h0->head.deltaxr;
        beam_map.map_step_y = h0->head.deltayr;
        beam_map.map_size_x = beam_map.map_w * beam_map.map_step_x;
        beam_map.map_size_y = beam_map.map_h * beam_map.map_step_y;

        st_report("Load Grid Map: %d X %d", beam_map.map_w, beam_map.map_h);

        /* Set plot extents */
        beam_map.plot_xb = 0;
        beam_map.plot_xe = beam_map.map_w;
        beam_map.plot_yb = 0;
        beam_map.plot_ye = beam_map.map_h;

        if(!append)
        {
          beam_map.fscan = h0->head.scan;
        }

        first = 0;
      }

      if(w != beam_map.map_w || h != beam_map.map_h)
      {
        st_fail("Scan %.2f, has incorrect gridmap size: %d x %d", h0->head.scan, w, h);

        return(1);
      }

      showScan();

      if(delay)
      {
        usleep(env->loopdelay1);
      }

      if(!gridmap_grid_data(h0->yplot1, &h0->head))
      {
        n++;
        beam_map.nscans++;

        beam_map.lscan = h0->head.scan;
      }
    }
    else
    {
      st_fail("Unable to read scan %.2f", p->scan);
    }
  }

  if(n)
  {
    st_report("Found %d Scan(s) out of %d Matching Criteria", n, env->numofscans);

    if(notSet(PF_N_FLAG))
    {
      mapBaseline();
    }
    else
    {
      st_report("Skipping baseline removal");
    }

    validate(CONT_MAP);
  }
  else
  {
    st_report("No Scans Matched Criteria or all already included");
  }

  return(0);
}


int gridMapGrab()
{
  char sysBuf[256];

  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  sprintf(sysBuf, "xgraphic load '/tmp/%s/bmap.dem'", session);

  parseCmd(sysBuf);

  return(0);
}



int mBaselines()
{
  int bl;

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

    if(bl > -10 && bl < 10)
    {
      if(verboseLevel & VERBOSE_GRIDMAP)
      {
        st_report("Setting Map Baseline Region to %d samples", bl);
      }

      env->mapbaselines = bl;
    }
    else
    {
      st_warn("Map Baselines of %d out of range: -10 to +10", bl);
    }
  }

  return(0);
}


int setMapAngle()
{
  if(!args->secondArg || !args->thirdArg)
  {
    st_warn("Missing Args for mapangle");
  }

  env->map_angle[0] = atof(args->secondArg);
  env->map_angle[1] = atof(args->thirdArg);

  if(verboseLevel & VERBOSE_GRIDMAP)
  {
    st_report("Setting Map Angle to %.1f %.1f", env->map_angle[0], env->map_angle[1]);
  }

  return(0);

}


int centerRow()
{
  double scan;
  double bscan, escan;

  if(!args->secondArg)
  {
    st_warn("Missing Args for center_row()");

    return(1);
  }

  beam_map.gmap_bw1 = beam_map.gmap_bw2 = beam_map.gmap_peak1 = beam_map.gmap_peak2 = 0.0;

  bscan = env->bscan;
  escan = env->escan;

  scan = floor(atof(args->secondArg));

  if(verboseLevel & VERBOSE_GRIDMAP)
  {
    st_report("Analysing Center Map Row %.0f", scan);
  }

  env->bscan = scan;
  env->escan = scan;

  gridMapGrid();

  mapGaussFit();

  env->bscan = bscan;
  env->escan = escan;

  return(0);
}


int centerAll()
{
  int i, n=0, pause = 0, ret, first = 1;
  struct HEADER *p;
  double bscan, escan, err, error;
  double bsep, bw1t, bw2t, peak1, peak2, ratio;
  struct BEAM_MAPS_SEP *bms;
  char *symbol;

  if(!env->numofscans)
  {
    st_warn("No data file loaded");

    return(1);
  }

  if(getYesNoQuit("Pause Between Scans: y/N ?"))
  {
    pause = 1;
  }

  bscan = env->bscan;
  escan = env->escan;

  p = &sa[0];
  bms = &beam_map_sep[0];

  for(i=0;i<env->numofscans;i++,p++)
  {
    if(!meetsCritiria(p, CHECK_ALL))
    {
      continue;
    }

    if(strncmp(p->obsmode, "CONTMAP", 7))
    {
      continue;
    }

    if(isScanIgnored(p->scan))
    {
      continue;
    }

    if(p->noypts > 1)
    {
      continue;
    }
    
    subArgs(2, "%.0f", p->scan);

    if(first)
    {
      clearScreen();
      first = 0;
    }

    gotoxy(1, 1);

    centerRow();

    if(pause)
    {
      mapGaussPlot();

      ret = getYesNoQuit("Any key to continue; 'q' to quit: [*,q]?");

      if(ret == 2)
      {
        st_report("Aborting Plots");
        pause = 0;
      }
    }

    bms->scan       = p->scan;
    bms->seperation = beam_map.gmap_results;
    bms->bw1        = beam_map.gmap_bw1;
    bms->bw2        = beam_map.gmap_bw2;
    bms->peak1      = beam_map.gmap_peak1;
    bms->peak2      = beam_map.gmap_peak2;
    bms->ratio      = bms->peak1 / bms->peak2;
    
    /* the above functions messes with the bscan and escan; reset them here */
    env->bscan = bscan;
    env->escan = escan;

    n++;
    bms++;
  }

  if(n)
  {
    st_report("Found %d Scan(s) out of %d Matching Criteria", n, env->numofscans);

    st_report("Results of %d map scans", n);
    st_report("");

    st_print(" Scan#    Seperation  Error    BW1    BW2      PK1    PK2     Ratio\n");

    error = 0.0;
    bw1t = bw2t = 0.0;
    peak1 = peak2 = ratio = 0.0;

    bms = &beam_map_sep[0];
    for(i=0;i<n;i++,bms++)
    {
      bw1t  += bms->bw1;
      bw2t  += bms->bw2;
      peak1 += bms->peak1;
      peak2 += bms->peak2;
      ratio += bms->ratio;
    }

    bw1t  /= (double)n;
    bw2t  /= (double)n;
    peak1 /= (double)n;
    peak2 /= (double)n;
    ratio /= (double)n;

    bsep = fabs(h0->head.pbeam[0] - h0->head.mbeam[0]);

    bms = &beam_map_sep[0];
    for(i=0;i<n;i++,bms++)
    {
      err = bms->seperation - bsep;

      if(bms->ratio < 0.995)
      {
        symbol = "->";
      }
      else
      if(bms->ratio > 1.005)
      {
        symbol = "<-";
      }
      else
      {
        symbol = "==";
      }

      st_print("%.2f      %.2f  %6.2f   %5.2f  %5.2f   %5.2f   %5.2f    %5.2f %s\n", 
                bms->scan, 
                bms->seperation,
                err,
                bms->bw1,
                bms->bw2,
                bms->peak1,
                bms->peak2,
                bms->ratio, 
                symbol);

      error += err;
    }

    error /= (double)n;

    st_print("---------------------------\n");

    st_print(" %3d Scans:  %.2f %7.3f   %5.2f  %5.2f   %5.2f   %5.2f    %5.2f\n\n",
                                n, bsep, error, bw1t, bw2t, peak1, peak2, ratio);
  }
  else
  {
    st_fail("No Scans Matched Criteria");
  }

  return(0);
}


int showMidRow()
{
  double scan;

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

    sprintf(middleRowFname, "/tmp/%s/middle_row.%04.2f.dat", session, scan);

    gridmap_showMiddleRow();
  }

  return(0);
}


int cont_removeBaselines(d, len, scan)
float *d;
int len;
double scan;
{
  int    i, j=0, nb, start, end;
  double xd[32], yd[32], sig[32];
  double a, b, siga, sigb, chi2, q;

  bzero((char *)&xd[0], sizeof(xd));
  bzero((char *)&yd[0], sizeof(yd));

  nb = env->mapbaselines;

  if(abs(nb) < 2)
  {
    st_report("Not enough baseline cells definded: %d", nb);
    return(0);
  }

  st_report("cont_removeBaselines(): Removing Cont Baseline On %d datapoints", len);

  if(nb < 0)
  {
    st_report("cont_removeBaselines(): Removing Linear Baseline Using center %d samples", nb*-1);
  }
  else
  {
    st_report("cont_removeBaselines(): Removing Linear Baseline Using %d samples on both ends", nb);
  }

  if(nb < 0)
  {
    nb *= -1;

    start = beam_map.map_w / 2 - nb / 2;
    end   = beam_map.map_w / 2 + nb / 2;

    for(i=start;i<end;i++)
    {
      xd[j]  = (float)i;
      yd[j]  = d[i];
      sig[j] = 0.0;

      j++;
    }
  }
  else
  {
    for(i=0;i<nb;i++)
    {
      xd[j]  = (float)i;
      yd[j]  = d[i];
      sig[j] = 0.0;

      j++;
    }

    for(i=len-nb;i<len;i++)
    {
      xd[j]  = (float)i;
      yd[j]  = d[i];
      sig[j] = 0.0;

      j++;
    }
  }

  gridmap_lsqfit(xd, yd, 0, j-1, sig, 0, &a, &b, &siga, &sigb, &chi2, &q);

  if(a != 0.0 && b != 0.0) /* Remove a baseline where there is data */
  {
    for(i=0;i<len;i++)
    {
      d[i] = d[i] -(a+b*(float)i);
    }
  }

  return(0);
}



int getPlotWinId(w)
int w;
{
  int id = 0;
  char pbuf[256], inbuf[256], *cp;
  FILE *fp;

  sprintf(pbuf, "/usr/bin/xwininfo -name 'GNUPLOT_%s %d' | grep 'Window id' | awk '{print $4}'", session, w);

  if( (fp = popen(pbuf, "r") ) <= (FILE *)0 )
  {
    st_fail("getPlotWinId(): Unable to open %s", pbuf);

    return(0);
  }

  bzero(inbuf, sizeof(inbuf));

  while((cp = fgets(inbuf, sizeof(inbuf), fp)) != NULL)
  {
    deNewLineIt(inbuf);

    id = atox(inbuf);
  }

  pclose(fp);

  return(id);
}


int gmapRotate()
{
  int i, win, record = 0, keep = 0;
  double angle;
  char name[256], sysBuf[256];

  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  if(isSet(PF_R_FLAG)) /* Record the rotation ?? */
  {
    record = 1;

    if(isSet(PF_K_FLAG)) /* Keep intermidiate frames ?? */
    {
      keep = 1;
    }
  }

  setFlag(PF_R_FLAG, 0); /* Unset so as not to mess with gridplot command */
  setFlag(PF_K_FLAG, 0); /* Unset so as not to mess with gridplot command */

  if(record)
  {
    sprintf(name, "%s.%d.seq", env->source, (int)(env->if_feed*100.0));

    win = getPlotWinId(20);

    if(!win)
    {
      st_fail("gmapRotate(): Unable to determine GNUPLOT windo ID");

      return(1);
    }

    st_report("gmapRotate(): Using Window ID: 0x%x, filenames: %s", win, name);
  }

  doGetMapAngle();

  angle = env->map_angle[1];  /* We will only rotate the second angle */

  /* Make the equal to the last time gmapplot was called */
  for(i=0;i<MAX_PARSE_FLAGS;i++)
  {
    if((gmap_pflags[i].set))
    {
      setFlag(i, gmap_pflags[i].set);
    }
  }

  if(record) /* Delete any old ones */
  {
    st_report("Deleting any Leftover Intermidiate Frames: /tmp/%s.???.dpx", name);

    sprintf(sysBuf, "rm -f /tmp/%s.???.dpx", name);
    system(sysBuf);
  }


  setFlag(PF_F_FLAG, 1);  /* I want gridMapPlot() to ignore the current map angle and use my new one */

  clearInterrupt();

  for(i=0;i<120;i++)
  {
    angle += 3.0;

    if(angle > 360.0)
    {
      angle -= 360.0;
    }

    env->map_angle[1] = angle;

    st_report("Using Angles: %2d: %f, %f", i, env->map_angle[0], env->map_angle[1]);

    gridMapPlot();

    usleep(200000);

    if(isInterrupt()) /* Check for Control-C input; any will cause reading to stop */
    {
      st_print("Control-C hit; Stoping\n");

      if(record && !keep) /* Delete intermidiate frames */
      {
        st_report("Deleting Intermidiate Frames: /tmp/%s.???.dpx", name);

        sprintf(sysBuf, "rm -f /tmp/%s.???.dpx", name);
        system(sysBuf);
      }

      return(0);
    }

    doGetMapAngle();

    if(record)
    {
      sprintf(sysBuf, "import -window 0x%x -silent /tmp/%s.%03d.dpx", win, name, i);
      system(sysBuf);
    }
  }

  if(record)
  {
    st_report("Deleting any Old Movie File: %s.m4v", name);

    sprintf(sysBuf, "rm -f %s.m4v", name);
    system(sysBuf);

    st_report("Creating New Movie File: %s.m4v", name);

    sprintf(sysBuf, "xterm -bg black -fg white -fn 10x20 -e /usr/bin/ffmpeg -r 8 \
                     -i \"/tmp/%s.%%3d.dpx\" -crf 10 -vcodec h264 %s.m4v; sleep 5 &", name, name);

    strcmprs(sysBuf);

    system(sysBuf);

    if(!keep) /* Delete intermidiate frames */
    {
      st_report("Deleting Intermidiate Frames: /tmp/%s.???.dpx", name);

      sprintf(sysBuf, "rm -f /tmp/%s.???.dpx", name);
      system(sysBuf);
    }
  }

  return(0);
}



int gridMapDive()
{
  int     i, lowerdb, maxdb, win, record = 0, keep = 0, indx = 0;
  double  old_lowerdb, old_mapangle[2], old_zmin, old_zmax;
  char name[256], sysBuf[256];

  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  if(isSet(PF_R_FLAG)) /* Record ?? */
  {
    record = 1;

    if(isSet(PF_K_FLAG)) /* Keep intermidiate frames ?? */
    {
      keep = 1;
    }
  }

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

    return(1);
  }

  maxdb = abs(atoi(args->secondArg));

  old_zmin          = env->gmap_zmin;
  old_zmax          = env->gmap_zmax;
  old_mapangle[0]   = env->map_angle[0];
  old_mapangle[1]   = env->map_angle[1];
  env->map_angle[0] = 180.0;
  env->map_angle[1] =  90.0;
  env->gmap_zmin    = 0.0;
  env->gmap_zmax    = 0.0;

  setFlag(PF_R_FLAG, 0); /* Unset so as not to mess with gridplot command */
  setFlag(PF_K_FLAG, 0); /* Unset so as not to mess with gridplot command */

  doingDive = 1;

  sprintf(name, "%s.%d.seq", env->source, (int)(env->if_feed*100.0));


  old_lowerdb = env->lowerdb;

  /* Make the equal to the last time gmapplot was called */
  for(i=0;i<MAX_PARSE_FLAGS;i++)
  {
    if((gmap_pflags[i].set))
    {
      setFlag(i, gmap_pflags[i].set);
    }
  }

  if(record) /* Delete any old ones */
  {
    st_report("Deleting any Leftover Intermidiate Frames: /tmp/%s.???.dpx", name);

    sprintf(sysBuf, "rm -f /tmp/%s.???.dpx", name);
    system(sysBuf);
  }

  setFlag(PF_F_FLAG, 1);  /* I want gridMapPlot() to ignore the current map angle and use my new one */
  setFlag(PF_D_FLAG, 1);  /* I want gridMapPlot() to plot in dBs */

  indx = 0;

  clearInterrupt();

  for(lowerdb=-1;lowerdb>-maxdb;lowerdb--)
  {
   st_report("Using LowerDB's: %2d", lowerdb);

    env->lowerdb = (double)lowerdb;

    gridMapPlot();

    if(record && !win) /* Can't find window before the first plot */
    {
      win = getPlotWinId(21);

      if(!win)
      {
        st_fail("gmapRotate(): Unable to determine GNUPLOT windo ID");

        doingDive = 0;

        env->map_angle[0] = old_mapangle[0];
        env->map_angle[1] = old_mapangle[1];
        env->gmap_zmin = old_zmin;
        env->gmap_zmax = old_zmax;
  
        return(1);
      }

      st_report("gmapRotate(): Using Window ID: 0x%x, filenames: %s", win, name);
    }

    usleep(200000);

    if(isInterrupt()) /* Check for Control-C input; any will cause reading to stop */
    {
      st_print("Control-C hit; Stoping\n");

      if(record && !keep) /* Delete intermidiate frames */
      {
        st_report("Deleting Intermidiate Frames: /tmp/%s.???.dpx", name);

        sprintf(sysBuf, "rm -f /tmp/%s.???.dpx", name);
        system(sysBuf);
      }

      env->lowerdb      = old_lowerdb;
      env->map_angle[0] = old_mapangle[0];
      env->map_angle[1] = old_mapangle[1];
      env->gmap_zmin    = old_zmin;
      env->gmap_zmax    = old_zmax;
      doingDive         = 0;

      return(0);
    }

    doGetMapAngle();

    if(record)
    {
      sprintf(sysBuf, "import -window 0x%x -silent /tmp/%s.%03d.dpx", win, name, indx);
      system(sysBuf);
    }

    indx++;
  }

  clearInterrupt();

  for(lowerdb=-maxdb;lowerdb<0;lowerdb++)
  {
   st_report("Using LowerDB's: %2d", lowerdb);

    env->lowerdb = (double)lowerdb;

    gridMapPlot();

    usleep(200000);
    if(isInterrupt()) /* Check for Control-C input; any will cause reading to stop */
    {
      st_print("Control-C hit; Stoping\n");

      if(record && !keep) /* Delete intermidiate frames */
      {
        st_report("Deleting Intermidiate Frames: /tmp/%s.???.dpx", name);

        sprintf(sysBuf, "rm -f /tmp/%s.???.dpx", name);
        system(sysBuf);
      }

      env->lowerdb      = old_lowerdb;
      env->map_angle[0] = old_mapangle[0];
      env->map_angle[1] = old_mapangle[1];
      env->gmap_zmin    = old_zmin;
      env->gmap_zmax    = old_zmax;
      doingDive         = 0;

      return(0);
    }

    doGetMapAngle();

    if(record)
    {
      sprintf(sysBuf, "import -window 0x%x -silent /tmp/%s.%03d.dpx", win, name, indx);
      system(sysBuf);
    }

    indx++;
  }

  if(record)
  {
    st_report("Deleting any Old Movie File: %s.m4v", name);

    sprintf(sysBuf, "rm -f %s.m4v", name);
    system(sysBuf);

    st_report("Creating New Movie File: %s.m4v", name);

    sprintf(sysBuf, "xterm -bg black -fg white -fn 10x20 -e /usr/bin/ffmpeg -r 8 \
                     -i \"/tmp/%s.%%3d.dpx\" -crf 10 -vcodec h264 %s.m4v; sleep 5 &", name, name);


    strcmprs(sysBuf);

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

    system(sysBuf);

    if(!keep) /* Delete intermidiate frames */
    {
      st_report("Deleting Intermidiate Frames: /tmp/%s.???.dpx", name);

      sprintf(sysBuf, "rm -f /tmp/%s.???.dpx", name);
      system(sysBuf);
    }
  }

  env->lowerdb      = old_lowerdb;
  env->map_angle[0] = old_mapangle[0];
  env->map_angle[1] = old_mapangle[1];
  env->gmap_zmin    = old_zmin;
  env->gmap_zmax    = old_zmax;
  doingDive         = 0;

  return(0);
}



int contMapFitGauss()
{
  int i, j, x, n, indx1, indx2, ret, hm, glen;
  double xdata[512], ydata[512];
  double arcsec, max1=0.0, max2=0.0;
  double xd[512], yd[512], sig[512];
  double aa1[4], aa2[4], peak1, peak2, beam1, beam2;
  double g1, g2, y, dyda[4];

  glen = 5;

  bzero((char *)&xd[0], sizeof(xd));
  bzero((char *)&yd[0], sizeof(yd));

  hm = cont_plain_send.dlen / 2;  /* Half row */

  arcsec = cont_plain_send.head.deltaxr;

  st_report( "Fitting Gaussian Using Scan(s) %.2f; Arcsec = %.1f", cont_plain_send.head.scan, arcsec);

  n = cont_plain_send.dlen;

  for(x=0;x<n;x++) /* Create the x,y-axis data; Find Min/Max too */
  {
    xdata[x] = cont_plain_send.xdata[x];
    ydata[x] = cont_plain_send.ydata[x];

    sig[x] = 1.0;

    if(ydata[x] > max1)
    {
      max1  = ydata[x];
      indx1 = x;
    }
    if(ydata[x] < max2)
    {
      max2  = ydata[x];
      indx2 = x;
    }
  }

  // Set up for gausian fitting

  if(max1 == 0.0 || max2 == 0.0)
  {
    st_fail("No data to do gausfit");
    return(1);
  }

  st_report("Guess, Plus  Peak: %7.2f @ Sample %d or %.2f\"", max1, indx1, xdata[indx1]);
  st_report("Guess, Minus Peak: %7.2f @ Sample %d or %.2f\"", max2, indx2, xdata[indx2]);
  st_report("Guess, Beam Sep:    %7.2f\"", fabs(indx1-indx2) * arcsec);

  /* Compute ideal gaussian 1/2 power width in samples */
  glen = (int)(cont_plain_send.head.bfwhm / arcsec) + 1;

// Load first peak
  for(i=indx1-glen,j=0;i<indx1+glen;i++,j++)
  {
    xd[j] = xdata[i];
    yd[j] = ydata[i];
  }

  st_report("Fitting a %d Point Gauss to Plus Peak", j+1);
  aa1[1] = max1;
  aa1[2] = xdata[indx1];
  aa1[3] = (double)j * arcsec; /* Beam width */

  ret = gauss_fit(&xd[0], &yd[0], sig, j, aa1, isSet(PF_V_FLAG));

  if(!ret)
  {
    st_fail("1st Half Gauss Fit Failed");

    for(i=0;i<j;i++)
    {
      st_print("%2d %4.1f %4.1f\n", i, xd[i], yd[i]);
    }

    return(1);
  }

  st_report("Found Plus  Beam @ %7.2f\", Peak = %7.2f, Beam Width = %7.2f", aa1[2], aa1[1], fabs(aa1[3]));

  peak1 = aa1[2];
  beam1 = fabs((aa1[2] - xdata[hm]));

// Load second peak
  for(i=indx2-glen,j=0;i<indx2+glen;i++,j++)
  {
    xd[j] = xdata[i];
    yd[j] = ydata[i];
  }

  st_report("Fitting a %d Point Gauss to Minus Peak", j+1);
  aa2[1] = max2;
  aa2[2] = xdata[indx2];
  aa2[3] = (double)j * arcsec;

  ret = gauss_fit(&xd[0], &yd[0], sig, j, aa2, isSet(PF_V_FLAG));

  if(!ret)
  {
    st_fail("2nd Half Gauss Fit Failed");

    for(i=0;i<j;i++)
    {
      st_print("%2d %4.1f %4.1f\n", i, xd[i], yd[i]);
    }

    return(1);
  }

  st_report("Found Minus Beam @ %7.2f\", Peak = %7.2f, Beam Width = %7.2f", aa2[2], aa2[1], fabs(aa2[3]));

  peak2 = aa2[2];
  beam2 = fabs((aa2[2] - xdata[hm]));

  st_report("Plus  Beam  +%7.2f\" From Center", beam2);
  st_report("Minus Beam  -%7.2f\" From Center", beam1);
  st_report("Beam Sep:    %7.2f\"", fabs(peak2-peak1));

  sprintf(cont_plain_send.glabels[0], "Plus  Beam  +%7.2f Arc-Seconds From Center", beam2);
  sprintf(cont_plain_send.glabels[1], "Minus Beam  -%7.2f Arc-Seconds From Center", beam1);
  sprintf(cont_plain_send.glabels[2], "Beam Sep:    %7.2f Arc-Seconds", fabs(peak2-peak1));

  /* Populate the gdata array */
  for(x=0;x<n;x++)
  {
    fgauss(cont_plain_send.xdata[x], aa1, &y, dyda);
    g1 = y;

    fgauss(cont_plain_send.xdata[x], aa2, &y, dyda);
    g2 = y;

    cont_plain_send.gdata[x] = g1 + g2;
  }

  cont_plain_send.ngauss = x;

  if(xgraphic)
  {
    sendCont(SEND_PLAIN);
  }

  return(0);
}


/* Remove a baseline from the loaded CONTMAP scane*/
int removeContMapBaselines(pf, nb)
int pf, nb;
{
  int    x, j, n, start, end;
  double xd[MAX_CONT_SAMPLES], yd[MAX_CONT_SAMPLES], sig[MAX_CONT_SAMPLES];
  double a, b, siga, sigb, chi2, q;

  if(!pf)
  {
    return(0);
  }

  if(nb == 0)
  {
    st_report("mapbaselines set to 0; skipping removal");

    return(0);
  }

  if(nb < 0)
  {
    st_report("removeContMapBaselines(): Removing Linear Baseline Using Center %d samples", abs(nb));
  }
  if(pf == 2)
  {
    st_report("removeContMapBaselines(): Removing Linear Baseline Using %d samples on both ends + middle", nb);
  }
  else
  {
    st_report("removeContMapBaselines(): Removing Linear Baseline Using %d samples on both ends", nb);
  }

  bzero((char *)&xd[0], sizeof(xd));
  bzero((char *)&yd[0], sizeof(yd));

  n = cont_plain_send.dlen;

  j = 0;

  if(nb < 0)                            /* Take them from the middle */
  {
    start = (n / 2) - (abs(nb) / 2);
    end   = (n / 2) + (abs(nb) / 2);

    for(x=start;x<end;x++)
    {
      xd[j] = cont_plain_send.xdata[x];
      yd[j] = cont_plain_send.ydata[x];

      j++;
    }
  }
  else
  if(pf == 2)                           /* Take them from the beginning middle and end */
  {
    for(x=0;x<abs(nb);x++)
    {
      xd[j] = cont_plain_send.xdata[x];
      yd[j] = cont_plain_send.ydata[x];

      j++;
    }

    start = (n / 2) - (abs(nb) / 2);
    end   = (n / 2) + (abs(nb) / 2);

    for(x=start;x<end;x++)
    {
      xd[j] = cont_plain_send.xdata[x];
      yd[j] = cont_plain_send.ydata[x];

      j++;
    }

    for(x=n-abs(nb);x<n;x++)
    {
      xd[j] = cont_plain_send.xdata[x];
      yd[j] = cont_plain_send.ydata[x];

      j++;
    }
  }
  else                                  /* Take them from the beginning and end */
  {
    for(x=0;x<nb;x++)
    {
      xd[j] = cont_plain_send.xdata[x];
      yd[j] = cont_plain_send.ydata[x];

      j++;
    }

    for(x=n-nb;x<n;x++)
    {
      xd[j] = cont_plain_send.xdata[x];
      yd[j] = cont_plain_send.ydata[x];

      j++;
    }
  }

  gridmap_lsqfit(xd, yd, 0, j-1, sig, 0, &a, &b, &siga, &sigb, &chi2, &q);

  if(a != 0.0 && b != 0.0) /* Remove a baseline where there is data */
  {
    for(x=0;x<n;x++)
    {
      cont_plain_send.ydata[x] = cont_plain_send.ydata[x] -(a + b * cont_plain_send.xdata[x]);
    }
  }

  if(xgraphic)
  {
    sendCont(SEND_PLAIN);
  }

  return(0);
}


int doNutator()
{
  char tbuf[256];
  double scanno;

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

    return(1);
  }

  /* Set these to '0' so x-axis is set properly */
  env->map_w = 0;
  env->map_h = 0;

  scanno = atof(args->secondArg);

  sprintf(tbuf, "get %.2f", scanno);
  parseCmd(tbuf);

  strcpy(tbuf, "show");
  parseCmd(tbuf);

  removeContMapBaselines(isSet(PF_B_FLAG), env->mapbaselines);

  if(isSet(PF_G_FLAG))
  {
    contMapFitGauss();
  }

  return(0);
}


int doGmapClear()
{
  int ret;

  if(isSet(PF_H_FLAG))
  {
    printH_FlagHelp(args->firstArg);

    return(0);
  }

  ret = isSet(PF_Y_FLAG);
  if(!ret)                      /* Ask */
  {
    ret = getYesNoQuit("Delete Internal Grid Map: y/N ?");
  }

  if(ret)
  {
    st_report("Deleting Internal Grid Map");

    bzero((char *)&beam_map, sizeof(beam_map));

    remScanHeaders();                   /* Zero out scan Headers */

    env->map_w = 0;
    env->map_h = 0;
  }
  else
  {
    st_report("Skipping Grid Map Deleting");
  }

  return(0);
}


/* Just a frontend to gmap functions */
/* CAUTION: Flags are overloaded here */
int doGmap()
{
  int    indx;
  struct PARSE_FLAGS *p;

  p = findFirstFlag();

  if(p)
  {
    indx = p->index;    /* Save value */

    shiftFlags();       /* Remove this first flag and decriment is 'set' count */

    switch(indx)
    {
      case PF_B_FLAG: strcpy(args->firstArg, "gmapmapbase"); mapBaseline();   break;            /* Remove Baseline from Gridded Map */
      case PF_C_FLAG: strcpy(args->firstArg, "gmapslice");   gridMapSlice();  break;            /* Cut a slice of map */
      case PF_D_FLAG: strcpy(args->firstArg, "gmapdive");    gridMapDive();   break;            /* Plot Grid Map with varing dB's */
      case PF_E_FLAG: strcpy(args->firstArg, "gmapclear");   doGmapClear();   break;            /* Clear out internal Gridmap data */
      case PF_F_FLAG: strcpy(args->firstArg, "gmapgauss");   mapGaussFit();   break;            /* Fit Gaussian to a Gridded Map Center*/
      case PF_G_FLAG: strcpy(args->firstArg, "gmapgrid");    gridMapGrid();   break;            /* Accumulate Selected CONTMAP Scans */
      case PF_L_FLAG: strcpy(args->firstArg, "gmaplimits");  gridMapLimits(); break;            /* Set / Reset Plotting Limits */
      case PF_P_FLAG: strcpy(args->firstArg, "gmapplot");    gridMapPlot();   break;            /* Plot 3D Gaussian Gridded Map */
      case PF_R_FLAG: strcpy(args->firstArg, "gmaprotate");  gmapRotate();    break;            /* Rotate Gridmap and Snap images */
      case PF_S_FLAG: strcpy(args->firstArg, "gmapstats");   gridMapStats();  break;            /* Display Stats of Grid Map */
      case PF_X_FLAG: strcpy(args->firstArg, "gmapstats");   gridMapGrab();   break;            /* Force Gridmap Window to be Active */

      default: printH_FlagHelp(args->firstArg); break;
    }

    return(0);
  }

  printH_FlagHelp(args->firstArg);

  return(0);
}