atm.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 SMT_ELEVATION (3.186)
#define KP_ELEVATION  (1.914)

#define D2R (M_PI / 180.0)

struct ALMA_BANDS {
        char  *name;
        char  *color;
        char  *frontback;
        double lower;
        double upper;
};

struct ALMA_BANDS alma_bands[] = {
        { " Band-1",  "light-coral",    "behind",  ALMA_BAND_1_LOWER,  ALMA_BAND_1_UPPER },
        { " Band-2",  "orange",         "behind",  ALMA_BAND_2_LOWER,  ALMA_BAND_2_UPPER },
        { " Band-3",  "yellow",         "behind",  ALMA_BAND_3_LOWER,  ALMA_BAND_3_UPPER },
        { " Band-4",  "light-green",    "behind",  ALMA_BAND_4_LOWER,  ALMA_BAND_4_UPPER },
        { " Band-5",  "light-blue",     "behind",  ALMA_BAND_5_LOWER,  ALMA_BAND_5_UPPER },
        { " Band-6",  "aquamarine",     "behind",  ALMA_BAND_6_LOWER,  ALMA_BAND_6_UPPER },
        { " Band-7",  "violet",         "behind",  ALMA_BAND_7_LOWER,  ALMA_BAND_7_UPPER },
        { " Band-8",  "light-cyan",     "behind",  ALMA_BAND_8_LOWER,  ALMA_BAND_8_UPPER },
        { " Band-9",  "light-magenta",  "behind",  ALMA_BAND_9_LOWER,  ALMA_BAND_9_UPPER },
        { "Band-10",  "coral",          "behind", ALMA_BAND_10_LOWER, ALMA_BAND_10_UPPER }
};

#define MAX_ALMA_BANDS (sizeof(alma_bands) / sizeof(struct ALMA_BANDS))

#define MAX_ATM_ELEMENTS 10000

float xdata[MAX_ATM_ELEMENTS];
float ydata[MAX_ATM_ELEMENTS];
float tdata[MAX_ATM_ELEMENTS];
float odata[MAX_ATM_ELEMENTS];
float bdata[MAX_ATM_ELEMENTS];

int    nAtm;
double mmh2oAtm, amAtm, tambAtm, pressAtm, elevAtm, fminAtm, fmaxAtm, fdeltaAtm;
char   siteAtm[256];

int atm_atmosp_(float *tamb, float *press, float *h0);
int atm_transm_(float *mmwat, float *am, float *freqs, float *tems, float *tatms, float *tauo, float *tauwat, float *tautots, int *ier);


/* Use more accurate calculations; for now */
int loadAtmFromServer(which, fmin, fmax, fstep, pwv, gtemp, gpress, altitude, ymin)        
char *which;
double fmin, fmax, fstep, pwv, gtemp, gpress, altitude, ymin;
{
  int indx = 0, n;
  char pbuf[1024], bbuf[256], vbuf[256], ebuf[256];
  char *string, line[256];
  FILE *fp;
  double freq, dry, wet, sky, trans;

  sprintf(bbuf, "ssh %s \"/usr/local/bin/absorption ", which);

  sprintf(vbuf, "--fmin %.2f --fmax %.2f --fstep %.2f --pwv %.2f --gtemp %.1f --gpress %.0f --altitude %.0f", 
                        fmin, fmax, fstep, pwv, gtemp, gpress, altitude);

  strcpy(ebuf, " \"");

  sprintf(pbuf, "%s %s %s", bbuf, vbuf, ebuf);

  if(isSet(PF_V_FLAG))
  {
    printf("%s\n", pbuf);
  }

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

    return(1);
  }

  st_report(" ");
  st_report("Using Server: %s ", which);
  st_report("Standby, as this may take a few seconds");
  st_report(" ");

  while((string = fgets(line, sizeof(line), fp)) != NULL)
  {
    n = sscanf(line, "%le,\t%lf,\t%lf,\t%lf", &freq, &dry, &wet, &sky);

    if(n == 4)
    {
      xdata[indx] = freq / 1e9;

      tdata[indx] = wet;
      odata[indx] = dry;
      bdata[indx] = wet+dry;

      trans = 1.0 - (wet+dry);

      if(trans < ymin || gsl_isnan(trans))
      {
        trans = ymin;
      }
     
      ydata[indx] = trans;

      indx++;

      if(indx >= MAX_ATM_ELEMENTS)
      {
        break;
      }
    }
  }

  pclose(fp);

  return(indx);
}


