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

/* LINE Gaussian Routines */
//struct GAUSSIAN_FIT {
//        int    xmode;
//        double peak;
//        double center;
//        double hp;
//        double sum;
//};

struct GAUSSIAN_FIT gfit;


int printGauss(image)
struct LP_IMAGE_INFO *image;
{
  int i, n;

  bzero((char *)&sat.gauss, sizeof(sat.gauss));

  sprintf(sat.gauss.results[0].string, "Results Of Last Gaussian Fit");

  if(env->jfact == 0.0)
  {
    if(env->issmt)
    {
      env->jfact = MG_J_FACTOR;
    }
    else
    {
      env->jfact = KP_J_FACTOR;
    }
  }

  switch(gfit.xmode)
  {
    case XMODE_VEL:
      sprintf(sat.gauss.results[1].string, "Peak                  (K): %9.3f", gfit.peak);
      sprintf(sat.gauss.results[2].string, "Peak Velocity      (KM/s): %9.3f", gfit.center);
      sprintf(sat.gauss.results[3].string, "Peak HP Width      (KM/s): %9.3f", gfit.hp);
      sprintf(sat.gauss.results[4].string, "Sum Under Gauss       (K): %9.3f", gfit.sumK);
      sprintf(sat.gauss.results[5].string, "Sum Under Gauss (Jk km/s): %9.3f", gfit.sumKm * env->jfact);
      n = 6;
      break;

    case XMODE_CHAN:
      sprintf(sat.gauss.results[1].string, "Peak                  (K): %9.3f", gfit.peak);
      sprintf(sat.gauss.results[2].string, "Peak Channel             : %9.3f", gfit.center);
      sprintf(sat.gauss.results[3].string, "Peak HP Width  (in chans): %9.3f", gfit.hp);
      sprintf(sat.gauss.results[4].string, "Sum Under Gauss       (K): %9.3f", gfit.sumK);
      sprintf(sat.gauss.results[5].string, "Sum Under Gauss (Jk km/s): %9.3f", gfit.sumKm * env->jfact);
      n = 6;
      break;

    case XMODE_FREQ:
      sprintf(sat.gauss.results[1].string, "Peak                  (K): %9.3f", gfit.peak);
      sprintf(sat.gauss.results[2].string, "Peak Frequency      (MHz): %9.3f", gfit.center);
      sprintf(sat.gauss.results[3].string, "Rest Frequency      (MHz): %12.6f", (image->head.restfreq+gfit.center) / 1000.0);
      sprintf(sat.gauss.results[4].string, "Peak HP Width       (MHz): %9.3f", gfit.hp);
      sprintf(sat.gauss.results[5].string, "Sum Under Gauss       (K): %9.3f", gfit.sumK);
      sprintf(sat.gauss.results[6].string, "Sum Under Gauss (Jk km/s): %9.3f", gfit.sumKm * env->jfact);
      n = 7;
      break;
  }

  if(image->xaxisotfoff)
  {
    sprintf(sat.gauss.results[1].string, "Peak                  (K): %9.3f", gfit.peak);
    sprintf(sat.gauss.results[2].string, "Peak Offset     (arcsecs): %9.3f", gfit.center);
    sprintf(sat.gauss.results[3].string, "Peak HP Width   (arcsecs): %9.3f", gfit.hp);
    sprintf(sat.gauss.results[4].string, "Sum Under Gauss       (K): %9.3f", gfit.sumK);

    n = 6;
  }

  for(i=0;i<n;i++)
  {
    st_report(sat.gauss.results[i].string);
    sat.gauss.results[i].used = 1;
  }

  return(0);
}


int clearGauss(image, clear)
struct LP_IMAGE_INFO *image;
int clear;
{
  if(clear)
  {
    bzero((char *)&image->ginfo, sizeof(struct GAUSSIAN)); /* Clear the gaussian struct */
  }

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

  bzero((char *)&sat.gauss, sizeof(sat.gauss));

  sendSat();

  return(0);
}


