focus.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 FOCUS_HISTORY_SEND focus_send;

int focus_axis = FOCUS_AXIS_AX;

int sendFocus()
{
  if(!focus_send.nfocs)
  {
    st_fail("No Focus Scans Loaded yet");
    return(1);
  }

  strcpy(focus_send.magic, "FOCUS_SEND");

  focus_send.plotwin = env->plotwin;
  focus_send.axis    = focus_axis;;

  if(last_mouse_x != 0.0 && last_mouse_y != 0.0)
  {
    focus_send.xlabel  = last_mouse_x;
    focus_send.ylabel  = last_mouse_y;
  }

  sock_write(xgraphic, (char *)&focus_send, sizeof(focus_send));

  st_report("Focus Struct send to Xgraphic");

  return(0);
}

int focusSortEl(p1, p2)
struct FOCUS_HISTORY_ELEM *p1, *p2;
{
  if(p1->el < p2->el)
  {
    return(-1);
  }
  else
  if(p1->el > p2->el)
  {
    return(1);
  }
  else
  {
    return(0);
  }
}


int focusSortTemp(p1, p2)
struct FOCUS_HISTORY_ELEM *p1, *p2;
{
  if(p1->tamb < p2->tamb)
  {
    return(-1);
  }
  else
  if(p1->tamb > p2->tamb)
  {
    return(1);
  }
  else
  {
    return(0);
  }
}


int focusSortScan(p1, p2)
struct FOCUS_HISTORY_ELEM *p1, *p2;
{
  if(p1->scan < p2->scan)
  {
    return(-1);
  }
  else
  if(p1->scan > p2->scan)
  {
    return(1);
  }
  else
  {
    return(0);
  }
}


int dofocusSort()
{

  if(isSet(PF_I_FLAG))
  {
    st_report("dofocusSort(): Sorting %d records by Scan Number", focus_send.nfocs);
    qsort(&focus_send.focuses, focus_send.nfocs, sizeof(struct FOCUS_HISTORY_ELEM), (__compar_fn_t)focusSortScan);
    focus_send.sort_order = FOCUS_SORTED_SCAN;
  }
  else
  if(isSet(PF_T_FLAG))
  {
    st_report("dofocusSort(): Sorting %d records by Ambient Temperature", focus_send.nfocs);
    qsort(&focus_send.focuses, focus_send.nfocs, sizeof(struct FOCUS_HISTORY_ELEM), (__compar_fn_t)focusSortTemp);
    focus_send.sort_order = FOCUS_SORTED_TAMB;
  }
  else
  {
    st_report("dofocusSort(): Sorting %d records by Elevation", focus_send.nfocs);
    qsort(&focus_send.focuses, focus_send.nfocs, sizeof(struct FOCUS_HISTORY_ELEM), (__compar_fn_t)focusSortEl);
    focus_send.sort_order = FOCUS_SORTED_ELEV;
  }

  return(0);
}