/* airmass calculations */
float computeAtmophers(elv)
float elv;
{
  float am = 1.0;
// only use doubles in the calc of am
  double airmass, hz, hz1, hz2, hz3, rzero, hzero, eps, Gamma;

  airmass = 1.0 / sin(elv * D2R);

  rzero = 6372.0;
  hzero = 7.8;

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("Telescope Elev: %.1f", elv);
  }

  eps = asin(rzero / (rzero + hzero) * cos(elv * D2R));

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("EPS = %f", eps);
  }

  Gamma = elv*D2R +eps;

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("Gamma = %f", Gamma);
  }

  hz1 = pow(rzero, 2.0);

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("hz1 = %f", hz1);
  }

  hz2 = pow(rzero+hzero, 2.0);

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("hz2 = %f", hz2);
  }

  hz3 = 2.0 * rzero * (rzero+hzero) * sin(Gamma);

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("hz3 = %f", hz3);
  }

  hz = hz1 + hz2 - hz3;

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("HZ double sums = %f", hz);
  }

  hz = sqrt(hz);

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("sqrt(HZ) = %f", hz);
  }

  am = (float)(hz/hzero);

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("AM = %.3f; (1.0 / sin(%.1f)) = %.3f", am, elv, airmass);
  }

  return(am);
}



int findMinMax(fmin, fmax)
double *fmin, *fmax;
{
  int set=0;
  double fn=0.0, fx=0.0;

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("findMinMax(): Enter %f %f", *fmin, *fmax);
  }

  if(isSet(PF_1_FLAG))
  {
    fn = ALMA_BAND_1_LOWER;
    fx = ALMA_BAND_1_UPPER;
    set = 1;
  }

  if(isSet(PF_2_FLAG))
  {
    if(fn == 0.0)
    {
      fn = ALMA_BAND_2_LOWER;
    }

    fx = ALMA_BAND_2_UPPER;
    set = 1;
  }

  if(isSet(PF_3_FLAG))
  {
    if(fn == 0.0)
    {
      fn = ALMA_BAND_3_LOWER;
    }

    fx = ALMA_BAND_3_UPPER;
    set = 1;
  }

  if(isSet(PF_4_FLAG))
  {
    if(fn == 0.0)
    {
      fn = ALMA_BAND_4_LOWER;
    }

    fx = ALMA_BAND_4_UPPER;
    set = 1;
  }

  if(isSet(PF_5_FLAG))
  {
    if(fn == 0.0)
    {
      fn = ALMA_BAND_5_LOWER;
    }

    fx = ALMA_BAND_5_UPPER;
    set = 1;
  }

  if(isSet(PF_6_FLAG))
  {
    if(fn == 0.0)
    {
      fn = ALMA_BAND_6_LOWER;
    }

    fx = ALMA_BAND_6_UPPER;
    set = 1;
  }

  if(isSet(PF_7_FLAG))
  {
    if(fn == 0.0)
    {
      fn = ALMA_BAND_7_LOWER;
    }

    fx = ALMA_BAND_7_UPPER;
    set = 1;
  }

  if(isSet(PF_8_FLAG))
  {
    if(fn == 0.0)
    {
      fn = ALMA_BAND_8_LOWER;
    }

    fx = ALMA_BAND_8_UPPER;
    set = 1;
  }

  if(isSet(PF_9_FLAG))
  {
    if(fn == 0.0)
    {
      fn = ALMA_BAND_9_LOWER;
    }

    fx = ALMA_BAND_9_UPPER;
    set = 1;
  }

  if(isSet(PF_10_FLAG))
  {
    if(fn == 0.0)
    {
      fn = ALMA_BAND_10_LOWER;
    }

    fx = ALMA_BAND_10_UPPER;
    set = 1;
  }

  if(set)
  {
    *fmin = fn;
    *fmax = fx;

    if(verboseLevel & VERBOSE_ATM)
    {
      st_report("findMinMax(): Leaving %f %f", *fmin, *fmax);
    }
  }

  return(0);
}