int copyGoodChars(out, in)
char *out, *in;
{
  char *cp, *cp2;

  cp2 = out;
  for(cp=in;*cp;cp++)
  {
    if(isalnum(*cp))
    {
      *cp2 = *cp;

       cp2++;
    }
  }

  *cp2 = '\0';

  return(0);
}



int findLineName(f, d, buf)
double f, d;
char *buf;
{
  int    try, j;
  char   tryBuf[256], line[256][256], *string;
  double fr;
  FILE *fp;

  fr = d;

  for(try=0;try<100;try++,fr *= 1.5)
  {
    sprintf(tryBuf, "slovas %f %f", f, fr);

    if(isSet(PF_V_FLAG))
    {
      st_print("TRY(%d): %s\n", try+1, tryBuf);
    }

    if((fp = popen(tryBuf, "r")) <= (FILE *)0)
    {
      st_fail("Unable to open: %s", tryBuf);
      return(1);
    }

    bzero(line, sizeof(line));

    j = 0;
    while((string = fgets(line[j], sizeof(line[j]), fp)) != NULL)
    {
      strcmprs(line[j]);

      if(isSet(PF_V_FLAG))
      {
        st_print("READ(%d) GOT: %s\n", j+1, line[j]);
      }

      j++;
    }

    pclose(fp);

    if(j)
    {
      if(j == 1)
      {
        copyGoodChars(buf, line[j-1]);
        st_report("findLineName(): Try: %d, Got: %s", try+1, buf);
      }

      return(0);
    }     
  }

  return(0);
}

int gaussList(header)
int header;
{
  int    i;
  struct LP_IMAGE_INFO *p;

  p = h0;
  
  if(p->ginfo.nfits)
  {
    if(header)
    {
      st_print( " Fit#   Peak          Offset           Rest-Freq        FWHM       Height-Error  Center-Error   FWHM-Error  Line-Name (Guess)\n\n");
    }

    for(i=0;i<p->ginfo.nfits;i++)
    {
      st_print("%2d   %9.6f     %12.6f      %12.8f    %9.6f      %9.6f      %9.6f    %9.6f    %s\n", 
                         i, 
                         p->ginfo.fits[i].peak, 
                         p->ginfo.fits[i].center, 
                        (p->ginfo.fits[i].center+p->head.restfreq) / 1000.0, 
                         p->ginfo.fits[i].hp, 
                         p->ginfo.fits[i].hghterr, 
                         p->ginfo.fits[i].cnterr, 
                         p->ginfo.fits[i].hwerr,
                         p->ginfo.fits[i].linename);
    }

    if(header)
    {
      st_print("\n");
    }
  }

  return(0);
}


int gaussRecord()
{
  char fname[256];
  FILE *fp;
  
  strcpy(fname, "gaussian.results");

  fp = fopen(fname, "a");

  if(!fp)
  {
    st_fail("gaussRecord();Unable to open %s for appending", fname);
    return(1);
  }

  fprintf(fp, "%8.2f - %8.2f  %10.6f  %9.3f  %9.3f  %9.3f  %9.3f  %9.6f  %9.6f  %9.6f\n", 
                         h0->ainfo.bscan, 
                         h0->ainfo.escan,
                        (env->gaussCenter + h0->head.restfreq) / 1000.0, 
                         env->gaussPeak, 
                         env->gaussCenter, 
                         env->gaussHP, 
                         env->gaussSum, 
                         env->hghterr, 
                         env->cnterr, 
                         env->hwerr);

  fclose(fp);

  st_report("gaussRecord(): Recorded Gaussian fit to %s", fname);

  return(0);
}


