planets.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 ONE_HOUR (60 * 60)  /* Seconds in one hour */

struct NEW_EPHEM   new_ephem;

static char eph_tok[MAX_TOKENS][MAX_TOKENS];

char marsDate[80];

long fluxTablesLoadTime = -999999999;

double  mars_freqs[] = { 30.0,  80.0, 115.0, 150.0, 200.0, 230.0,  260.0, 300.0,
                        330.0, 360.0, 425.0, 650.0, 800.0, 950.0, 1000.0
};

double planet_value;

struct TBRIGHT venus_tbright[MAX_VENUS_ENTRIES];
struct TBRIGHT mars_tbright[MAX_MARS_ENTRIES];
struct TBRIGHT jupiter_tbright[MAX_JUPITER_ENTRIES];
struct TBRIGHT uranus_tbright[MAX_URANUS_ENTRIES];
struct TBRIGHT neptune_tbright[MAX_NEPTUNE_ENTRIES];
struct TBRIGHT vesta_tbright[MAX_VESTA_ENTRIES];

struct EPH_OBJECT save_planet_values[MAX_PLANETS];

struct FLUX_INFO flux_info;


struct PLANET_LOOKUP {
        char *vname;
        char *fmt;
        int   type;
        void *offset;
};

#define ZRO ((struct EPH_OBJECT *)0)

struct PLANET_LOOKUP planet_lookup[] = {
        { "name",        "%s",          CHARSTAR,       &ZRO->name },
        { "dbase",       "%s",          CHARSTAR,       &ZRO->dbase },
        { "ra",          "%.10f",       DOUBLE,         &ZRO->ra },
        { "dec",         "%.10f",       DOUBLE,         &ZRO->dec },
        { "illum",       "%.10f",       DOUBLE,         &ZRO->illum },
        { "dmajor",      "%.10f",       DOUBLE,         &ZRO->dmajor },
        { "dminor",      "%.10f",       DOUBLE,         &ZRO->dminor },
        { "sdist",       "%.10f",       DOUBLE,         &ZRO->sdist },
        { "rdot",        "%.10f",       DOUBLE,         &ZRO->rdot },
        { "dist",        "%.10f",       DOUBLE,         &ZRO->dist },
        { "vel",         "%.10f",       DOUBLE,         &ZRO->vel },
        { "sangle",      "%.10f",       DOUBLE,         &ZRO->sangle },
        { "t_src",       "%.10f",       DOUBLE,         &ZRO->tsrc },
        { "trj",         "%.10f",       DOUBLE,         &ZRO->trj },
        { "flux",        "%.10f",       DOUBLE,         &ZRO->flux },
        { "flux2",       "%.10f",       DOUBLE,         &ZRO->flux2 },
        { "tmb",         "%.10f",       DOUBLE,         &ZRO->tmb },
        { "csize", "%.10f",     DOUBLE,         &ZRO->csize },

        { NULL,           NULL,              0,          NULL }
};


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

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

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

  return(c);
}

int parseLine(obj, name)
struct EPH_OBJECT *obj;
char *name;
{
  double hr, min, sec, deg;

  strcpy(obj->name, name);

  hr  = atof(eph_tok[2]);
  min = atof(eph_tok[3]);
  sec = atof(eph_tok[4]);
  obj->ra = hr + min / 60.0 + sec / 3600.0;

  deg = atof(eph_tok[5]);
  min = atof(eph_tok[6]);
  sec = atof(eph_tok[7]);
  obj->dec = fabs(deg) + min / 60.0 + sec / 3600.0;
  if(deg < 0.0)
  {
    obj->dec *= -1.0;
  }

  obj->illum     = atof(eph_tok[8]);
  obj->dmajor    = atof(eph_tok[9]);
  obj->sdist  = atof(eph_tok[10]);
  obj->rdot      = atof(eph_tok[11]);
  obj->dist  = atof(eph_tok[12]);
  obj->vel       = atof(eph_tok[13]);
  obj->sangle = atof(eph_tok[14]);

  return(0);
}


/* Interalate between to sets of positions and params */
int eph_interp(ephem, p1, p2)
struct NEW_EPHEM *ephem;
struct EPH_OBJECT *p1, *p2;
{
  double frac;
  struct EPH_OBJECT *p;
  double eccent[] = { 1.0, 1.0, 0.9952, 0.93512, 0.892467, 0.97118, 0.97247, 1.0, 1.0, 1.0, 1.0 };

  p = &ephem->planets[ephem->planet];

  frac = (ephem->time.hour + (ephem->time.min / 60.0)) / 24.0;

  strcpy(p->name, p1->name);

  p->ra     = p1->ra     + frac * (p2->ra     - p1->ra);
  p->dec    = p1->dec    + frac * (p2->dec    - p1->dec);
  p->illum  = p1->illum  + frac * (p2->illum  - p1->illum);
  p->dmajor = p1->dmajor + frac * (p2->dmajor - p1->dmajor);
  p->sdist  = p1->sdist  + frac * (p2->sdist  - p1->sdist);
  p->rdot   = p1->rdot   + frac * (p2->rdot   - p1->rdot);
  p->dist   = p1->dist   + frac * (p2->dist   - p1->dist);
  p->vel    = p1->vel    + frac * (p2->vel    - p1->vel);
  p->sangle = p1->sangle + frac * (p2->sangle - p1->sangle);
  p->dminor = p->dmajor  * eccent[ephem->planet];
  p->frac   = frac;

  return(0);
}


int parseNewEphem(ephem)
struct NEW_EPHEM *ephem;
{
  int n, indx=0;
  FILE *fp;
  char fname[256], *string, line[256], dstr[256];
  extern char *months[12];
  struct EPH_OBJECT p1, p2;

  bzero((char *)&p1, sizeof(p1));
  bzero((char *)&p2, sizeof(p2));

  sprintf(fname, "%s/share/%s-Ephem.tab", pkg_home, eph_source_names[ephem->planet]);