int generateAtm(fmin, fmax, mmwat, tamb, press, elev)
double fmin, fmax;
float   mmwat, tamb, press, elev;
{
  int    indx=0, ier;
  double delta;
  float  am=1.0, freqs;               /* Inputs */
  float  tems, tatms, tauo, tauwat, tautots; /* Returns */
  float  trans;

  bzero((char *)&xdata, sizeof(xdata));
  bzero((char *)&ydata, sizeof(ydata));
  bzero((char *)&tdata, sizeof(tdata));
  bzero((char *)&odata, sizeof(odata));
  bzero((char *)&bdata, sizeof(bdata));

  if(mmwat < 0.01)
  {
    st_fail("mm H2O must be greater than 0.01");
    return(0);
  }

  st_report("Computing atm profile: f = %.1f to %.1f", fmin, fmax);

  if(!strcmp(site_name, "TUCSON") && notSet(PF_F_FLAG)) /* Use ATM server unless fast flag passed */
  {
    indx = 0;

    delta = (fmax - fmin) / 8192.0;
    indx = loadAtmFromServer("localhost", fmin, fmax, delta, mmwat, tamb, press, elev*1e3, 0.0);

    if(indx == 0)
    {
      st_fail("No ATM elements found");
      return(1);
    }

    st_report("Loaded %d ATM Elements", indx);
  }
  else
  if(!strcmp(site_name, "KITT_PEAK") && notSet(PF_F_FLAG)) /* Use ATM server unless fast flag passed */
  {
    indx = 0;

    delta = (fmax - fmin) / 8192.0;
    indx = loadAtmFromServer("modelo", fmin, fmax, delta, mmwat, tamb, press, elev*1e3, 0.0);

    if(indx == 0)
    {
      st_fail("No ATM elements found");
      return(1);
    }

    st_report("Loaded %d ATM Elements", indx);
  }
  else
  if(!strcmp(site_name, "MT_GRAHAM") && notSet(PF_F_FLAG)) /* Use ATM server unless fast flag passed */
  {
    indx = 0;

    delta = (fmax - fmin) / 8192.0;
    indx = loadAtmFromServer("smtoast", fmin, fmax, delta, mmwat, tamb, press, elev*1e3, 0.0);

    if(indx == 0)
    {
      st_fail("No ATM elements found");
      return(1);
    }

    st_report("Loaded %d ATM Elements", indx);
  }
  else                                  /* Use internal calculations */
  {
    atm_atmosp_(&tamb, &press, &elev); // One call keeps model in memory

    indx  = 0;
    freqs = fmin;
    delta = (fmax - fmin) / 8192.0;
    while(freqs < fmax)
    {
      atm_transm_(&mmwat, &am, &freqs, &tems, &tatms, &tauo, &tauwat, &tautots, &ier);

      trans = 1.0 - tautots;

      if(trans < 0.0 || gsl_isnan(trans))
      {
        trans = 0.0;
      }

      xdata[indx] = freqs;
      tdata[indx] = tauwat;
      odata[indx] = tauo;
      bdata[indx] = tautots;
      ydata[indx] = trans;

      indx++;

      if(indx >= MAX_ATM_ELEMENTS)
      {
        break;
      }

      freqs += delta;
    }

    st_report("Loaded %d ATM Elements", indx);
  }

  nAtm      = indx;
  mmh2oAtm  = mmwat;
  amAtm     = am;
  tambAtm   = tamb;
  pressAtm  = press;
  elevAtm   = elev;
  fminAtm   = fmin;
  fmaxAtm   = fmax;
  fdeltaAtm = delta;

  return(0);
}


