beam.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   MAX_BKNDS  128
#define   H          6.626e-34          /* Planck's constant */
#define   K          1.381e-23          /* Boltzmann's constant */
#define   C          0.299792458        /* speed of light (m/ns) */
#define   LN2        0.693147181        /* ln(2.0) */
#define   CMB        2.725              /* Cosmic Microwave Bkgd */
#define   R0         1.524              /* Mars mean heliocentric distance */
#define   T_R0       206.8              /* Mars temp at R0 */
#define   P          29.530589          /* Period of Moon */


struct RESULTS {
        double avg;
        double rms;
        double eta_mb;
        double err;
        double rej;
        struct HEADER head;
};

struct RESULTS results[MAX_BKNDS][MAX_BKNDS];
struct RESULTS averages[MAX_BKNDS];
double tmpBeamEff[4];

char *months[12] = {"Jan", "Feb", "Mar", "Apr", "May", "Jun",
                    "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};

char *planet_names[MAX_PLANETS] = {"Mercury", "Venus", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", "Moon"};

double parseBeamEff(buf)
char *buf;
{
  double val;

  val = atof(buf);

  if(val < 0.01)
  {
    st_warn("Beam Eff %f out of range: 0.0 - 1.0, or 1.0 - 100.0", val);
    return(1.0);
  }

  if(val> 1.0)
  {
    val /= 100.0;
  }

  if(val < 0.01 || val > 1.00)
  {
    return(1.0);
  }

  return(val);
}


int setBeamEff()
{
  int i;

  if(strlen(args->secondArg))
  {
    if(strncmp(args->secondArg, "-", 1))
    {
      env->bm_eff[0] = parseBeamEff(args->secondArg);
    }
  }
  else
  {
    st_print("\nCurrent Beam Eff Values: ");

    for(i=0;i<4;i++)
    {
      st_print("%.2f ", env->bm_eff[i]);
    }

    st_warn("");

    return(1);
  }

  if(strlen(args->thirdArg))
  {
    if(strncmp(args->thirdArg, "-", 1))
    {
      env->bm_eff[1] = parseBeamEff(args->thirdArg);
    }
  }

  if(strlen(args->fourthArg))
  {
    if(strncmp(args->fourthArg, "-", 1))
    {
      env->bm_eff[2] = parseBeamEff(args->fourthArg);
    }
  }

  if(strlen(args->fifthArg))
  {
    if(strncmp(args->fifthArg, "-", 1))
    {
      env->bm_eff[3] = parseBeamEff(args->fifthArg);
    }
  }

  return(0);
}


int setRejections()
{
  int i;

  if(strlen(args->secondArg))
  {
    if(strncmp(args->secondArg, "-", 1))
    {
      env->rejections[0] = atoi(args->secondArg);
    }
  }
  else
  {
    st_print("Current Rejection Values: ");

    for(i=0;i<MAX_BEAMS;i++)
    {
      st_print("%.2f ", env->rejections[i]);
    }

    st_print("\n");

    return(1);
  }

  if(strlen(args->thirdArg))
  {
    if(strncmp(args->thirdArg, "-", 1))
    {
      env->rejections[1] = atoi(args->thirdArg);
    }
  }

  if(strlen(args->fourthArg))
  {
    if(strncmp(args->fourthArg, "-", 1))
    {
      env->rejections[2] = atoi(args->fourthArg);
    }
  }

  if(strlen(args->fifthArg))
  {
    if(strncmp(args->fifthArg, "-", 1))
    {
      env->rejections[3] = atoi(args->fifthArg);
    }
  }

  return(0);
}


int doAcceptEff()
{
  int i, tot=0;

  for(i=0;i<4;i++)
  {
    if(tmpBeamEff[i] > 0.0)
    {
      st_print("Accepting Beam Eff Value %.1f%% for Chan #%d\n", tmpBeamEff[i], i+1);
      env->bm_eff[i] = tmpBeamEff[i] / 100.0;
      tot++;
    }
  }

  if(!tot)
  {
    st_warn("No Saved Beam Eff Values to Accept");
  }

  return(0);
}


struct tm *breakTime(burnT)
double burnT;
{
  static struct tm *tmn1;
  long sec;

  sec = (long)burnT + 7 * 3600;
  tmn1 = localtime(&sec);

  return(tmn1);
}


int findBackend(scan)
double scan;
{
  int indx = -1;
  double fscan;

  fscan = scan - (int)scan;

  fscan *= 100.0;

  if(fabs(fscan-1.0) < 0.001)
  {
    indx = 0;
  }
  else
  if(fabs(fscan-2.0) < 0.001)
  {
    indx = 1;
  }
  else
  if(fabs(fscan-3.0) < 0.001)
  {
    indx = 2;
  }
  else
  if(fabs(fscan-4.0) < 0.001)
  {
    indx = 3;
  }

  return(indx);
}

int beam_eff(first)
int first;
{
  int    planet;
  double f_sky, rej, T_ant, ta_s, Tmb;
  double T_src, f_rej, theta_p, R_mars, incl;
  double lambda, eta_mb, eta_hi, eta_lo, theta_mb;
  double T0, chi, xi, ts_s, denom, err, diameter;
  STATE state, *st;
  EPH eph;

  eph       = global.eph;
  state     = global.state;
  st        = &global.state;
  planet    = state.planet; /* Planet indx */
  f_sky     = state.f_sky;  /* Observed freq in GHz */
  rej       = state.rej;    /* Sideband rejection in dB */
  T_ant     = state.T_ant;  /* AVG temp of observations */
  ta_s      = state.ta_s;   /* RMS temp of observation */

     /* compute T_src and set theta_p depending on planet */
  switch (planet)
  {
    case JUPITER:

               /* from Shirley: best linear fit to Mangum Table 1 */
               /* with additional data points from WMAP & MASM1   */
      T_src = 178.309 - 0.036 * f_sky;
      ts_s = 10.0;
      theta_p = eph.JupDiam;
      break;

    case VENUS:

               /* from Shirley: logistic fit to brightness temperature data */
      T_src = 266.4 + (245.1 / (1.0 + (f_sky*f_sky / (64.7*64.7))));
      ts_s = 10.0;
      theta_p = eph.VenDiam;
      break;

    case MARS:

               /* from Ulich (1981) as discussed by Mangum */
      R_mars = eph.MarDist;
      T_src = T_R0 * sqrt(R0 / R_mars);
      ts_s = 6.0;
      theta_p = eph.MarDiam;
      break;

    case SATURN:

               /* from Shirley: interpolated planetary flux densities */
               /* (Table 6.2 posted in 12-meter control room          */
      incl = eph.SatIncl;
      T_src = (f_sky - 90.0) * (149.0 - 140.0) / (90.0 - 227.0);
      T_src += 149.0;
      ts_s = 10.0;
      theta_p = eph.SatDiam;
      break;

    case MOON:
 
               /* from Krotikov and Troitskii (1964) and Linsky (1966,1973) */
               /* as discussed by Mangum                                    */

      st_print("*** CAUTION:  Moon option for test purposes only! ***\n");
       
      lambda = 1000.0 * C / f_sky;

      chi = (eph.MoonAge / P) * 360.0 - 180.0;
      chi *= PI / 180.0;

      xi = atan(0.24 * lambda / (1.0 + 0.24 * lambda));
      xi *= PI / 180.0;

               /* from Begam: linear fit to Table 2 data */

      T0 = 220.02682849 - 0.000348260 * f_sky;  

               /* Mangum eqn (2) */

      T_src = T0 + 0.77 * T0 / sqrt(1.0 + 0.48*lambda + 0.11*lambda*lambda);
      T_src *= cos(chi - xi);
    
               /* token value for ts_s -- not actually determined! */

      ts_s = 20.0;
      break;

    default:  

      st_fail("invalid planet selection");
      return(0);
  }
 
  
     /* telescope main-beam FWHM -- assumes 10 dB edge taper */

  lambda = C / f_sky;

  if(!strncmp(h0->head.telescop, "ARO-HHT", 7))
  {
    diameter = MG_DIAMETER;
  }
  else
  {
    diameter = KP_DIAMETER;
  }

  if(first)
  {
    st_report("Using a Telescope Diameter of %.1f Meters", diameter);
  }

  theta_mb = 1.17 * (lambda / diameter) * (180.0 / PI) * 3600.0;

     /* "rejection factor" per Schlingman et al 2011,ApJS,195,14  */
     /*  replaces factor of 1/2 (for dual SB) in Mangum eqn (1)   */

  f_rej = 1.0 / (1.0 + pow(10.0, (-rej/10.0)));

     /* Scale at the 12-meter is actually T*_R.  Multiply this */
     /* by Eta_fss to convert back to T*_A                     */
     /*                                                        */
     /* T_ant *= 0.68;                                         */
     /*                                                        */
     /* No longer the case as of 2014                          */
     /*                                                        */

     /* main-beam efficiency from Mangum eqn (1) */

  eta_mb = f_rej * T_ant / (planck(f_sky, T_src) - planck(f_sky, CMB)); 
  eta_hi = f_rej * (T_ant + ta_s) / (planck(f_sky, T_src - ts_s) - planck(f_sky, CMB)); 
  eta_lo = f_rej * (T_ant - ta_s) / (planck(f_sky, T_src + ts_s) - planck(f_sky, CMB)); 

     /* convolve with source disk except in the case of the Moon */

  denom = (1.0 - exp((-LN2 * theta_p * theta_p) / (theta_mb * theta_mb))); 

  Tmb = T_src * denom;

  if(planet != MOON)
  { 
    eta_mb /= denom;
    eta_hi /= denom;
    eta_lo /= denom;
  }

  st->eta_mb = eta_mb;

     /* compute the uncertainty */

  err = (eta_hi - eta_lo) / 2.0;
  st->err = err;

     /* outputs */

  if(first)
  {
    st_report("Brightness Temperature for %s = %6.2lf K", planet_names[planet], T_src);
  }

  if(planet == MOON)
  {
    st_report("Days Since New Moon = %6.2lf", eph.MoonAge);
  }
  else 
  {
    if(first)
    {
      st_report("Projected Diameter for %s = %6.2lf arcsecs", planet_names[planet], theta_p);
    }
  }

  if(planet == MARS)
  {
    if(first)
    {
      st_report("Heliocentric Distance for Mars = %5.3lf a.u.", R_mars);
    }
  }

  if(planet == SATURN)
  {
    if(first)
    {
      st_report("Inclination of Rings = %4.1lf degs", incl);
    }
  }

  if(first)
  {
    st_report("Telescope Main-Beam Temp = %4.1lf K", Tmb);
    st_report("Telescope Main-Beam FWHM = %4.1lf arcsecs", theta_mb);
  }

  st_report("Rejection Factor = %7.5lf", f_rej);

  st_report("Main-Beam Efficiency = %5.2lf +/- %.2lf%%", 
                    (st->eta_mb * 100.0), (st->err * 100.0));

  return(0);
}


/* Planck Function */

double planck(double freq, double T_bright)
{
     double p;

     p = (H * freq / K) / (exp((H * freq) / (K * T_bright)) - 1.0);

//     printf("Planck(): %f, %f, %f\n", freq, T_bright, p);
     return p;
}



int processResults(n, indx, fp)
int n, indx;
FILE *fp;
{
  int i, nif;
  struct RESULTS *a, *p;
  static struct tm *tmn3;

  if(n)
  {
    bzero((char *)&averages, sizeof(averages));

    a = &averages[0];

    for(i=0;i<n;i++,p++)
    {
      p = &results[indx][i];

      switch(n) /* Copy the appropriate header */
      {
        case 3:
          if(i == 1)
          {
            bcopy((char *)&p->head, (char *)&a->head, sizeof(struct HEADER));
            a->rej = p->rej;
          }
          break;

        case 2:
          if(i == 0)
          {
            bcopy((char *)&p->head, (char *)&a->head, sizeof(struct HEADER));
            a->rej = p->rej;
          }
          break;

        default:
          bcopy((char *)&p->head, (char *)&a->head, sizeof(struct HEADER));
          a->rej = p->rej;
          break;
      }
      a->avg    += p->avg;
      a->rms    += p->rms;
      a->eta_mb += p->eta_mb;
      a->err    += p->err;
    }

    a->avg    /= n;
    a->rms    /= n;
    a->eta_mb /= n;
    a->err    /= n;

    nif = getFeed(a->head.scan);

    tmn3 = breakTime(a->head.burn_time);

/* Write out the log file beam_eff.dat */
    fprintf(fp, "%02d %3s %4d   %02d:%02d  (%d) %7.2f  %3.3s  %7.7s %10.6lf   %02d   %6.2lf  %4.1lf  %.1f  %.0f\n",

        tmn3->tm_mday, months[tmn3->tm_mon], tmn3->tm_year+1900,
        tmn3->tm_hour, tmn3->tm_min, n, a->head.scan,
        a->head.obsid,  a->head.object,
        restFreq(&a->head) / 1000.0, nif,
        a->eta_mb*100.0, a->err*100.0,
        a->head.el, a->rej);

    st_report("IF%2d: Avg = %.2f, RMS = %.3f, Eff = %.2f, Err = %.2f",
                nif, a->avg, a->rms, a->eta_mb*100.0, a->err*100.0);

    tmpBeamEff[nif-1] = a->eta_mb*100.0;
  }

  return(0);
}



int writeResults(indx)
int indx[4];
{
  FILE *fp;
  int i;

  st_report(" ");
  st_report("Averages: %d %d %d %d", indx[0], indx[1], indx[2], indx[3]);
  st_report(" ");

  if((fp = fopen("beam_eff.dat", "a")) == NULL)
  {
    printf("BEAM_EFF: Unable to open beam_eff.dat file\n");
    exit(1);
  }

  fprintf(fp, "#\n");
  fprintf(fp, "Analysed: %s\n", getcTime(0));
  fprintf(fp, "#\n");

  for(i=0;i<4;i++)
  {
    processResults(indx[i], i, fp);
  }

  fclose(fp);

  return(0);
}


int doBeamEff()
{
  int first = 1, i, n, bknd, bkindx[MAX_BKNDS], found = 0, saveCalMode;
  double ret;
  struct HEADER *p;
  struct RESULTS *p2;
  static struct tm *tmn2;
  struct ACCUM_RESULTS *ar;

  bzero((char *)&tmpBeamEff[0], sizeof(tmpBeamEff));
  bzero((char *)&bkindx,        sizeof(bkindx));
  bzero((char *)&accum_results, sizeof(accum_results));

  saveCalMode = env->calmode;

  if(env->issmt)
  {
    if(env->rejections[0] == 0.0 || env->rejections[1] == 0.0 || 
       env->rejections[2] == 0.0 || env->rejections[3] == 0.0)
    {
      st_fail("One or More Rejections are not set");
      setRejections();

      return(1);
    }
  }
  else
  if(env->rejections[0] == 0.0 || env->rejections[1] == 0.0)
  {
    st_fail("One or More Rejections are not set");
    setRejections();

    return(1);
  }

  p = &sa[0];
  for(i=0;i<env->numofscans;i++,p++)
  {
    if(meetsCritiria(p, CHECK_ALL))
    {
      if(strncmp(p->obsmode, "CONTSEQ", 7)) /* We only process Seq */
      {
        continue; /* No good */
      }

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

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

      if(ret > 0.0)
      {
        found++;

        tmn2 = breakTime(h0->head.burn_time);

        global.time.year       = tmn2->tm_year+1900;
        global.time.mon        = tmn2->tm_mon+1;
        global.time.day        = tmn2->tm_mday;
        global.time.hour       = tmn2->tm_hour;
        global.time.min        = tmn2->tm_min;
        global.time.UseCurrent = 0; /* Use these values */

        if(first)
        {
          st_report("Obs-Time: Mon %d, Day %d, Year %d, Hour %d, Min %d",
                tmn2->tm_mon+1,
                tmn2->tm_mday,
                tmn2->tm_year+1900,
                tmn2->tm_hour,
                tmn2->tm_min);

/* Note to self: Remember *tmnow is a pointer to a static structure 
                 inside time(), and localtime(); Once you call them
                 you must use them imediately befor calling other
                 functions that may or not use time() and localtime(),
                 thus resetting their internal static structure.
 */
        }

        getephem(first);

        if(!strncmp(h0->head.object, "Jupiter", 7))
        {
          global.state.planet = JUPITER;
        }
        else
        if(!strncmp(h0->head.object, "Venus", 5))
        {
          global.state.planet = VENUS;
        }
        else
        if(!strncmp(h0->head.object, "Mars", 4))
        {
          global.state.planet = MARS;
        }
        else
        if(!strncmp(h0->head.object, "Saturn", 6))
        {
          global.state.planet = SATURN;
        }
        else
        if(!strncmp(h0->head.object, "Moon", 4))
        {
          global.state.planet = MOON;
        }
        else
        {
          st_fail("Unknown Object %s", h0->head.object);
          continue;
        }

        n = h0->head.datalen / sizeof(float);

        env->calmode     = 0;
        env->accum       = 0;

        bknd = findBackend(h0->head.scan);

        if(bknd == -1)
        {
          st_fail("Unable to determine backend %8.8s", h0->head.backend);
          env->calmode = saveCalMode;

          return(1);
        }

        ar = &accum_results[bknd];

        accum_v1(h0->yplot1, &h0->head, env->invert, ar);

        avg = env->accumtr;
        rms = 0.0;

        st_report("Nchan = %d, AVG = %f, RMS = %f", n, avg, rms);

        st_report("Backend %8.8s = %d, Rej = %.1f",
                                        h0->head.backend, bknd, env->rejections[bknd]);

        global.state.f_sky  = restFreq(&h0->head) / 1000.0;
        global.state.rej    = env->rejections[bknd];
        global.state.T_ant  = avg;
        global.state.ta_s   = rms;
        strncpy(global.state.obs, h0->head.obsid, 3);
        global.state.obs[3] = '\0';

        beam_eff(first);

        /* Save the results */
        if(bknd < MAX_BKNDS)
        {
          if(bkindx[bknd] < MAX_BKNDS)
          {
            p2 = &results[bknd][bkindx[bknd]];

            bcopy((char *)&h0->head, (char *)&p2->head, sizeof(struct HEADER));
            p2->avg    = avg;
            p2->rms    = rms;
            p2->rej    = env->rejections[bknd];
            p2->eta_mb = global.state.eta_mb;
            p2->err    = global.state.err;
          }
          else
          {
            printf("Number of Scans %d Out Of Range\n", bkindx[bknd]);
            env->calmode = saveCalMode;

            return(1);
          }

          bkindx[bknd]++;
        }
        else
        {
          printf("Backend Indx %d Out Of Range\n", bknd);
          env->calmode = saveCalMode;

          return(1);
        }

        st_report(" ");
        first = 0;
      } 
    }
  }


  if(found)
  {
    writeResults(bkindx);
    st_report("Results written to beam_eff.dat");
  }

  st_report(" ");
  st_report("Found %d scans matching your criteria", found);

  if(found)
  {
    st_report("\n");

    if(getYesNoQuit("Use beam eff value of %.3f %.3f. %.3f %.3f: y/N ?",
                        tmpBeamEff[0], tmpBeamEff[1],
                        tmpBeamEff[2], tmpBeamEff[3]))
    {
      doAcceptEff();
    }
  }

  env->calmode = saveCalMode;

  return(0);
}