  if((fp = fopen(fname, "r") ) <= (FILE *)0)
  {
    st_fail("Unable to open file %s", fname);

    return(0);
  }

  sprintf(dstr, "%4d-%s-%02d", ephem->time.year, months[ephem->time.mon], ephem->time.day);

//  st_report("Searching for Date: %s", dstr);

  while((string = fgets(line, sizeof(line), fp)) != NULL)
  {
    indx++;

    if(strlen(line) < 4)
    {
      continue;
    }

    if(line[0] == '#')
    {
      continue;
    }

    n = eph_get_tok(line);

    if(n == 16)
    {
      if(!strcmp(eph_tok[0], dstr))
      {
        parseLine(&p1, eph_source_names[ephem->planet]);

//        st_report("Found Date on Line: %d", indx);

        /* Now read the nextline for interp */
        if((string = fgets(line, sizeof(line), fp)) != NULL)
        {
          n = eph_get_tok(line);

          if(n == 16)
          {
            parseLine(&p2, eph_source_names[ephem->planet]);
            eph_interp(ephem, &p1, &p2);

            if(!strcmp(eph_source_names[ephem->planet], "Mars"))
            {
              st_report("Mars Ephem line = %.2f", (double)indx + ephem->planets[ephem->planet].frac);
            }

            fclose(fp);

            return(0);
          }
        }
      }
    }
    else
    {
      st_report("Parsed %d tokens from %s", n, basename(fname));
    }
  }

  st_fail("Date: %s, not found in %s", dstr, basename(fname));

  fclose(fp);

  return(1);
} 


double rjot(t, f)
double t, f;
{
  double rjt;
  double a;

        // Calculation of the Brightness-Temperature from the 
        // Physical Planck Temperature:
//  bolt  = 1.38054e-23       ! J/K
//  plan  = 6.6256e-34        ! W s**2

  a    = plancK * f / boltzmann;
  rjt  = a / (exp( a / t) -1.0);    // = T_b [K]

  return(rjt);
}

        //Derive a linear equation from two points and calculate the y (=mx+b)
        //corresponding to a x
double tb_lin(f, f1, tb1, f2,tb2)
double f, f1, tb1, f2, tb2;
{
  double tb, xslope, yaxis;

  xslope = (tb1-tb2) / (f1-f2);
  yaxis  = tb1 - xslope * f1;
  tb     = xslope * f + yaxis;

  return(tb);
}



int loadPlantFluxTables(t)
struct EPH_TIME   *t;
{
  int i, n, len, indx=0;
  FILE *fp;
  char fname[256], *string, line[256], dstr[256];
  struct TBRIGHT *p;
  double scale = 1.0;

  if(t->year < 2012) /* Tables only go back to 2012 */
  {
    return(0);
  }

  st_report("Loading Planet Brightness Tables");

  for(i=1;i<7;i++)
  {
    if(i==4 || i == 2) /* Skip Saturn; and Mars for now */
    {
      continue;
    }

    sprintf(fname, "%s/share/%s_Tb.dat", pkg_home, eph_source_names[i]);

    if((fp = fopen(fname, "r") ) <= (FILE *)0)
    {
      st_fail("Unable to open file %s", fname);

      continue;
    }

    switch(i)
    {
      case 1: p = &venus_tbright[0];   scale = 1.0;    break;
      case 3: p = &jupiter_tbright[0]; scale = 1.0562; break;
      case 5: p = &uranus_tbright[0];  scale = 1.0;    break;
      case 6: p = &neptune_tbright[0]; scale = 1.0;    break;
      default:
        st_fail("Unknown Planet index %d", i);
        continue;
    }

    while((string = fgets(line, sizeof(line), fp)) != NULL)
    {
      if(line[0] == '#')
      {
        continue;
      }

      n = eph_get_tok(line);

      if(n == 2)
      {
        p->freq = atof(eph_tok[0]);
        p->tb = atof(eph_tok[1]) * scale;

        p++;
      }
    }
  }

        // Mars is a different animal;
        // Only load the entry for todays date.

  sprintf(fname, "%s/share/Mars_Tb.dat", pkg_home);
  sprintf(dstr, "%4d %02d %02d %02d", t->year, t->mon+1, t->day, t->hour);

  st_report("Looking for Mars date of: %s", dstr);

  len = strlen(dstr);

  if((fp = fopen(fname, "r") ) <= (FILE *)0)
  {
    st_fail("Unable to open file %s", fname);

    fluxTablesLoadTime = time(NULL);

    return(1);
  }

  while((string = fgets(line, sizeof(line), fp)) != NULL)
  {
    indx++;

    if(line[0] == '#')
    {
      continue;
    }

    if(!strncmp(line, dstr, len))
    {
//      st_print("MARS: %s", line);

      strcpy(marsDate, dstr);

      n = eph_get_tok(line);

      if(n == 20) /* Good Read */
      {
        for(i=0;i<MAX_MARS_ENTRIES;i++)
        {
          mars_tbright[i].freq =   mars_freqs[i];
          mars_tbright[i].tb   =   atof(eph_tok[5+i]);
        }

        fclose(fp);

        st_report("Mars date line  = %d", indx);

        fluxTablesLoadTime = time(NULL);

        return(0);
      }
      else
      {
        st_fail("Got %d token in file %s", n, fname);
      }
    }
  }

  /* If we got this far, we missed something */

  st_fail("Unable to locate %s, in %s", dstr, fname);

  fclose(fp);

  fluxTablesLoadTime = time(NULL);

  return(0);
}