int doATM()
{
  int    i, j, ret, abox, yscale = 0;;
  char   sysbuf[1024], demname[256], datname[256];
  char   *cp, *site, plotfname[256], tbuf[256], pdfname[512];
  double min = 1e9, max = -1e9, freq, fmin, fmax, ymin, ymax, 
         freqDiffMin, freqdiff, freqVal, wl, x1;

  double ofreq, ofreqDiffMin, ofreqdiff, ofreqVal, bw;

  float  tamb, press, elev;
  float  mmwat=0.0;               /* Inputs */
  float  *f;
  FILE   *fp;
  struct ALMA_BANDS *p;

  if(isSet(PF_L_FLAG))                  /* List what we have used to generate the model */
  {
    if(nAtm)
    {
      st_report(" ");
      st_report("Number of Elements in Model:   %6d", nAtm);
      st_report("Starting Freq of Model:       %7.1f GHz", fminAtm);
      st_report("Ending Freq of Model:         %7.1f GHz", fmaxAtm);
      st_report("Delta Freq of Model:            %5.3f GHz", fdeltaAtm);
      st_report("Water Vapor Content:            %.3f mm", mmh2oAtm);
      st_report("Number of Atmophers:            %.3f", amAtm);
      st_report("Surface Temperature:          %.3f K", tambAtm);
      st_report("Surface Barometric Pressure:  %.3f mbar", pressAtm);
      st_report("Site:                %16.16s", siteAtm);
      st_report("Site Elevation:                 %.3f km", elevAtm);
      st_report(" ");
    }
    else
    {
      st_report("ATM Model not generated yet");
    }

    return(0);
  }

  if(env->issmt)
  {
    elev = SMT_ELEVATION;
    site = "Mt Graham";
  }
  else
  {
    elev = KP_ELEVATION;
    site = "Kitt Peak";
  }

  if(nAtm == 0 || isSet(PF_G_FLAG))             /* Alway generate the full 10 bands of the model */
  {

    if(isSet(PF_A_FLAG))                                /* Fetch from scan header ?? */
    {
      if(h0->head.scan == 0.0)
      {
        st_fail("Must load scan first before using Auto Mode");
        return(1);
      }

      mmwat = h0->head.mmh2o;
      tamb  = h0->head.tamb + 273.4;
      press = h0->head.pressure * 1.33222;      /* Convert to millibars from torr */
    }
    else
    if(isSet(PF_K_FLAG))                                /* Test using Mauma Kea values ?? */
    {
      if(args->ntokens < 2)                     /* If you are not fetching from header; Look for passed arg */
      {
        printH_FlagHelp(args->firstArg);

        return(1);
      }

      mmwat = atof(args->secondArg);
      tamb  = 270.0;
      press = 560.0; 
      elev  = 5.000;
      site  = "Mauna Kea";
    }
    else
    if( !nAtm && args->ntokens < 2)             /* If you are not fetching from header; Look for passed arg */
    {
      printH_FlagHelp(args->firstArg);

      return(1);
    }
    else                                                /* Use some defaults */
    {
      mmwat = atof(args->secondArg);

      if(env->issmt)
      {
        press = 696.0; 
        tamb  = 0.0 + 273.4;                    /* ATP Set to some defaults */
      }
      else
      {
        press = 810.0;                          /* millibars not torr */
        tamb  = 15.0 + 273.4;
      }
    }

    if(mmwat && mmwat != mmh2oAtm)              /* User passed a different water vapor amount; force regen of model */
    {
      nAtm = 0;
    }

    if(strcmp(siteAtm, site))                   /* User passwd a different site; force regen of model */
    {
      nAtm = 0;
    }

    fmin = ALMA_BAND_1_LOWER;
    fmax = ALMA_BAND_10_UPPER;

    generateAtm(fmin, fmax, mmwat, tamb, press, elev);

    strcpy(siteAtm, site);                      /* Remember fro which site this model was generated for */

    setFlag(PF_L_FLAG, 1);                      /* Force the 'list' flag and call myself */

    doATM();

    return(0);
  }
                                                /* User may only want to display part of the model */
  if(args->ntokens > 2)                         /* Possible range of args  */
  {                                             /* 1              2    3   */
    if(args->ntokens == 3)                      /* atm [options] fmin fmax */
    {
      fmin = atof(args->secondArg);
      fmax = atof(args->thirdArg);
    }
    else                                        /* 1              2    3    4    5 */
    if(args->ntokens == 5)                      /* atm [options] fmin fmax ymin ymax */
    {
      fmin = atof(args->secondArg);
      fmax = atof(args->thirdArg);

      ymin = atof(args->fourthArg);
      ymax = atof(args->fifthArg);

      if(ymin > 0.0 || ymax > 0.0)
      {
        yscale = 1;
      }
    }

    if(fmin > fmax || fmin < 1.0 || fmin > 1000.0 || fmax < 0.0 || fmax > 1000.0)
    {
      st_fail("fmin: %f > fmax: %f or out of range [1.0 - 1000.0]", ymin, ymax);

      printH_FlagHelp(args->firstArg);

      return(1);
    }

    if(ymin > ymax)
    {
      st_fail("ymin: %f > ymax: %f", ymin, ymax);

      printH_FlagHelp(args->firstArg);

      return(1);
    }
  }
  else                                          /* Look for range limits in flags /1, /2, /3 ... */
  {
    fmin =   35.0;
    fmax = 1000.0;

    findMinMax(&fmin, &fmax);
  }

  /* Now plot it out */
  if(isSet(PF_W_FLAG))                          /* water contribution ?? */
  {
    f = &tdata[0];
  }
  else
  if(isSet(PF_O_FLAG))                          /* Oxygen contribution ?? */
  {
    f = &odata[0];
  }
  else
  if(isSet(PF_B_FLAG))                          /* Both ?? */
  {
    f = &bdata[0];
  }
  else
  {
    f = &ydata[0];                              /* Else, transmission of atmosphere */
  }
                                                /* Go write the dat file */
  sprintf(datname, "/tmp/%s/atm.dat", session);

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("Opening ATM Dat File: %s", datname);
  }

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

    return(1);
  }

  if(h0->head.bw)
  {
    bw = h0->head.bw / 1e3;
  }
  else
  {
    if(env->issmt && env->mode == MODE_SPEC)
    { 
      bw = 1.0;
    }
    else
    {
      bw = 4.0;
    }
  }

  freq = h0->head.obsfreq / 1e3;
   freqDiffMin = 1e9;
  ofreqDiffMin = 1e9;

  /* determine other sideband freq */
  if(h0->head.sideband == 2)
  {
    ofreq = (h0->head.obsfreq / 1e3) - ((h0->head.firstif / 1e3) * 2.0);
  }
  else
  {
    ofreq = (h0->head.obsfreq / 1e3) + ((h0->head.firstif / 1e3) * 2.0);
  }

  for(i=0;i<nAtm;i++,f++)
  {
    fprintf(fp, "%f %f\n", xdata[i], *f*100.0);

    if(*f*100.0 < min)
    {
      min = *f*100.0;
    }

    if(*f*100.0 > max)
    {
      max = *f*100.0;
    }

    freqdiff = fabs(xdata[i] - freq);           /* Find the atm transmission at observed freq */
    if(freqdiff < freqDiffMin)
    {
      freqDiffMin = freqdiff;

      if(isSet(PF_W_FLAG))
      {
        freqVal = tdata[i];
      }
      else
      if(isSet(PF_O_FLAG))
      {
        freqVal = odata[i];
      }
      else
      if(isSet(PF_B_FLAG))
      {
        freqVal = bdata[i];
      }
      else
      {
        freqVal = ydata[i];
      }
    }

    ofreqdiff = fabs(xdata[i] - ofreq);           /* Find the atm transmission at observed freq */
    if(ofreqdiff < ofreqDiffMin)
    {
      ofreqDiffMin = ofreqdiff;

      if(isSet(PF_W_FLAG))
      {
        ofreqVal = tdata[i];
      }
      else
      if(isSet(PF_O_FLAG))
      {
        ofreqVal = odata[i];
      }
      else
      if(isSet(PF_B_FLAG))
      {
        ofreqVal = bdata[i];
      }
      else
      {
        ofreqVal = ydata[i];
      }
    }
  }

  fclose(fp);
                                                /* Go write the dem file */
  sprintf(demname, "/tmp/%s/atm.dem", session);

  if(verboseLevel & VERBOSE_ATM)
  {
    st_report("Opening ATM Dem File: %s", demname);
  }

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

    return(1);
  }

  fprintf(fp, "reset\n");

  if(isSet(PF_W_FLAG))
  {
    fprintf(fp, "set title 'Tau h2o @ %.2f mm H2O' tc rgb \"red\"\n", mmh2oAtm);
    fprintf(fp, "set ylabel 'Atmospheric Tau (Water)' tc rgb \"red\"\n");
  }
  else
  if(isSet(PF_O_FLAG))
  {
    fprintf(fp, "set title 'Tau O2 @ %.2f mm H2O' tc rgb \"red\"\n", mmh2oAtm);
    fprintf(fp, "set ylabel 'Atmospheric Tau (O2)' tc rgb \"red\"\n");
  }
  else
  if(isSet(PF_B_FLAG))
  {
    fprintf(fp, "set title 'Tau O2+Water @ %.2f mm H2O' tc rgb \"red\"\n", mmh2oAtm);
    fprintf(fp, "set ylabel 'Atmospheric Tau (O2+Water)' tc rgb \"red\"\n");
  }
  else
  {
    fprintf(fp, "set title 'Atmospheric Transmission: %s                          PWV: %.2f mm' tc rgb \"red\"\n", siteAtm, mmh2oAtm);
    fprintf(fp, "set ylabel 'Transmission (%%)' tc rgb \"red\"\n");
  }

  fprintf(fp, "set xlabel 'Frequency (GHz)' tc rgb \"red\"\n");

  fprintf(fp, "set xrange [%.2f:%.2f]\n", fmin, fmax);

  if(yscale)
  {
    fprintf(fp, "set yrange [%.2f:%.2f]\n", ymin*100.0, ymax*100.0);
  }
  else
  if(notSet(PF_B_FLAG) && notSet(PF_O_FLAG) && notSet(PF_W_FLAG))
  {
    fprintf(fp, "set yrange [-5.0:100.0]\n");  /* Give a little space below 0.0 for visibility */
  }
  /* Otherwise, let gnuplot figure it out */

  fprintf(fp, "set grid\n");
  fprintf(fp, "set key off\n");

  if(isSet(PF_U_FLAG) && freq)
  {
    st_report("Atm Transparency @ %f: %f", freq, freqVal);

    fprintf(fp, "set object 8 rect front from %f,%f to  %f,%f fill transparent pattern 2 fc rgb 'green'\n", 
                        freq - (bw/2.0), (freqVal-0.02)*100.0, freq + (bw/2.0), (freqVal+0.02)*100.0);
  }

  if(isSet(PF_U_FLAG) && ofreq)
  {
    st_report("Atm Transparency @ %f: %f", ofreq, ofreqVal);

    fprintf(fp, "set object 9 rect front from %f,%f to  %f,%f fill transparent pattern 2 fc rgb 'magenta'\n", 
                        ofreq - (bw/2.0), (ofreqVal-0.02)*100.0, ofreq + (bw/2.0), (ofreqVal+0.02)*100.0);
  }

  if(isSet(PF_M_FLAG))
  {
    p = &alma_bands[0];
    for(j=0;j<MAX_ALMA_BANDS;j++,p++)
    {
      abox = 0;

      if(p->lower >= fmin && p->upper <= fmax) /* Fully visible */
      {
        fprintf(fp, "set object %d rect %s from %f,%f to  %f,%f fill transparent solid 0.6 fc rgb '%s'\n", 
                                j+10, p->frontback, p->lower, min, p->upper, max, p->color);

        abox = 1;
      }
      else                                      /* partially visiable */
      if(p->lower >= fmin)
      {
        fprintf(fp, "set object %d rect %s from %f,%f to  %f,%f fill transparent solid 0.6 fc rgb '%s'\n", 
                                j+10, p->frontback, p->lower, min, fmax, max, p->color);

        abox = 1;
      }
      else
      if(p->upper <= fmax)
      {
        fprintf(fp, "set object %d rect %s from %f,%f to  %f,%f fill transparent solid 0.6 fc rgb '%s'\n", 
                                j+10, p->frontback, fmin, min, p->upper, max, p->color);

        abox = 1;
      }

      if(abox)
      {
        x1 = p->lower + (p->upper-p->lower) / 2.0;
        if(x1 > fmin && x1 < fmax) /* Trap labels outside the plotting area */
        {
          wl = 300e3 / (( ((p->upper-p->lower) / 2.0) + p->lower) * 1e3);
          sprintf(tbuf, "%s %3.1fmm", p->name, wl);
          fprintf(fp, "set label '%s' at %.6f, %.6f rotate by 90 tc rgb 'black' \n", tbuf, x1, (max-min) / 2.0);
        }
      }
    }
  }

  fprintf(fp, "#\n");

  if(isSet(PF_P_FLAG) || isSet(PF_S_FLAG))      /* print or save ?? */
  {
    fprintf(fp, "set terminal postscript color solid enhanced landscape\n");

    sprintf(plotfname, "/tmp/%s/ccal_plot.ps", session);

    fprintf(fp, "set out '%s'\n", plotfname);
  }
  else
  {
    fprintf(fp, "set terminal x11 enhanced %d\n", env->plotwin+1 );
  }

  fprintf(fp, "plot '%s' with lines lc rgb \"red\"", datname);

  fclose(fp);

  sprintf(sysbuf, "xgraphic load '/tmp/%s/atm.dem'", session);

  parseCmd(sysbuf);

  if(isSet(PF_P_FLAG))
  {
    usleep(1000000); /* Give Xgraphics time to plot */
    sprintf(sysbuf, "lpr /tmp/%s/ccal_plot.ps", session);
    system(sysbuf);
  }
  else
  if(isSet(PF_S_FLAG))
  {
    usleep(1000000); /* Give Xgraphics time to plot */
    sprintf(pdfname, "atm-%s-%.2fmm.pdf", siteAtm, mmh2oAtm);

    cp = strchr(pdfname, ' ');
    if(cp)
    {
      *cp = '_';
    }
    sprintf(sysbuf, "cat /tmp/%s/ccal_plot.ps | /usr/bin/epstopdf -f --outfile=%s", session, pdfname);
    st_report("%s", sysbuf);

    ret = system(sysbuf);

    if(!ret)
    {
      st_report("AMT Plot saved to %s", pdfname);
    }
    else
    {
      st_fail("Error creating PDF file");
    }
  }

  return(0);
}