int findLoadFocuses()
{
  int    i, j, ret, good, oc;
  double F0;
  struct HEADER *p;
  struct FOCUS_HISTORY_ELEM *e;

  resetRejections();
  env->last_scan_count = 0;

  bzero((char *)&focus_send, sizeof(focus_send));
  strcpy(focus_send.magic, "FOCUS_SEND");

  focus_send.lines = MODE_LINES;

  j = 0;
  e = &focus_send.focuses[j];
  p = &sa[j];

  clearInterrupt();

  for(i=j;i<env->numofscans;i++,p++)
  {
    if(isInterrupt())
    {
      env->last_scan_count = j;
      focus_send.nfocs     = j;

      if(j > 5)
      {
        dofocusSort();
      }

      return(1);
    }

    good = 0;

    if(!meetsCritiria(p, CHECK_ALL))
    {
      continue;
    }

    if(isSet(PF_N_FLAG))        /* Load NSFC insteads */
    {
      focus_axis = FOCUS_AXIS_NS;

      if(strncmp(p->obsmode, "CONTNSFC", 8))
      {
        incrementRejMode();
        continue;
      }

      F0 = p->focusv;
    }
    else
    if(isSet(PF_E_FLAG))        /* Load EWFC insteads */
    {
      focus_axis = FOCUS_AXIS_EW;

      if(strncmp(p->obsmode, "CONTEWFC", 8))
      {
        incrementRejMode();
        continue;
      }

      F0 = p->focusv;
    }
    else
    {
      focus_axis = FOCUS_AXIS_AX;

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

      F0 = p->focusr;
    }

    if(isScanIgnored(p->scan))
    {
      oc = oldColor;
      setColor(YELLOW);
      st_report("Ignoring scan %.2f", p->scan);
      setColor(oc);

      incrementRejIgnored();
      continue;
    }

    /* Check if tsys is bad */
    if((env->tsyslimit && (fabs(p->stsys) > env->tsyslimit)) || gsl_isnan(p->stsys))
    {
      oc = oldColor;
      setColor(YELLOW);
      st_report("Rejecting scan %.2f, Tsys Greater then %.1f: %.1f", p->scan, env->tsyslimit, p->stsys);
      setColor(oc);

      incrementBogusTsys();

      continue;
    }

    if(env->focus_freq_upper != 0.0 && env->focus_freq_lower != 0.0)
    {
      if((p->restfreq / 1e3) < env->focus_freq_lower || (p->restfreq / 1e3) > env->focus_freq_upper)
      {
        oc = oldColor;
        setColor(YELLOW);
        st_report("Rejecting scan %.2f, Out of Band", p->scan);
        setColor(oc);

        incrementRejFreq();

        continue;
      }
    }

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

    if(ret > 0.0)
    {
      showScan();

      strncpy(e->source, p->object, 16);
      strcmprs(e->source);

        /* Only save if they meet limits */
      if(env->gaussPeak > env->focus_lower_limit)
      {
        good++;
      }
      else
      {
        oc = oldColor;
        setColor(YELLOW);
        st_report("Rejecting scan %.2f for being too weak", p->scan);
        setColor(oc);

        incrementBogusData();

        continue;
      }

      if(env->focus_cnterr_limit)
      { 
        if(cnterr < env->focus_cnterr_limit)
        {
          good++;
        }
        else
        {
          oc = oldColor;
          setColor(YELLOW);
          st_report("Rejecting scan %.2f for being too noisy", p->scan);
          setColor(oc);

          incrementBogusData();

          continue;
        }
      }
      else
      {
        good++;
      }

      if(env->focus_center_diff)
      { 
        if(fabs(F0-env->gaussCenter) < env->focus_center_diff)
        {
          good++;
        }
        else
        {
          oc = oldColor;
          setColor(YELLOW);
          st_report("Rejecting scan %.2f for Fit being too diff", p->scan);
          setColor(oc);

          incrementBogusData();

          continue;
        }
      }
      else
      {
        good++;
      }

      if(good == 3)
      {
        e->scan   = (float)p->scan;
        e->freq   = (float)p->restfreq / 1000.0;
        e->el     = (float)p->el;

        if(isSet(PF_N_FLAG))    /* Load NSFC insteads */
        {
          e->F0     = (float)p->focusv;
        }
        else
        {
          e->F0     = (float)p->focusr;
        }
        e->fit    = (float)env->gaussCenter;
        e->peak   = (float)env->gaussPeak;
        e->hp     = (float)env->gaussHP;
        e->cnterr = (float)cnterr;
        e->tamb   = (float)p->tamb;

        j++;

        if(j >= MAX_FOCUS_HISTORY)
        {
          st_report("Max Num of Focuses, %d,  Loaded", MAX_FOCUS_HISTORY);

          break;
        }

        e++;
      }

      if(notSet(PF_X_FLAG))             /* if /x flag set don't show the plot; goes faster on remote screens */
      {
        if(env->xdelay)
        {
          usleep(env->xdelay);
        }
        else
        {
          usleep(100000);
        }
      }
    }
  }

  env->last_scan_count = j;

  focus_send.nfocs   = j;

  dofocusSort();

  st_report("findLoadFocuses(): Wrote %d records to: focus_send", j);

  sprintf(focus_send.blabels[0], "Focus Results: %d Scans", focus_send.nfocs);
  sprintf(focus_send.blabels[1], "Scans: %.2f to %.2f", 
                        focus_send.focuses[0].scan, 
                        focus_send.focuses[focus_send.nfocs-1].scan);

  return(0);
}