int autoNGauss(n)
int n;
{
  int    i, j, ret, len;
  int    pgauss = 0, h1gauss = 0, h2gauss = 0;
  double aa[4], y;
  double xd[MAX_CHANS], yd[MAX_CHANS], sig[MAX_CHANS];
  char   cmdBuf[256];
  float  ydata[MAX_CHANS];
  struct LP_IMAGE_INFO *image = h0;

  st_report("Auto Computing %d Gaussian Curves", n);

  /* First save a copy the plot data to a local array
     because I'm going to destroy the yplot1 array in the 
     process of fitting all the peaks; 
     copy it back when down */

  bcopy((char *)&image->yplot1[0], (char *)&ydata[0], sizeof(image->yplot1));
  len = image->head.noint;

  for(i=0;i<n;i++)
  {
    line_determine_min_max(image);

    st_report("Peak %d: Chan: %4d, Temp: %f", i, image->peakChan, image->datamax);

    pgauss = env->peakChan;

        /* Find the channel to the left of the line that's < 1/2 peak value */
    for(j=0;j<100;j++)
    {
      if(image->yplot1[pgauss-j] < (image->yplot1[env->peakChan] / 2.0)) /* Is it 1/2 the power? */
      {
        break;
      }
    }

    if(j > 99)
    {
      st_fail("Autofind Gauss fit Lower HP; failed");
      bcopy((char *)&ydata[0], (char *)&image->yplot1[0], sizeof(image->yplot1));
      line_determine_min_max(image);

      return(1);
    }

    h1gauss = pgauss-j;

    if(h1gauss < 0)
    {
      h1gauss = 0;
    }
    if(h2gauss >= len)
    {
      h1gauss = len-1;
    }

        /* Find the channel to the right of the line that's < 1/2 peak value */
    for(j=0;j<100;j++)
    {
      if(image->yplot1[pgauss+j] < (image->yplot1[env->peakChan] / 2.0)) /* Is it 1/2 the power? */
      {
        break;
      }
    }
    if(j > 99)
    {
      st_fail("Autofind Gauss fit Upper HP; failed");
      bcopy((char *)&ydata[0], (char *)&image->yplot1[0], sizeof(image->yplot1));
      line_determine_min_max(image);

      return(1);
    }
    h2gauss = pgauss+j;

    st_report("With a initial guess of: Peak %f @ Chan %d", image->yplot1[pgauss], pgauss);
    st_report("And Half-Power Points of: %d & %d", h1gauss, h2gauss);

    clearGauss(image, 0); /* clear out the local gfit */

    /* Copy data into double arrays for fitting */
    for(j=0;j<len;j++)
    {
      xd[j]  = image->x1data[j];
      yd[j]  = image->yplot1[j];
      sig[j] = 1.0;
    }

    aa[0] = image->yplot1[pgauss];
    aa[1] = image->x1data[pgauss];
    aa[2] = fabs(image->x1data[h1gauss]-image->x1data[h2gauss]);

    st_report("Gauss Guess:   %f %f %f", aa[0], aa[1], aa[2]);

    /* Use older Unipops version of gauss_fit */
    if(!env->usegslgauss)
    {
      st_report("Using Unipops Gaussian Routines");
      /* shift up one element */
      aa[3] = aa[2];
      aa[2] = aa[1];
      aa[1] = aa[0];
      aa[0] = 0.0;

      ret = gauss_fit(&xd[0], &yd[0], &sig[0], len-1, &aa[0], isSet(PF_V_FLAG));

      if(ret) /* return value is oppisite of gsl version */
      {
        ret = 0;
      }
      else
      {
        ret = 1;
      }

      aa[0] = aa[1];
      aa[1] = aa[2];
      aa[2] = aa[3];
      aa[3] = 0.0;
    }
    else
    {
      st_report("Using GSL Gaussian Routines");
      ret = gaussian_fit_gsl(&xd[0], &yd[0], len, &aa[0]);
      aa[2] *= 2.0; /* gsl functions return 1/2 the proper value; double it here */
    }

    if(!ret)
    {
      st_report("Gauss Results: %f %f %f", aa[0], aa[1], aa[2]);

      gfit.xmode   = env->x1mode;
      gfit.peak    = aa[0];
      gfit.center  = aa[1];
      gfit.hp      = aa[2];
      gfit.sumK    = 0.0;
      gfit.hghterr = hghterr;
      gfit.cnterr  = cnterr;
      gfit.hwerr   = hwerr;

      for(j=0;j<len;j++)
      {
        y = gaussian(aa[0], aa[1], aa[2], xd[j]);
        gfit.sumK += y;
      }

      env->gaussPeak   = gfit.peak;
      env->gaussCenter = gfit.center;
      env->gaussHP     = gfit.hp;
      env->hghterr     = gfit.hghterr;
      env->cnterr      = gfit.cnterr;
      env->hwerr       = gfit.hwerr;

      findLineName((gfit.center+image->head.restfreq) / 1000.0, fabs(image->head.freqres) / 1000.0, gfit.linename);

      st_report("%4d Channels @ deltax = %f", len, image->head.deltax);

      gfit.sumKm    = gfit.sumK * fabs(image->head.deltax);
      env->gaussSum = gfit.sumK;

      printGauss(image);

      bcopy((char *)&gfit, (char *)&image->ginfo.fits[image->ginfo.nfits], sizeof(struct GAUSSIAN_FIT));
      image->ginfo.nfits++;

      xx();

      if(isSet(PF_S_FLAG))
      {
        sprintf(cmdBuf, "splat /l");
        parseCmd(cmdBuf);
      }
    }
    else
    {
      st_fail("Gaussian Fit Failed");
      bcopy((char *)&ydata[0], (char *)&image->yplot1[0], sizeof(image->yplot1));
      line_determine_min_max(image);

      return(1);
    }

    /* Zero out these channels so next round doesn't see them */
    h1gauss -= 5;
    h2gauss += 5;

    for(j=h1gauss;j<h2gauss;j++)
    {
      image->yplot1[j] = 0.0;
    }

  }

  bcopy((char *)&ydata[0], (char *)&image->yplot1[0], sizeof(image->yplot1));
  line_determine_min_max(image);

  gaussList(1);

  return(0);
}


