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

#define TESTING 1

#define MAX_RESULTS 200
#define SPL         0.299792458
#define MAX_DECON_SIZE 1024

#define MAX_GOOD_SOURCES 5
char *good_sources[] = { "Mercury", "Venus", "Mars", "Saturn", "Jupiter"};

double b_array[MAX_RESULTS];
double s_array[MAX_RESULTS];
double r_array[MAX_RESULTS];
double x_array[MAX_RESULTS];
double y_array[MAX_RESULTS];

double decon_beam_size;

struct RESULTS {
        int i;
        double beam;
        double source;
        double fit;
};


struct RESULTS results[MAX_RESULTS];


struct DECONVOL_STRUCT {
        char   planet[32];
        int    initialized;
        int    telescope;
        int    count;
        double source_size;
        double freq_array[MAX_DECON_SIZE];
        double fw_array[MAX_DECON_SIZE];
        double bs_array[MAX_DECON_SIZE];
};

struct DECONVOL_STRUCT dcs;

int mon;

char convol_tok[MAX_TOKENS][MAX_TOKENS];

int convol_get_tok(edit)
char *edit;
{
  int c = 0;
  char *token;

  token = strtok(edit, " \n\t");
  while(token != NULL)
  {
    sprintf(convol_tok[c++], "%s", token);
    token = strtok(NULL, " \n\t");
  }

  return(c);
}

double fwhm(freq)
double freq;
{
  double lambda, theta, hm, dish_size, taper;

  if(env->issmt)
  {
    dish_size = 10.0;
    taper     = 1.28;
  }
  else
  {
    dish_size = 12.0;
    taper     = 1.12; /* 10 dB edge taper; derived Jan-05-2016; glauria/tfolkers */
                      /* Was 1.22 then changed to 1.28; now 1.15 */
//    taper     = 1.28;
  }

  if(env->dish_size)
  {
    dish_size = env->dish_size;
    st_report("fwhm(): Using Non Standard Dish Size: %.2f M", dish_size);
  }

  if(env->taper)
  {
    taper = env->taper;
    st_report("fwhm(): Using Non Standard Taper: %.2f", taper);
  }

  st_report("fwhm(): Dish Size     = %5.2f M", dish_size);
  st_report("fwhm(): Taper Factor  = %5.2f", taper);

  lambda = SPL / freq;
  theta = taper * (lambda / dish_size);
  hm  = theta * (180.0 / M_PI) * 1800.0;

  st_report("fwhm(): lambda        = %5.2f mm",   lambda * 1000.0);

  return(hm);
}


double findObj(name)
char *name;
{
  FILE *fp;
  int day, mon, n;
  int yr, curyea, curmon, curday;
  double ven_d, mar_d, jup_d, sat_d;
  double ow;
  char tbuf[256], ename[256];

  curyea = getYear(0);
  curmon = getMonth(0);
  curday = getDom(0);

  sprintf(ename, "%s/share/%s", pkg_home, EPHEM);

  if((fp = fopen(ename, "r")) == NULL)
  {
    st_fail("Can't Open Ephemeris File: %s", ename);

    return(-1.0);
  }

     /* read line matching input date and time plus the one after that */
     /* (skip the sun-angle data fields in the ephemeris)              */
/*
 YR  M  D    Vd     ?      Mdia  Mdist  ?        Jdia   ?        Sdia  Sinc  ?
2015 1  21  10.76  21.2    4.53  1.392  36.1    43.45  160.8    15.11  24.8  57.9
  0  1  2    3      4       5      6     7        8     9        10     11   12
 */

  while(fgets(tbuf, sizeof(tbuf), fp) != NULL)
  {
    n = convol_get_tok(tbuf);

    if(n == 13)
    {
      yr  = atoi(convol_tok[0]);
      mon = atoi(convol_tok[1]);
      day = atoi(convol_tok[2]);

      if(yr == curyea && mon == curmon && day == curday)
      {
        ven_d = atof(convol_tok[ 3]);
        mar_d = atof(convol_tok[ 5]);
        jup_d = atof(convol_tok[ 8]);
        sat_d = atof(convol_tok[10]);

        break;
      }
    }
  }

  fclose(fp);

  st_report("Venus         = %.2f, Mars = %.2f, Jup = %.2f, Sat = %.2f",
        ven_d, mar_d, jup_d, sat_d);

  if(!strcmp(name, "Mars"))
  {
    ow = mar_d;
  }
  else
  if(!strcmp(name, "Jupiter"))
  {
    ow = jup_d;
  }
  else
  if(!strcmp(name, "Venus"))
  {
    ow = ven_d;
  }
  else
  if(!strcmp(name, "Saturn"))
  {
    ow = sat_d;
  }
  else
  {
    ow = 0.0;
  }

  return(ow);
}