double tsrcInterp(tb, f, n)
struct TBRIGHT *tb;
double f;
int n;
{
  double diff, min, t=0.0, slope;
  int i, len, indx=0;
  struct TBRIGHT *p;
  double f1, f2, t1, t2, frac, dy, dx;

  len = n / sizeof(struct TBRIGHT);

  if(interpl8_debug)
  {
    st_report("Size of array: %d", len);
    st_report("Looking for %f", f);
  }

  min = 999999.9;
  p = tb;
  for(i=0;i<len;i++,p++)
  {
    diff = fabs(f - p->freq);

    if(diff < min)
    {
      min = diff;
      indx = i;
    }
  }

  t  = tb[indx].tb;

  f1 = tb[indx].freq;

  if(interpl8_debug)
  {
    st_report("Found Min freq @ indx %d; min = %f, freq = %f, tb = %f", indx, min, f1, t);
  }

  /* Now interpl8 */
  if(f1 < f) /* We found the entry before our request freq */
  {
    f1 = tb[indx].freq;
    f2 = tb[indx+1].freq;

    t1 = tb[indx].tb;
    t2 = tb[indx+1].tb;
  }
  else /* We found the entry after our request freq */
  {
    f1 = tb[indx-1].freq;
    t1 = tb[indx-1].tb;

    f2 = tb[indx].freq;
    t2 = tb[indx].tb;
  }


  dy = t2 - t1;
  dx = f2 - f1;
  slope = dy / dx;
  frac = (f - f1) * slope;

  if(interpl8_debug)
  {
    st_report("T1 %6.2f  T2 %6.2f", t1, t2);
    st_report("F1 %6.2f  F2 %6.2f", f1, f2);
    st_report("DX = %f", dx);
    st_report("DY = %f", dy);
    st_report("Slope = %f", slope);
    st_report("Frac = %f", frac);
  }

  t = t1 + frac;

  return(t);
}


int fluxTableOld(t)
long t;
{
  if(abs(fluxTablesLoadTime-t) > ONE_HOUR)
  {
    return(1);
  }

  return(0);
}


int computePlanetFlux(planet, obj, t, s, logit)
int                planet, logit;
struct EPH_OBJECT *obj;
struct EPH_TIME   *t;
struct EPH_STATE  *s;
{
  int    tmb_color = BLACK;
  double diameter, alpha, taper;
  double clight, frequenz, xlambda;
  double radius, omega_source, fluss;
  double A, x, y/*, tbright*/;
  double hnuk, expterm, jnuplan, jnucmb;
  extern char *eph_source_names[];
  long tt;
#ifdef USE_OLD_WAY
  double R_mars;
  double  a0,  a1,  a2,  a3; // Uranus calculations
  double aa0, aa1, aa2, aa3; // Neptune calculations
#endif


  tt = time(NULL);

  if(fluxTableOld(tt))
  {
    loadPlantFluxTables(t);
  }

  obj->tsrc = 0.0;

  if(s->f_sky)
  {
    s->lambda = SOL / s->f_sky;
  }
  else
  {
//    st_fail("Freq set to 0.0");
    return(0);
  }

  if(env->issmt)
  {
    diameter = MG_DIAMETER;
    taper    = MG_TAPER;
  }
  else
  {
    diameter = KP_DIAMETER;
    taper    = KP_TAPER;
  }

  if(env->dish_size)
  {
    diameter = env->dish_size;
  }

  if(env->taper)
  {
    taper = env->taper;
  }

//  alpha = taper * 0.0996; /* was 0.115; 04/30/2016; twf */
//  alpha = 1.17;  /* Shirley & Folkers, Mangum 1993; for 10 db taper */
  alpha = 1.18;  /* Shirley & Folkers, for 12 db taper Oct 17th, 2016 */

  alpha = taper * (1.18 / 12.0);

  s->theta_mb = alpha * (s->lambda / diameter) * (180.0 / M_PI) * 3600.0;

  if(env->theta_mb > 0.0)       /* user supplied theta_mb */
  {
    s->theta_mb = env->theta_mb;
    tmb_color = RED;
  }