int do_line_gauss()
{
  int    i, k, ret, len, n, rn;
  int    pgauss = 0, h1gauss = 0, h2gauss = 0;
  double aa[4], y;
  double xd[MAX_CHANS], yd[MAX_CHANS], sig[MAX_CHANS];
  char   cmdBuf[256];
  struct LP_IMAGE_INFO *image = h0;

  if(isSet(PF_C_FLAG))
  {
    st_report("Clearing out all previous gaussian fits");

    clearGauss(image, 1);


    return(0);
  }

  if(isSet(PF_L_FLAG))
  {
    gaussList(1);

    return(0);
  }

  xx();                          /* Repaint with any previous gaussian curves cleared */

  st_report("Gaussian Fitting Routine");

  rn = getNset();

  if(isSet(PF_A_FLAG) && rn)
  {
    return(autoNGauss(rn));
  }
  else
  if(isSet(PF_A_FLAG) && !rn)
  {
    return(autoNGauss(1));
  }
  else
  {
    do
    {
      n = getNumericResponse("Enter Number of Gaussians to fit [ Up to %d]: ", MAX_GAUSSIAN_FITS);
    }
    while(n < 0);

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

    if(n > MAX_GAUSSIAN_FITS)
    {
      st_fail("Too many fit requests: max = %d", MAX_GAUSSIAN_FITS);
      return(1);
    }

    st_report("You will need to make 3 Mouse Clicks per Gaussian");
  }

  for(k=0;k<n;k++)
  {
    if(n > 1)
    {
      st_report(" ");
      st_report("Fit #%d", k);
      st_report(" ");
    }

    ret = _doCCur("Mouse-Click plot on PEAK of Gaussian Area");
    if(ret == 1)
    {
      st_fail("Gauss fit Cancelled");
      return(1);
    }

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

    if(pgauss == 0)
    {
      st_fail("Gauss fit Cancelled");
      return(1);
    }

    st_report("gauss center channel: %d, temp = %f", pgauss, image->yplot1[pgauss]);

    ret = _doCCur("Mouse-Click plot on LEFT edge of Gaussian Half-Power Point");

    if(ret == 1)
    {
      st_fail("Gauss fit Cancelled");
      return(1);
    }

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

    st_report("gauss left channel: %d, temp = %f", h1gauss, image->yplot1[h1gauss]);

    if(h1gauss == 0)
    {
      st_fail("Gauss fit Cancelled");
      return(1);
    }

    ret = _doCCur("Mouse-Click plot on RIGHT edge of Gaussian Half-Power Point");

    if(ret == 1)
    {
      st_fail("Gauss fit Cancelled");
      return(1);
    }

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

    st_report("gauss right channel: %d, temp = %f", h2gauss, image->yplot1[h2gauss]);

    if(h2gauss == 0)
    {
      st_fail("Gauss fit Cancelled");
      return(1);
    }

    st_report("With a initial guess of: Peak %f @ Chan %d", image->yplot1[pgauss], pgauss);
    st_report("And Half-Power Points of: %d & %d", h1gauss, h2gauss);

    clearGauss(image, 0); /* clear out the local gfit */

    len = image->head.noint;

    for(i=0;i<len;i++)
    {
      xd[i]  = image->x1data[i];
      yd[i]  = image->yplot1[i];
      sig[i] = 1.0;
    }

    aa[0] = yd[pgauss];
    aa[1] = xd[pgauss];
    aa[2] = fabs(xd[h1gauss]-xd[h2gauss]);

    st_report("Gauss Guess:   %f %f %f", aa[0], aa[1], aa[2]);

  /* computer gaussian using Temp / Channel */
    /* Use older Unipops version of gauss_fit */
    if(!env->usegslgauss)
    {
      st_report("Using Unipops Gaussian Routines");
      /* shift up one element for Fortran */
      aa[3] = aa[2];
      aa[2] = aa[1];
      aa[1] = aa[0];
      aa[0] = 0.0;

      ret = gauss_fit(&xd[0], &yd[0], &sig[0], len-1, &aa[0], isSet(PF_V_FLAG));

      if(ret) /* return value is oppisite of gsl version */
      {
        ret = 0;
      }
      else
      {
        ret = 1;
      }

      /* shift back down one element for Unix */
      aa[0] = aa[1];
      aa[1] = aa[2];
      aa[2] = aa[3];
      aa[3] = 0.0;
    }
    else
    {
      st_report("Using GSL Gaussian Routines");
      ret = gaussian_fit_gsl(&xd[0], &yd[0], len, &aa[0]);
      aa[2] *= 2.0; /* gsl functions return 1/2 the proper value; double it here */
    }

    if(!ret)
    {
      st_report("Gauss Results: %f %f %f", aa[0], aa[1], aa[2]);

      gfit.xmode   = env->x1mode;
      gfit.peak    = aa[0];
      gfit.center  = aa[1];
      gfit.hp      = aa[2];
      gfit.sumK    = 0.0;
      gfit.hghterr = hghterr;
      gfit.cnterr  = cnterr;
      gfit.hwerr   = hwerr;

      for(i=0;i<len;i++)
      {
        y = gaussian(aa[0], aa[1], aa[2], xd[i]);
        gfit.sumK += y;
      }

      env->gaussPeak   = gfit.peak;
      env->gaussCenter = gfit.center;
      env->gaussHP     = gfit.hp;
      env->hghterr     = gfit.hghterr;
      env->cnterr      = gfit.cnterr;
      env->hwerr       = gfit.hwerr;

      findLineName((gfit.center+image->head.restfreq) / 1000.0, fabs(image->head.freqres) / 1000.0, gfit.linename);

      st_report("%4d Channels @ deltax = %f", len, image->head.deltax);

      gfit.sumKm    = gfit.sumK * fabs(image->head.deltax);
      env->gaussSum = gfit.sumK;

      printGauss(image);

      bcopy((char *)&gfit, (char *)&image->ginfo.fits[image->ginfo.nfits], sizeof(struct GAUSSIAN_FIT));
      image->ginfo.nfits++;

      xx();

      if(isSet(PF_S_FLAG))
      {
        sprintf(cmdBuf, "splat /l");
        parseCmd(cmdBuf);
      }
    }
    else
    {
      st_fail("Gaussian Fit Failed");

      return(1);
    }

  }

  gaussList(1);

  return(0);
}