int doFocusGracefit()
{
  int    i, k, n;
  float  *f;
  char   tbuf[256], cmdBuf[512];
  char   formulaBuf[256], form[256];
  double xd[MAX_FOCUS_HISTORY], yd[MAX_FOCUS_HISTORY], coeff[32], avg;
  struct FOCUS_HISTORY_SEND *foc;
  struct FOCUS_HISTORY_ELEM *e;

  st_report("doFocusGracefit(): Fitting a %d order LSQ to Loaded Focus Scans", env->nfit);

  if(focus_send.nfocs < 5)
  {
    st_fail("Little or No Focus Scans Loaded");
    return(1);
  }

  bzero((char *)xd,    sizeof(xd));
  bzero((char *)yd,    sizeof(yd));
  bzero((char *)coeff, sizeof(coeff));


  if(env->nfit == 0)
  {
    foc = &focus_send;
    e   = &foc->focuses[0];

    n = foc->nfocs;

    avg  = 0.0;

    for(i=0;i<n;i++,e++)
    {
      avg += e->fit;
    }

    avg /= (double)n;

    st_report("doFocusGracefit(): AVERAGE: %f", avg);

    f = &foc->bplot[0];
    for(i=0;i<n;i++,f++)
    {
      *f = avg;
    }

    sprintf(focus_send.blabels[2], "Focus Fit: Average: %f", avg);

    foc->bfit = n;
  }
  else
  if(env->nfit > 0 && env->nfit < MAX_NFIT)
  {
    foc = &focus_send;
    e   = &foc->focuses[0];

    n = foc->nfocs;

    for(i=0;i<n;i++,e++)
    {
      if(foc->sort_order == FOCUS_SORTED_SCAN)
      {
        xd[i] = e->scan;
      }
      else
      if(foc->sort_order == FOCUS_SORTED_ELEV)
      {
        xd[i] = e->el;
      }
      else
      {
        xd[i] = e->tamb;
      }

      yd[i] = e->fit;
    }

    fitcurve(xd, yd, n, env->nfit, coeff);

    st_print("Linuxpops(C) >>  # ");

    for(k=0;k<env->nfit+1;k++)
    {
      st_print("a%d = %.4e ", k, coeff[k]);
    }
    st_print("\n");

    /* now make a formula based on fit */

    if(env->nfit < 5)
    {
      strcpy(formulaBuf, "(");

      for(k=0;k<env->nfit+1;k++)
      {
        switch(k)
        {
          case 0: sprintf(form, "%e ", coeff[k]); break;
          case 1: sprintf(form, "+ (%e * el) ", coeff[k]); break;
          case 2: sprintf(form, "+ (%e * (el * el)) ", coeff[k]); break;
          case 3: sprintf(form, "+ (%e * (el * el * el)) ", coeff[k]); break;
          case 4: sprintf(form, "+ (%e * (el * el * el * el)) ", coeff[k]); break;
        }

        strcat(formulaBuf, form);
      }

      strcat(formulaBuf, ")");

      st_report("Formula: %s", formulaBuf);

      if(isSet(PF_A_FLAG)) /* auto flag; save to 'nfit' formula slot */
      {
        sprintf(cmdBuf, "formula /d /%d %s", env->nfit, formulaBuf);
        parseCmd(cmdBuf);
      }
    }

    f = &foc->bplot[0];
    e = &foc->focuses[0];

    for(i=0;i<n;i++,e++,f++)
    {
      if(foc->sort_order == FOCUS_SORTED_SCAN)
      {
        *f = (float)leasev(coeff, env->nfit, e->scan);
      }
      else
      if(foc->sort_order == FOCUS_SORTED_ELEV)
      {
        *f = (float)leasev(coeff, env->nfit, e->el);
      }
      else
      {
        *f = (float)leasev(coeff, env->nfit, e->tamb);
      }
    }

    if(env->nfit < 12) /* any more and it won't fit in the blabel */
    {
      sprintf(focus_send.blabels[2], "Focus Fit: ");

      for(k=0;k<env->nfit+1;k++)
      {
        sprintf(tbuf, "a%d = %.4e ", k, coeff[k]);

        strcat(focus_send.blabels[2], tbuf);
      }
    }
    else
    {
      strcpy(focus_send.blabels[2], "");
    }

    foc->bfit = n;
  }
  else
  {
    st_fail("doFocusGracefit(): Error: nfit out of range: %d, [0 - %d]", env->nfit, MAX_NFIT-1);
    return(1);
  }

  return(0);
}