  if(logit)
  {
    st_report("Lambda:         = %6.2fmm", s->lambda*1000.0);

    if(env->dish_size) /* Flag unusual dish_size */
    {
      setColor(RED);
      st_report("Diameter:       = %6.2fm", diameter);
      setColor(BLACK);
    }
    else
    {
      st_report("Diameter:       = %6.2fm", diameter);
    }

    if(env->taper) /* Flag unusual taper */
    {
      setColor(RED);
      st_report("Taper:          = %6.2fdB", taper);
      setColor(BLACK);
    }
    else
    {
      st_report("Taper:          = %6.2fdB", taper);
    }

    st_report("Alpha:          = %6.2f", alpha);

    if(tmb_color != BLACK)
    {
      setColor(tmb_color);
      st_report("Theta_mb:       = %6.2f\"", s->theta_mb);
      setColor(BLACK);
    }
    else
    {
      st_report("Theta_mb:       = %6.2f\"", s->theta_mb);
    }
  }

  
  switch(planet)
  {
    case MERCURY:
        obj->tsrc = 450.0;
        strcpy(obj->dbase, "Fixed");
        break;

    case VENUS:
#ifdef USE_OLD_WAY
#warning Using Shirley

      /* Old Shirley: logistic fit to brightness temperature data */
      obj->tsrc = 266.4 + (245.1 / (1.0 + (s->f_sky*s->f_sky / (64.7*64.7))));
      strcpy(obj->dbase, "Shirley");
#else
//#warning Using ALMA
//      if(s->f_sky > 116.0)
//      {
        obj->tsrc = tsrcInterp(&venus_tbright[0], s->f_sky, sizeof(venus_tbright));
        strcpy(obj->dbase, "ALMA");
//      }
//      else /* From Shirley; March 2016 */
//      {
//        obj->tsrc = 362.92 * pow(s->f_sky / 90.0, -0.4055);
//        strcpy(obj->dbase, "Shirley");
//      }
#endif
      break;

    case MARS:
#ifdef USE_OLD_WAY
      /* from Ulich (1981) as discussed by Mangum */
      R_mars = obj->sdist;
      obj->tsrc = T_R0 * sqrt(R0 / R_mars);
      strcpy(obj->dbase, "Ulich");
#else
      obj->tsrc = tsrcInterp(&mars_tbright[0], s->f_sky, sizeof(mars_tbright));
      strcpy(obj->dbase, "ALMA");
#endif
      break;

   case JUPITER:
#ifdef USE_OLD_WAY
      /* from Shirley: best linear fit to Mangum Table 1 */
      /* with additional data points from WMAP & MASM1   */
      obj->tsrc = 178.309 - 0.036 * s->f_sky;

     // From Glidas/Astro: iram-planet.f90
     // 179 K at 90 GHz, 173 K at 150 GHz, 171 K at 227 GHz, 174 K at337GHz
     // inter- and extrapolate linearily:
      if(s->f_sky < 150.0)
      {
         obj->tsrc = tb_lin(s->f_sky, 90.0, 179.0, 150.0, 173.0);
      }
      else
      if(s->f_sky <= 227.0)
      {
         obj->tsrc = tb_lin(s->f_sky, 150.0, 173.0, 227.0, 171.0);
      }
      else
      if(s->f_sky > 227.0)
      {
         obj->tsrc = tb_lin(s->f_sky, 227.0, 171.0, 337.0, 174.0);
      }
#else
//      if(s->f_sky > 183.0)
//      {
        obj->tsrc = tsrcInterp(&jupiter_tbright[0], s->f_sky, sizeof(jupiter_tbright));
        strcpy(obj->dbase, "Alma");
//      }
//      else /* From Shirley; March 2016 */
//      {
//        obj->tsrc = 173.67677 * pow(s->f_sky / 90.0, 0.0071986);
//        strcpy(obj->dbase, "Shirley");
//      }
#endif
      break;

    case SATURN:
               /* from Shirley: interpolated planetary flux densities */
               /* (Table 6.2 posted in 12-meter control room          */
      obj->tsrc  = (s->f_sky - 90.0) * (149.0 - 140.0) / (90.0 - 227.0);
      obj->tsrc += 149.0;
      strcpy(obj->dbase, "Shirley");
      break;

    case URANUS:
#ifdef USE_OLD_WAY
                //  Griffin & Orton 1993, ICARUS, 105, 537 
      clight   = 2.9979e+10;         // [cm/s]
      frequenz = s->f_sky * 1.e9;    // [Hz]
      xlambda  = clight / frequenz;  // [cm]
      xlambda  = xlambda * 1.e4;     // [micrometers]

      a0 = -795.694;
      a1 =  845.179;
      a2 = -288.946;
      a3 =   35.200;
      obj->tsrc = a0 + a1 * log10(xlambda) + a2 * pow(log10(xlambda), 2.0) + a3 * pow(log10(xlambda), 3.0);
#else
      obj->tsrc = tsrcInterp(&uranus_tbright[0], s->f_sky, sizeof(uranus_tbright));
      strcpy(obj->dbase, "ALMA");
#endif
      break;

    case NEPTUNE:
#ifdef USE_OLD_WAY
                //  Griffin & Orton 1993, ICARUS, 105, 537 
      clight   = 2.9979e+10;        // [cm/s]
      frequenz = s->f_sky * 1.e9;   // [Hz]
      xlambda  = clight / frequenz; // [cm]
      xlambda  = xlambda * 1.e4;    // [micrometers]

      aa0 = -598.901;
      aa1 =  655.681;
      aa2 = -229.545;
      aa3 =   28.994;

      obj->tsrc = aa0 + aa1 * log10(xlambda) + aa2 * pow(log10(xlambda), 2.0) + aa3 * pow(log10(xlambda), 3.0);
#else
      obj->tsrc = tsrcInterp(&neptune_tbright[0], s->f_sky, sizeof(neptune_tbright));
      strcpy(obj->dbase, "ALMA");
#endif
      break;
  }

  clight   = 2.9979e+10;        // [cm/s]
  frequenz = s->f_sky * 1.e9;   // [Hz]
  xlambda  = clight / frequenz; // [cm]

#ifdef USE_RJ
  obj->trj = rjot(obj->tsrc, frequenz);
#else
  obj->trj = obj->tsrc;
#endif

        // Flux of the planet, independent of the telescope or the atmosphere:
        // Derived from the Rayleigh-Jeans corrected Brightness Temperatures!
  radius       = obj->dmajor / 2.0;                                                            // ["]
  omega_source = M_PI * pow(radius * M_PI / (180.0 * 3600.0), 2.0);                           // [sterad]
  fluss        = 2.0 * boltzmann / (pow(xlambda/100.0, 2.0)) * obj->trj * omega_source;// [Joule]
  obj->flux    = fluss * 1.e26;                                                               // [Jansky]
  A            = (M_PI * pow(diameter, 2.0) / 4.0);
  env->jfact   = ((2.0 * boltzmann / 1.0e-26) / A);

  if(logit)
  {
    st_report("Dish Area       = %6.2fm^2", A);
    st_report("JFact           = %6.2f", env->jfact);
  }

  obj->flux2  = obj->dmajor*obj->dminor*1158.1
                 * pow(s->f_sky, 3.0)
                 * (1.0 / (exp(0.04801 * s->f_sky / obj->trj)-1.0) 
                 -  1.0 / (exp(0.04801 * s->f_sky / 2.7)-1.0));

  x           = (obj->dmajor * obj->dminor) / pow(s->theta_mb / sqrt(log(2.0)), 2.0);
  y           = 1.0 - exp(x * -1.0);
  obj->flux2  = obj->flux2 * y / x;
  obj->flux2 /= 42e9;
  obj->tmb    = obj->trj * y; 

/*
 * Shirley Tmb Calculations (CGS)
 */
  hnuk       = (6.626e-27 * frequenz / 1.381e-16);
  jnucmb     = hnuk / (exp(6.626e-27 * frequenz / (1.381e-16 * 2.725)) - 1.0);
  jnuplan    = hnuk / (exp(6.626e-27 * frequenz / (1.381e-16 * obj->tsrc)) - 1.0);
  expterm    = 1.0 - exp(-1.0 * log(2.0) * obj->dmajor*obj->dminor / pow(s->theta_mb, 2.0));
  obj->tmb   = expterm * (jnuplan - jnucmb);
  obj->csize = sqrt(log(2.0) / 2.0 * obj->dmajor * obj->dminor + pow(s->theta_mb, 2.0));