/* Continuum Gaussian Routines */


int cont_gauss(y1, y2, y3, x1, x2, x3, a, b, c)
float  y1, y2, y3;                                        /* Y coordinates   */
double x1, x2, x3;                                         /* X coordinates   */
float  *a, *b, *c;                                         /* curve constants */
{
  float l2, l3, t;

  if( (t = y2/y1) > 0.0)
  {
    l2 = log(t);
  }
  else
  {
    return(0);
  }

  if( (t = y3/y1) > 0.0)
  {
    l3 = log(t);
  }
  else
  {
    return(0);
  }

  *c = (-l2 * x1 * x1 + l3 * x1 * x1 - l3 * x2 * x2 + l2 * x3 * x3)/
       (2 * l2 * x1 - 2 * l3 * x1 + 2 * l3 * x2 - 2 * l2 * x3);

  *b = (l2*x1 - l3*x1 + l3*x2 - l2*x3)/
       (x1*x1*x2 - x1*x2*x2 - x1*x1*x3 + x2*x2*x3 + x1*x3*x3 - x2*x3*x3);

  *a =  y1 *exp(*b * (*c + x1) * (*c + x1));

  if(*b < 0.0)
  {
    return(0);
  }

  return(1);
}


int doGauss()
{
  int  i, n;
  char buf[256], bbuf[512];

  if(h0->head.scan < 1.0)
  {
    st_warn("No Scan loaded; Use 'get scan#'");
    return(0);
  }

//  if(isSet(PF_R_FLAG)) /* Record the results ?? */
//  {
//    gaussRecord(0);
//
//    return(0);
//  }

  if(isSet(PF_N_FLAG))
  {
    if(args->ntokens > 2)
    {
      n = atoi(args->secondArg);

      if(n < 0 || n >= MAX_GAUSSIAN_FITS)
      {
        st_print("gauss /n: arg out of range: 0 to %d\n", MAX_GAUSSIAN_FITS-1);

        printH_FlagHelp(args->firstArg);
        return(1);
      }

      bzero(bbuf, sizeof(bbuf));

      for(i=2;i<args->ntokens;i++)
      {
        sprintf(buf, args->tok[i]);

        strcat(bbuf, buf);
      }

      st_report("Setting Gaussian fit[%d] name to %s", n, bbuf);

      strcpy(h0->ginfo.fits[n].linename, bbuf);

      return(0);
    }
    else
    {
      st_print("Missing arg for gauss /n\n");

      printH_FlagHelp(args->firstArg);

      return(1);
    }
  }

  if(env->mode == MODE_CONT)
  {
    do_cont_gauss(h0->yplot1, &h0->head);
  }
  else /* Works for both OTF and SPEC */
  {
    do_line_gauss();
    xx();
  }

  return(0);
}