void convolve(const double Signal[/* SignalLen */], size_t SignalLen,
              const double Kernel[/* KernelLen */], size_t KernelLen,
              double Result[/* SignalLen + KernelLen - 1 */])
{
  size_t n;

  for (n = 0; n < SignalLen + KernelLen - 1; n++)
  {
    size_t kmin, kmax, k;

    Result[n] = 0.0;

    kmin = (n >= KernelLen - 1) ? n - (KernelLen - 1) : 0;
    kmax = (n < SignalLen - 1) ? n : SignalLen - 1;

    for (k = kmin; k <= kmax; k++)
    {
      Result[n] += Signal[k] * Kernel[n - k];
    }
  }
}

int goodSource(s)
char *s;
{
  int i;

  for(i=0;i<MAX_GOOD_SOURCES;i++)
  {
    if(!strcasecmp(s, good_sources[i]) )
    {
      return(1);
    }
  }

  return(0);
}


int convol(source, first, freq, fw, source_size, beam_size)
int first;
char *source;
double freq, *fw, *source_size, *beam_size;
{
  int i;
  double beam_width, half_obj;
  double a[4], dyda[4], y, maxg=0.0;
  struct RESULTS *p;
  static double obj_width;

  st_report("convol(): Object        = %s",   source);
  st_report("convol(): Freq          = %10.6f GHz",   freq);

  beam_width = fwhm(freq) * 2.0;
  *fw = beam_width;

  st_report("convol(): FWHM          = %5.2f arcsecs",   beam_width);

  if(!goodSource(source))
  {
   st_report("convol(): Skipping %s", source);
   st_report("convol(): Setting planet_hp to %.2f", 3600.0 * beam_width / 2.0);

    return(1);
  }

  if(first)
  {
    obj_width    = findObj(source);
  }

  *source_size = obj_width;

  if(obj_width < 0.0) /* Something failed */
  {
    return(1);
  }

  st_report("convol(): Object Width  = %5.2f arcsecs",   obj_width);

  half_obj   = obj_width / 2.0;

  a[1] = 1.0;
  a[2] = 0.0;
  a[3] = beam_width;

// Compute a perfect gauss and object width
  p = &results[0];
  for(i=-100;i<100;i++, p++)
  {
    fgauss((double)i, a, &y, dyda);

    p->i = i;
    p->beam = y;

    if(i < (int)-half_obj || i > (int)half_obj)
    {
      p->source = 0.0;
    }
    else
    {
      p->source = 1.0;
    }
  }

// load the convol arrays
  p = &results[0];
  for(i=0;i<MAX_RESULTS;i++, p++)
  {
    b_array[i] = p->beam;
    s_array[i] = p->source;
    r_array[i] = 0.0;
  }

  convolve(b_array, MAX_RESULTS, s_array, MAX_RESULTS, r_array);

// Find the peak value for scaling later
  for(i=0;i<MAX_RESULTS;i++, p++)
  {
    if(r_array[i] > maxg)
    {
      maxg = r_array[i];
    }
  }

// Fill the fit array with the results of convolve()
  p = &results[0];
  for(i=0;i<MAX_RESULTS;i++, p++)
  {
    if(i <= 100)
    {
      p->fit = r_array[i+100];
    }
    else
    {
      p->fit = r_array[MAX_RESULTS-(i-100)];
    }

// Scale to ~1.0
    p->fit = p->fit / maxg;
  }

// Find the peak value again as there may be points slightly
// greater then 1.0; for scaling later
  maxg = 0.0;
  p = &results[0];
  for(i=0;i<MAX_RESULTS;i++, p++)
  {
    if(p->fit > maxg)
    {
      maxg = p->fit;
    }
  }

// Fit a gauss to the resulting beam to find useful FWHM for pointing
  a[1] = 1.0;
  a[2] = 100.0;
  a[3] = 10.0;

  p = &results[0];
  for(i=0;i<MAX_RESULTS;i++, p++)
  {
    x_array[i] = (double)i;
    y_array[i] = p->fit;
    s_array[i] = 1.0;
  }

  if(!gauss_fit(x_array, y_array, s_array, MAX_RESULTS-1, a, 0)) /* last flag is verbose */
  {
    st_fail("gauss_fit() returned failure");
    return(1);
  }

  st_report("convol(): GaussFit      = %.2f %.2f %.2f", a[1], a[2], a[3]);

  *beam_size = a[3];

  return(0);
}