  if(logit)
  {
    st_report("JFact(Real)     = %6.2f", obj->flux2 / obj->tmb);
  }

  return(0);
}


#define NOW_DATE  1
#define USER_DATE 2
#define SCAN_DATE 3
#define ENV_DATE  4

int getPlanet()
{
  int    i, saz, rptflag = 0, user_freq = 0, wcolor = 0;
  double limit, ddate, dyr, dmon, dday;
  char   raStr[80], decStr[80], rptstr[256];
  struct tm *tmn2;
  struct EPH_OBJECT *p;
  struct EPH_STATE  *s;
  struct EPH_TIME   *t;
  static int first = 1;
  char   cmdBuf[256];

  if(first)
  {
    new_ephem.time.year = getYear(7);
    new_ephem.time.mon  = getMonth(7) -1;
    new_ephem.time.day  = getDom(7);
    new_ephem.time.hour = getHour(7);
    new_ephem.time.min  = getMinute(7);
    wcolor              = NOW_DATE;

    st_report(" ");
    st_report("Initializing Ephem Tables");

    loadPlantFluxTables(&new_ephem.time);

    st_report(" ");

    first = 0;
  }

  if(isSet(PF_R_FLAG)) /* "reload" */
  {
    new_ephem.time.year = getYear(7);
    new_ephem.time.mon  = getMonth(7) -1;
    new_ephem.time.day  = getDom(7);
    new_ephem.time.hour = getHour(7);
    new_ephem.time.min  = getMinute(7);
    wcolor              = NOW_DATE;

    loadPlantFluxTables(&new_ephem.time);

    return(0);
  }

  if(isSet(PF_N_FLAG)) /* "now" */
  {
    new_ephem.time.year = getYear(7);
    new_ephem.time.mon  = getMonth(7) -1;
    new_ephem.time.day  = getDom(7);
    new_ephem.time.hour = getHour(7);
    new_ephem.time.min  = getMinute(7);
    wcolor              = NOW_DATE;

    st_report("Using 'now' for time");
    loadPlantFluxTables(&new_ephem.time);

  }
  else
  if(isSet(PF_D_FLAG)) /* Use Environmental 'date' */
  {
    if(env->date > 0.0)
    {
      dyr = round(env->date);

      dmon = (env->date - dyr) * 100.0;
      dmon = round(dmon);

      dday = (env->date - dyr) * 100.0;
      dday -= dmon;
      dday *= 100.0;
      dday = round(dday);

      new_ephem.time.year = (int)dyr;
      new_ephem.time.mon  = (int)(dmon - 1.0);
      new_ephem.time.day  = (int)dday;

      new_ephem.time.hour = 0;
      new_ephem.time.min  = 0;
      wcolor              = ENV_DATE;

      st_report("Using ENV:Date %.4f for time", env->date);
    }
    else
    {
      new_ephem.time.year = getYear(7);
      new_ephem.time.mon  = getMonth(7) -1;
      new_ephem.time.day  = getDom(7);
      new_ephem.time.hour = getHour(7);
      new_ephem.time.min  = getMinute(7);
      wcolor              = NOW_DATE;

      st_report("Date NOT Set: Using 'now' for time");
    }

    loadPlantFluxTables(&new_ephem.time);
  }
  else
  if(isSet(PF_U_FLAG)) /* User passed 'date' */
  {
    if(args->ntokens > 1)
    {
      ddate = atof(args->secondArg);
    }
    else
    {
      ddate = 0.0;
    }

    if(ddate > 2012.01 && ddate < 2025.13) /* Sanity check */
    {
      dyr = round(ddate);

      dmon = (ddate - dyr) * 100.0;
      dmon = round(dmon);

      dday = (ddate - dyr) * 100.0;
      dday -= dmon;
      dday *= 100.0;
      dday = round(dday);

      new_ephem.time.year = (int)dyr;
      new_ephem.time.mon  = (int)(dmon - 1.0);
      new_ephem.time.day  = (int)dday;

      new_ephem.time.hour = 0;
      new_ephem.time.min  = 0;
      wcolor              = USER_DATE;

      st_report("Using User Date: %.4f for time", ddate);
    }
    else
    {
      new_ephem.time.year = getYear(7);
      new_ephem.time.mon  = getMonth(7) -1;
      new_ephem.time.day  = getDom(7);
      new_ephem.time.hour = getHour(7);
      new_ephem.time.min  = getMinute(7);
      wcolor              = NOW_DATE;

      st_report("User Date: %.4f, NOT Valid: Using 'now' for time", ddate);
      st_report("Valid User Date Range: 2012.0101 - 2025.1231");
    }

    loadPlantFluxTables(&new_ephem.time);
  }
  else
  {
    if(h0->head.scan < 1.0)
    {
      st_report("No Scan Loaded for getPlanet(); Using 'now' for time");

      new_ephem.time.year = getYear(7);
      new_ephem.time.mon  = getMonth(7) -1;
      new_ephem.time.day  = getDom(7);
      new_ephem.time.hour = getHour(7);
      new_ephem.time.min  = getMinute(7);
      wcolor              = NOW_DATE;
    }
    else
    {
      tmn2 = breakTime(h0->head.burn_time);

      new_ephem.time.year = tmn2->tm_year+1900;
      new_ephem.time.mon  = tmn2->tm_mon;
      new_ephem.time.day  = tmn2->tm_mday;
      new_ephem.time.hour = tmn2->tm_hour;
      new_ephem.time.min  = tmn2->tm_min;
      wcolor              = SCAN_DATE;
    }

    if(fluxTableOld(h0->head.burn_time))
    {
      loadPlantFluxTables(&new_ephem.time);
    }
  }

  switch(wcolor)
  {
    case USER_DATE:  setColor(RED);   break;
    default:         setColor(BLACK); break;
  }

  st_report("UT-Date:        = %02d/%02d/%04d", 
                new_ephem.time.mon+1, new_ephem.time.day, new_ephem.time.year);

  setColor(BLACK);

  if(new_ephem.time.year >= 2012)
  {
    for(i=MERCURY;i<MAX_PLANETS-1;i++)
    {
      new_ephem.planet = i;
      parseNewEphem(&new_ephem);
    }
  }

  if(env->issmt)
  {
    limit = 45.0;
  }
  else
  {
    limit = 5.0;
  }

  /* Is there a passed freq flag */
  if(isSet(PF_F_FLAG))
  {
    /* First check if the /u flag is set too; that means there are 2 passwd args */
    if(isSet(PF_U_FLAG))
    {
      if(args->ntokens > 2) /* There better be 1+2 args out there */
      {
        new_ephem.state.f_sky = atof(args->thirdArg);   /* Freq arg must be second */
        user_freq = 1;
      }
      else
      {
        st_report("Missing arg for Freq; using ENV:Freq of %f", env->freq);

        new_ephem.state.f_sky = env->freq;
      }
    }
    else
    {
      if(args->ntokens > 1) /* There better be 1+ args out there */
      {
        new_ephem.state.f_sky = atof(args->secondArg);
        user_freq = 1;
      }
      else
      {
        st_report("Missing arg for Freq; using ENV:Freq of %f", env->freq);

        new_ephem.state.f_sky = env->freq;
      }
    }
  }
  else
  {
    new_ephem.state.f_sky = env->freq;
  }

  if(user_freq)
  {
    setColor(RED);
    st_report("User Freq:      = %6.2f", new_ephem.state.f_sky);
    setColor(BLACK);
  }
  else
  {
    st_report("System Freq     = %6.2f", new_ephem.state.f_sky);
  }

  p = &new_ephem.planets[0];
  s = &new_ephem.state;
  t = &new_ephem.time;

  for(i=MERCURY;i<MAX_PLANETS-1;i++,p++)
  {
    sexagesimal(p->ra,   raStr, DDD_MM_SS_S100);
    sexagesimal(p->dec, decStr, DDD_MM_SS_S100);

    saz = 1;
    if(p->sangle < limit)
    {
      saz = 1;
    }
    else
    {
      saz = 0;
    }

    computePlanetFlux(i, p, t, s, i==MERCURY);

    bcopy((char *)p, (char *)&save_planet_values[i], sizeof(struct EPH_OBJECT));

    if(i == 0)
    {
      st_print("\n [ Tbright: A = ALMA; S = Shirley; F = Fixed ]\n\n");
      st_print("  Name       RA            Dec        Vel  Sun  SAZ  Dist  Major   Minor   Tbright   Freq  FWHM  Csize    Tmb    Flux  Flux2\n");
    }

    st_print("%8.8s %s   %s  %5.1f %5.1f  %d  %5.2f  %5.2f\"  %5.2f\" ", 
                p->name, raStr, decStr, p->vel, 
                p->sangle, saz, 
                p->dist, p->dmajor, p->dminor);

    st_print("%6.2f(%1.1s) %5.1f %5.1f\" %5.2f\" %6.2f %6.1f %6.1f\n",
                p->tsrc, p->dbase, env->freq, s->theta_mb, p->csize, p->tmb, p->flux, p->flux2);

    switch(i)
    {
      case MERCURY: env->mercury_tmb = p->tmb; break;
      case VENUS:   env->venus_tmb   = p->tmb; break;
      case MARS:    env->mars_tmb    = p->tmb; break;
      case JUPITER: env->jupiter_tmb = p->tmb; break;
      case SATURN:  env->saturn_tmb  = p->tmb; break;
      case URANUS:  env->uranus_tmb  = p->tmb; break;
      case NEPTUNE: env->neptune_tmb = p->tmb; break;
    }

    if(!strcasecmp(p->name, env->source))
    {
      env->expectedtmb  = p->tmb;
      env->expectedflux = p->flux2;
      sprintf(rptstr, "Auto Setting Expected Tmb for %s to %f", env->source, p->tmb);
      rptflag = 1;
    }
  }

  if(rptflag && env->mode == MODE_SPEC)
  {
    st_report("%s", rptstr);

    sprintf(cmdBuf, "hline /c 31 %f /Cmg", env->expectedtmb);
    st_report(cmdBuf);
    parseCmd(cmdBuf);
  }

  return(0);
}