int do_cont_gauss(d, h)
float *d;
struct HEADER *h;
{
  float  ydata[MAX_CONT_SAMPLES];
  int    i, dlen, gstart;
  double aa1[4];
  double xd[256], yd[256], s[256], deltax;
  char sysBuf[256];

  st_report("do_cont_gauss(): Enter");

  if(processContinuum(h, d, &ydata[0], &dlen, "do_cont_gauss"))
  {
    return(1);
  }

  deltax = findDeltaX(h);

  /* Copy into (double)gauss array */
  for(i=0;i<dlen;i++)
  {
    xd[i] = (double)i * deltax - ((double)dlen * deltax / 2.0) + deltax / 2.0;
    yd[i] = ydata[i];
    s[i]  = 1.0;
  }

  aa1[0] = 0.0;
  aa1[1] = env->datamax;
  aa1[2] = xd[env->peakChan];
  aa1[3] = h->bfwhm;

  st_report("do_cont_gauss(): Guess: Center = %f @ %f: width = %f", aa1[1], aa1[2], aa1[3]);

  if(env->glength <= 0)
  {
    env->glength = 15;
  }
  else
  if(env->glength > dlen)
  {
    env->glength = dlen - 2;
  }

  gstart = dlen / 2 - 5;

  gstart = env->peakChan - env->glength / 2;

  gauss_fit(&xd[gstart], &yd[gstart], &s[gstart], env->glength, &aa1[0], isSet(PF_V_FLAG));

  st_report("do_cont_gauss(): %d-Point Gaussian Fit: Y = %.2f, X = %.2f, BW = %.2f",
                                                env->glength, aa1[1], aa1[2], aa1[3]);
  beam_map.gmap_peak1 = aa1[1];
  beam_map.gmap_bw1   = aa1[3];