int printFocusHeader()
{
  st_print("\n");
  st_print("Num        Source        Scan      Freq       Elev       Peak       F0       Fit      Diff    Center-Err  Tamb  Scaling\n");
  st_print("\n");

  return(0);
}


int listLoadedFocuses()
{
  int    i, n;
  char   tbuf[512];
  struct FOCUS_HISTORY_ELEM *p;

  n = focus_send.nfocs;

  if(!n)
  {
    st_fail("No Focus Scans Loaded yet");
    return(1);
  }


  p = &focus_send.focuses[0];

  for(i=0;i<n;i++,p++)
  {
    if(!(i%30))
    {
      printFocusHeader();
    }

    sprintf(tbuf, "%3d: %16.16s %9.2f %9.3f %9.2f %9.2f %9.3f %9.3f %9.3f %10.6f %7.1f %7.3f", 
                i, 
                p->source, 
                p->scan, 
                p->freq, 
                p->el, 
                p->peak, 
                p->F0, 
                p->fit, 
                p->F0-p->fit,
                p->cnterr,
                p->tamb,
                p->scaling);

    st_print("%s\n", tbuf);
  }

  return(0);
}


int listFileFocuses()
{
  int    i;
  struct HEADER *p;

  resetRejections();
  env->last_scan_count = 0;

  bzero((char *)&focus_send, sizeof(focus_send));
  strcpy(focus_send.magic, "FOCUS_SEND");

  st_print("%s\n\n", listHeader);

  p = &sa[0];
  for(i=0;i<env->numofscans;i++,p++)
  {                                                   /* If we want to auto reduce tips */
    if(!meetsCritiria(p, CHECK_ALL))
    {
      continue;
    }

    if(strncmp(p->obsmode, "CONTFOC", 7) && strncmp(p->obsmode, "CONTNSFC", 8) )
    {
      incrementRejMode();
      continue;
    }

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

    /* Check if tsys is bad */
    if((env->tsyslimit && (fabs(p->stsys) > env->tsyslimit)) || gsl_isnan(p->stsys))
    {
      printf("Scan %.2f Tsys Greater then %.1f: %.1f\n", p->scan, env->tsyslimit, p->stsys);
      incrementBogusTsys();

      continue;
    }

    printScan(p, 0);
    env->last_scan_count++;
  }

  st_report("Found %d Focus scans meeting criteria", env->last_scan_count);

  return(0);
}