int dopFlux(planet, flag, bfreq, efreq)
int planet, flag;
double bfreq, efreq;
{
  int i;
  double f;
  struct tm *tmn2;
  struct EPH_OBJECT *p;
  struct EPH_STATE  *s;
  struct EPH_TIME   *t;
  long tt;

  if(h0->head.scan < 1.0)
  {
    st_warn("No Scan Loaded for doFlux(); Using 'now' for time");

    new_ephem.time.year = getYear(7);
    new_ephem.time.mon  = getMonth(7) -1;
    new_ephem.time.day  = getDom(7);
    new_ephem.time.hour = getHour(7);
    new_ephem.time.min  = getMinute(7);
  }
  else
  {
    tmn2 = breakTime(h0->head.burn_time);

    new_ephem.time.year = tmn2->tm_year+1900;
    new_ephem.time.mon  = tmn2->tm_mon;
    new_ephem.time.day  = tmn2->tm_mday;
    new_ephem.time.hour = tmn2->tm_hour;
    new_ephem.time.min  = tmn2->tm_min;
  }

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

  new_ephem.planet = planet;
  parseNewEphem(&new_ephem); /* Load appropriate ephem terms */

  p = &new_ephem.planets[planet];
  s = &new_ephem.state;
  t = &new_ephem.time;

  tt = time(NULL);
  if(fluxTableOld(tt))
  {
    loadPlantFluxTables(t);
  }

  strcpy(flux_info.magic, "FLUX_INFO");
  flux_info.planet = planet;

  if(planet == VENUS)
  {
    flux_info.bfreq = venus_tbright[0].freq;
    for(i=0;i<MAX_VENUS_ENTRIES;i++)
    {
      flux_info.values[i].freq = venus_tbright[i].freq;
      flux_info.values[i].tb   = venus_tbright[i].tb;

      flux_info.efreq = flux_info.values[i].freq;
    }
    strcpy(flux_info.name, "Venus");
    flux_info.num = i;
  }
  else
  if(planet == MARS)
  {
    flux_info.bfreq = mars_tbright[0].freq;
    for(i=0;i<MAX_MARS_ENTRIES;i++)
    {
      flux_info.values[i].freq = mars_tbright[i].freq;
      flux_info.values[i].tb   = mars_tbright[i].tb;

      flux_info.efreq = flux_info.values[i].freq;
    }
    strcpy(flux_info.name, "Mars");
    flux_info.num = i;
  }
  else
  if(planet == JUPITER)
  {
    flux_info.bfreq = jupiter_tbright[0].freq;
    for(i=0;i<MAX_JUPITER_ENTRIES;i++)
    {
      flux_info.values[i].freq = jupiter_tbright[i].freq;
      flux_info.values[i].tb   = jupiter_tbright[i].tb;

      flux_info.efreq = flux_info.values[i].freq;
    }
    strcpy(flux_info.name, "Jupiter");
    flux_info.num = i;
  }
  else
  if(planet == URANUS)
  {
    flux_info.bfreq = uranus_tbright[0].freq;
    for(i=0;i<MAX_URANUS_ENTRIES;i++)
    {
      flux_info.values[i].freq = uranus_tbright[i].freq;
      flux_info.values[i].tb   = uranus_tbright[i].tb;

      flux_info.efreq = flux_info.values[i].freq;
    }
    strcpy(flux_info.name, "Uranus");
    flux_info.num = i;
  }
  else
  if(planet == NEPTUNE)
  {
    flux_info.bfreq = neptune_tbright[0].freq;
    for(i=0;i<MAX_NEPTUNE_ENTRIES;i++)
    {
      flux_info.values[i].freq = neptune_tbright[i].freq;
      flux_info.values[i].tb   = neptune_tbright[i].tb;

      flux_info.efreq = flux_info.values[i].freq;
    }
    strcpy(flux_info.name, "Neptune");
    flux_info.num = i;
  }
  else /* Only Meurcury & Saturn remain to be modeled */
  {
    f = 61.0;

    flux_info.bfreq = f;

    i = 0;
    while(f<850.0)
    {
      s->f_sky = (double)f;

      computePlanetFlux(planet, p, t, s, 0);

      flux_info.values[i].freq = f;
      flux_info.values[i].tb   = p->trj;

      flux_info.efreq = flux_info.values[i].freq;

      f += 0.2;
      i++;
    }

    strcpy(flux_info.name, eph_source_names[planet]);
    flux_info.num = i;
  }

  if(flag == 2)
  {
    if(bfreq <= 60.0)
    {
      flux_info.bfreq = 60.0;
    }
    else
    {
      flux_info.bfreq = bfreq;
    }

    if(efreq >= 850.0)
    {
      flux_info.efreq = 850.0;
    }
    else
    {
      flux_info.efreq = efreq;
    }

    flux_info.manFreq = 1;
  }

  if(xgraphic)
  {
    st_report("Sending Flux Plot to XGraphic");
    flux_info.plotwin = env->plotwin;

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

    last_plot = LAST_PLOT_FLUX_INFO;
  }

  return(0);
}