  beam_map.gmap_peak2 = 0.0;
  beam_map.gmap_bw2   = 0.0;

  doDeConvol(); /* Deconvol out the planet */

  if(isSet(PF_P_FLAG))
  {
    st_report("do_cont_gauss(): Displaying Plot");
    plotGauss(&xd[0], dlen, &aa1[0], h, 0);
  }

  if(isSet(PF_S_FLAG))
  {
    setFlag(PF_C_FLAG, 1); /* Set for capture */

    printPlot();

    usleep(1000000);

    st_report("do_cont_gauss(): Saving Plot");
    sprintf(sysBuf, "cp /tmp/%s/ccal_plot.ps Scan-%.2f.ps", session, h->scan);
    system(sysBuf);
  }

  return(0);
}


int plotGauss(x, len, a, h, mode)
double *x, *a;
int len, mode;
struct HEADER *h;
{
  int i;
  char  tbuf[512];
  double yd, dyda[4], max = -9999999.9;

  st_report("plotGauss(): Enter");

  if(cont_plain_send.dlen == 0)
  {
    st_fail("plotGauss(): Must load a scan and do a 'show' first");
    return(1);
  }
  /* Copy gaussian data to Struct*/

  for(i=0;i<len;i++)
  {
    fgauss(x[i], a, &yd, &dyda[0]);

    cont_plain_send.gdata[i] = yd;
  }

  if(env->yaxis_scale == LOG_PLOT)
  {
    for(i=0;i<len;i++)
    {
      if(cont_plain_send.gdata[i] > max)
      {
        max = cont_plain_send.gdata[i];
      }
    }

    if(max != 0.0)
    {
      for(i=0;i<len;i++)
      {
        cont_plain_send.gdata[i] = 10.0 * log10(cont_plain_send.gdata[i] / max);
      }
    }
  }


  cont_plain_send.ngauss = len;

  sprintf(tbuf, "%d Point Gaussian Fit",      env->glength);
  strcpy(cont_plain_send.glabels[0], tbuf);

  sprintf(tbuf, "Gaussian Peak:  %10.2f K",   a[1]);
  strcpy(cont_plain_send.glabels[1], tbuf);

  if(mode == 2)
  {
    sprintf(tbuf, "Gaussian     BW:    %6.2f mm", a[3]);
  }
  else
  {
    sprintf(tbuf, "Gaussian     BW:    %6.2f sec", a[3]);
  }

  strcpy(cont_plain_send.glabels[2], tbuf);

  if(mode != 2)
  {
    sprintf(tbuf, "De-Convolved BW:    %6.2f sec", decon_beam_size);
    strcpy(cont_plain_send.glabels[3], tbuf);
  }

  if(mode == 2)
  {
    sprintf(tbuf, "AX Focus Input: %6.2f mm, Results:  %6.2f mm",  h->focusr, a[2]);
  }
  else
  if(mode == 1)
  {
    sprintf(tbuf, "AX Focus Input: %6.2f mm, Results:  %6.2f mm",  h->focusr, a[2]);
  }
  else
  {
    sprintf(tbuf, "AX Focus Setting:   %6.2f mm",  h->focusr);
  }
  strcpy(cont_plain_send.glabels[4], tbuf);

  sprintf(tbuf, "NS Focus Setting:   %6.2f mm",  h->focusv);
  strcpy(cont_plain_send.glabels[5], tbuf);

  sprintf(tbuf, "EW Focus Setting:   %6.2f mm",  h->focush);
  strcpy(cont_plain_send.glabels[6], tbuf);

  if(mode != 2)
  {
    sprintf(tbuf, "Map Angle:          %6d Degs",  env->mapangle);
    strcpy(cont_plain_send.glabels[7], tbuf);
  }

  /* Save latest in global */
  env->gaussPeak   = a[1];
  env->gaussCenter = a[2];
  env->gaussHP     = a[3];

  if(xgraphic && notSet(PF_X_FLAG))
  {
    sendCont(SEND_PLAIN);
  }

  return(0);
}