int initDeConStruct()
{
  int i;
  double freq, fw, source_size, beam_size;
  struct DECONVOL_STRUCT *p;
  FILE *fp;

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

  st_report("initDeConStruct(): Initializing Convole Table for %s on planet: %s", 
                           env->issmt == 1 ? "SMT" : "12M", env->source);

  dcs.telescope = env->issmt;
  strcpy(dcs.planet, env->source);

  freq = 850.0;

  i = 0;
  stquiet = 1; /* Quiet from now on */

  while(freq >= 30.0)
  {
    if(fabs(freq-(restFreq(&h0->head)/1000.0)) < 0.5)
    {
      stquiet = 0;
    }

    if(convol(dcs.planet, i == 0, freq, &fw, &source_size, &beam_size))
    {
      stquiet = 0;
      return(1);
    }

    stquiet = 1; /* Quiet back on */

    dcs.freq_array[i] = freq;
    dcs.fw_array[i]   = fw;
    dcs.bs_array[i]   = beam_size;

    freq -= 1.0;
    i++;
    p++;
  }

  stquiet = 0;

  dcs.count       = i;
  dcs.source_size = source_size;

  st_report("initDeConStruct(): Table Size    = %4d", dcs.count);

  if((fp = fopen("convol.table", "w")) == NULL)
  {
    st_fail("initDeConStruct(): Can't Open Convol Table File: %s", "convol.table");

    return(1);
  }

  for(i=0;i<dcs.count;i++)
  {
    fprintf(fp, "%10.5f %10.6f %10.6f %10.6f\n", 
                dcs.freq_array[i], dcs.fw_array[i], dcs.source_size, dcs.bs_array[i]);
  }

  fclose(fp);

  dcs.initialized = 1;

  return(0);
}


int doDeConvol()
{
  int next, fail = 0;
  double formula_bs, source_size2, con2;

  if(!dcs.initialized || strcmp(env->source, dcs.planet) || dcs.telescope != env->issmt)
  {
    initDeConStruct();
  }

  st_report("doDeConvol(): Searching table BS =   %.2f", beam_map.gmap_bw1);

  decon_beam_size = interpl8( &dcs.bs_array[0], &dcs.fw_array[0], dcs.count, beam_map.gmap_bw1, &next );

  if(next < dcs.count)
  {
    st_report("doDeConvol(): Found @ Table Index: = %d", next);
    st_report("doDeConvol(): Decon Table Result   = %.2f", decon_beam_size);
  }
  else
  {
    st_fail("doDeConvol(): Decon Function ran off end of table");
    fail = 1;
  }

  /* Now compute deconvolved beam using a formula */
  con2 = beam_map.gmap_bw1*beam_map.gmap_bw1;
  source_size2 = dcs.source_size*dcs.source_size;

  formula_bs = sqrt(con2 - (log(2.0) / 2.0) * source_size2);
  st_report("doDeConvol(): Decon Formula Result = %.2f", formula_bs);

  if(fail)
  {
    decon_beam_size = formula_bs;
  }

  return(0);
}


/* Return the interpolated value in the table of "no" points, (xtab, ytab), 
 * that corresponds to the xtab value, "x".  "xtab" must be stored in increasing
 * order. "next" is returned such that xtab[next] < x < xtab[next+1].
 */

double interpl8( xtab, ytab, no, x, next )
double *xtab;
double *ytab;
int no;
double x;
int *next;
{
  int i;

  if(x >= xtab[0])
  {
    for(i=0;i<no-1 && xtab[i+1]<x;i++)
    {
      continue;
    }

    if(i >= no-1)
    {
      i = no - 2;
    }
  }
  else
  {
    for(i=0;i>0 && xtab[i]>x;i--)
    {
      continue;
    }

    if(i < 0)
    {
      i = 0;
    }
  }

  *next = i;

  return(ytab[i] + (ytab[i+1] - ytab[i]) * (x - xtab[i]) /
                   (xtab[i+1] - xtab[i]));
}