int showPlanetFlux()
{
  int flag = 0;
  double bfreq, efreq;

  if(strlen(args->secondArg) < 2)
  {
    st_fail("No Argument for showPlanetFlux()");
    return(1);
  }

  if(strlen(args->thirdArg))
  {
    bfreq = atof(args->thirdArg);
    flag++;
  }

  if(strlen(args->fourthArg))
  {
    efreq = atof(args->fourthArg);
    flag++;
  }

  if(flag == 1)
  {
    st_fail("Usage: pflux planet_name [begin_freq end_freq]");
    return(0);
  }

  if(!strncasecmp(args->secondArg, "Mercury", 7))
  {
    dopFlux(MERCURY, flag, bfreq, efreq);
  }
  else
  if(!strncasecmp(args->secondArg, "Venus", 5))
  {
    dopFlux(VENUS, flag, bfreq, efreq);
  }
  else
  if(!strncasecmp(args->secondArg, "Mars", 4))
  {
    st_report("Using Mars Flux for Date: %s", marsDate);
    dopFlux(MARS, flag, bfreq, efreq);
  }
  else
  if(!strncasecmp(args->secondArg, "Jupiter", 7))
  {
    dopFlux(JUPITER, flag, bfreq, efreq);
  }
  else
  if(!strncasecmp(args->secondArg, "Saturn", 6))
  {
    dopFlux(SATURN, flag, bfreq, efreq);
  }
  else
  if(!strncasecmp(args->secondArg, "Neptune", 7))
  {
    dopFlux(NEPTUNE, flag, bfreq, efreq);
  }
  else
  if(!strncasecmp(args->secondArg, "Uranus", 6))
  {
    dopFlux(URANUS, flag, bfreq, efreq);
  }
  else
  {
    st_fail("Unknown planet: %s", args->secondArg);
    return(1);
  }

  return(0);
}


#ifdef TOMMIE