int doFocusScaling()
{
  int    i, n;
  char   cmdBuf[256];
  struct FORMULA *f;
  struct FOCUS_HISTORY_ELEM *p;

  n = getNumericFlag();

  if(n<0)
  {
    st_fail("No numeric flag passed to 'focus /c'");
    return(0);
  }

  sprintf(cmdBuf, "double tmp");
  parseCmd(cmdBuf);

  f = &formulas.formulas[n];

  if(f->used)                   /* There is a formula here; create a scaling factor for each element */
  {
    n = focus_send.nfocs;

    if(!n)
    {
      st_fail("No Focus Scans Loaded yet");
      return(1);
    }

    p = &focus_send.focuses[0];

        /* Here comes the fun part; in order for the parse to properly
           substitute header variables in it's normal functions, we need
           to re-fetch each scan so the header is properly populated */
    for(i=0;i<n;i++,p++) 
    {
      sprintf(cmdBuf, "get %.2f /q", p->scan);
      parseCmd(cmdBuf);

      sprintf(cmdBuf, "tmp = %s /q", f->formula);
      parseCmd(cmdBuf);

      p->scaling = env->rst;
    }
  }
  else
  {
    st_report("Formula %d, not defined; removing correction", n);
    focus_send.scale = FOCUS_SCALE_NONE;
  }

  if(isSet(PF_ADD_FLAG))
  {
    focus_send.scale = FOCUS_SCALE_ADD;
  }
  else
  if(isSet(PF_SUB_FLAG))
  {
    focus_send.scale = FOCUS_SCALE_SUB;
  }
  else
  if(isSet(PF_MUL_FLAG))
  {
    focus_send.scale = FOCUS_SCALE_MULTI;
  }
  else
  if(isSet(PF_DIV_FLAG))
  {
    focus_send.scale = FOCUS_SCALE_DIVIDE;
  }
  else
  {
    st_report("No Function Flag; removing correction");
    focus_send.scale = FOCUS_SCALE_NONE;
  }

  sendFocus();

  return(0);
}



int doFocus()
{
  if(isSet(PF_R_FLAG))          /* Read them In and Analize */
  {
    st_report("doFocus(): Reading in Focus Scans");

    findLoadFocuses();

    sendFocus();

    return(0);
  }

  if(isSet(PF_C_FLAG))
  {
    return(doFocusScaling());
  }

  /* the '/n' is used in the '/r' function above, but there
     is a return so this flag does double duty */

  if(isSet(PF_N_FLAG))  /* Suppress Lines */
  {
    focus_send.lines = MODE_NOLINES;
  }
  else
  {
    focus_send.lines = MODE_LINES;
  }

  if(isSet(PF_F_FLAG))          /* Find and List */
  {
    st_report("doFocus(): Listing Available Focus Scans in file");

    listFileFocuses();

    return(0);
  }

  if(isSet(PF_L_FLAG))          /* List */
  {
    st_report("doFocus(): Listing Loaded Focus Scans");

    listLoadedFocuses();

    return(0);
  }

  if(isSet(PF_S_FLAG))          /* Sort */
  {
    dofocusSort();
    sendFocus();

    return(0);
  }

  if(isSet(PF_G_FLAG) == 1)             /* Fit a Guassian */
  {
    doFocusGracefit();
    sendFocus();

    return(0);
  }

  if(isSet(PF_G_FLAG) == 2)             /* Erase any Guassian */
  {
    focus_send.bfit = 0;
    bzero((char *)&focus_send.blabels[0], sizeof(focus_send.blabels));

    sendFocus();

    return(0);
  }

  if(isSet(PF_D_FLAG))          /* Re-Display Results */
  {
    sendFocus();
  }

  if(isSet(PF_M_FLAG))
  {
    sendFocus();                    /* Force a repaint first */

    getLabelPos("Click Mouse to Place Labels");

    focus_send.xlabel = last_mouse_x;
    focus_send.ylabel = last_mouse_y;

    sendFocus();

    return(0);
  }

  /* Default */
  sendFocus();

  return(0);
}