int showPlanet()
{
  FILE *fp;
  char *string, line[256], ebuf[512], timeStr[256], *site;

  if(env->check_freq == 0 || env->check_date == 0)
  {
    st_warn("Must select freq and date first");

    return(1);
  }

  if(env->issmt)
  {
    site = "SMT";
  }
  else
  {
    site = "KITT_PEAK";
  }

  sexagesimal(env->ptime, timeStr, HH_MM_SS);

  st_report("Computing Planet Fluxes for Time/Date: %s %9.4f @ site: %s", timeStr, env->date, site);

  sprintf(ebuf, "cd %s/astro; ./astro_planets %f %s %9.4f %s", pkg_home, env->freq, timeStr, env->date, site);
//  sprintf(ebuf, "cd ../astro; ./astro_planets %f %s %9.4f %s", /*pkg_home, */env->freq, timeStr, env->date, site);

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

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

    return(1);
  }

  st_report("");

  while((string = fgets(line, sizeof(line), fp)) != NULL)
  {
    st_print("%s", line);
  }

  st_report("");

  pclose(fp);

  return(0);
}
#endif


int isPlanetVar(name)
char *name;
{
  int i;
  struct EPH_OBJECT *p;

  if(verboseLevel & VERBOSE_PARSE)
  {
    st_report("Looking if arg(%s) is planet", name);
  }

  p = &save_planet_values[0];

  if(strlen(p->name) < 2) /* Have Planets been Loaded? */
  {
    getPlanet();
  }

  for(i=0;i<MAX_PLANETS;i++,p++)
  {
    if(verboseLevel & VERBOSE_PARSE)
    {
      st_report("isPlanetVar(): Comparing %s to %s", name, p->name);
    }

    if(p->name && !strncasecmp(p->name, name, strlen(name)))
    {
      if(verboseLevel & VERBOSE_PARSE)
      {
        st_report("isPlanetVar(): Found %s", p->name);
      }

      return(0);
    }
  }

  return(1);
}


int isPlanet(name)
char *name;
{
  int i;
  struct EPH_OBJECT *p;

  if(verboseLevel & VERBOSE_PARSE)
  {
    st_report("Looking if arg(%s) is planet", name);
  }

  p = &save_planet_values[0];

  if(strlen(p->name) < 2) /* Have Planets been Loaded? */
  {
    stquiet = 1;
    getPlanet();
    stquiet = 0;
  }

  for(i=0;i<MAX_PLANETS;i++,p++)
  {
    if(verboseLevel & VERBOSE_PARSE)
    {
      st_report("isPlanetVar(): Comparing %s to %s", name, p->name);
    }

    if(p->name && !strncasecmp(p->name, name, strlen(name)))
    {
      if(verboseLevel & VERBOSE_PARSE)
      {
        st_report("isPlanetVar(): Found %s", p->name);
      }

      return(1);
    }
  }

  return(0);
}



int lookUpPlanetVal(p, val, print)
struct EPH_OBJECT *p;
char *val;
int print;
{
  void *offset;
  char tbuf[256];
  struct PLANET_LOOKUP *pl;

  for(pl=&planet_lookup[0];pl->vname;pl++)
  {
    if(!strcmp(pl->vname, val))
    {
      offset = (long)pl->offset + (char *)p;

      if(pl->type == CHARSTAR)
      {
        if(print)
        {
          sprintf(tbuf, pl->fmt, (char *)offset);
          st_print("PLANET: %s.%s = %s\n", p->name, val, tbuf);
        }
        planet_value = 0.0;     
      }
      else
      if(pl->type == DOUBLE)
      {
        if(print)
        {
          sprintf(tbuf, pl->fmt, *((double *)offset));
          st_print("PLANET: %s.%s = %s\n", p->name, val, tbuf);
        }
        planet_value = *((double *)offset);     
      }

      return(0);
    }
  }

  return(1);
}



int printPlanetVar(name)
char *name;
{
  int i;
  char *cp;
  struct EPH_OBJECT *p;

  if(verboseLevel & VERBOSE_PARSE)
  {
    st_report("printPlanetVar(): Looking up planet variable: %s", name);
  }

  p = &save_planet_values[0];

  for(i=0;i<MAX_PLANETS;i++,p++)
  {
    if(!strncasecmp(p->name, name, strlen(p->name)))
    {
      if(verboseLevel & VERBOSE_PARSE)
      {
        st_report("printPlanetVar(): Found %s", p->name);
      }

      if(strlen(name) == strlen(p->name)) /* No .arg so print whole thing */
      {
        st_report("Name:   %s",      p->name);
        st_report("Dbase:  %s",      p->dbase);
        st_report("RA:     %15.10f", p->ra);
        st_report("Dec:    %15.10f", p->dec);
        st_report("Illum:  %15.10f", p->illum);
        st_report("Dmajor: %15.10f", p->dmajor);
        st_report("Dminor: %15.10f", p->dminor);
        st_report("Sdist:  %15.10f", p->sdist);
        st_report("Rdot:   %15.10f", p->rdot);
        st_report("Dist:   %15.10f", p->dist);
        st_report("Vel:    %15.10f", p->vel);
        st_report("Sangle: %15.10f", p->sangle);
        st_report("Tsrc:   %15.10f", p->tsrc);
        st_report("TRJ:    %15.10f", p->trj);
        st_report("Flux:   %15.10f", p->flux);
        st_report("Flux2:  %15.10f", p->flux2);
        st_report("TMB:    %15.10f", p->tmb);
        st_report("Csize:  %15.10f", p->csize);

        return(0);
      }
      else /* Print element of struct */
      {
        cp = strstr(name, ".");

        if(cp)
        {
          cp++;

          if(!lookUpPlanetVal(p, cp, 1))
          {
            return(0);
          }
        }
        else
        {
          return(1);
        }
      }
    }
  }

  return(1);
}


/* Return True if Variable Found */
int findPlanetVal(name, dret)
char *name;
double *dret;
{
  int ret, oldqt;

  oldqt = stquiet;


  ret = isPlanetVar(name);
  if(!ret)
  {
    stquiet = 1;
    if(!printPlanetVar(name))
    {
      *dret = planet_value;
      stquiet = oldqt;

      return(1);
    }
  }

  stquiet = oldqt;

  return(0);
}