Initial version to code.arm.gov. Does not include the latest development for...

Initial version to code.arm.gov. Does not include the latest development for mfrsr head id that Annette is working on.

src/Alexandrov_cloud_screen.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+
+/***************************************************/
+int isleapyear( int y )
+{
+  /* Check. */
+  if ( ( ( y % 4 == 0 ) && ( y % 100 != 0 ) ) ||
+       ( ( y % 100 == 0 ) && ( y % 400 == 0 ) ) ) 
+    {
+      /* Is leap year. */
+      return 1;
+    }
+  /* Not leap year. */
+  return 0;
+}
+
+/***************************************************
+
+
+/***************************************************
+ * function: epoch2julian
+ *
+ * DESCRIPTION:
+ *   This procedure will convert epoch time into
+ *    julian time.
+ *
+ *   The epoch time is seconds since 00:00:00 UTC
+ *   January 1, 1970.
+ *
+ *   The julian time is the number of days since
+ *   00:00:00 UTC  January 1 of the appropriate
+ *   year.
+ *
+ * INPUT VARIABLES:
+ *     time_t t - Seconds since 00:00:00 UTC January 1, 1970.
+ *
+ *
+ * OUTPUT VARIABLES:
+ *     double *jt - The number of days since 00:00:00 UTC
+ *                  January 1, for the appropriate year.
+ *
+ * RETURN VARIABLE:
+ *    int -
+ *        EXIT_SUCCESS - If the procedure completed successfully.
+ *        EXIT_FAILURE - If there was an error in the procedure.
+ *
+ *
+ **************************************************
+ */
+int epoch2julian( time_t t, double *jt )
+{
+  struct tm *gmt, jan1;
+  time_t     jan1_secs;
+  int        diff_secs;
+  time_t     timezone;
+  time_t     tt;
+
+  /* Error check. */
+  if ( jt == ( double * )NULL ) 
+    {
+      /* Return failure. */
+      return EXIT_FAILURE;
+    }
+
+  /* Initialize. */
+  tt = 0;
+
+  /* Determine the time offset for the current time zone. */
+  timezone = mktime( gmtime( &tt ) );
+
+  /* Convert seconds since 1970 to something useful. */
+  gmt = gmtime( &t );
+
+  /*
+   * Populate the struct with info
+   * for January 1st of current year.
+   */
+  jan1.tm_year = gmt->tm_year;
+  jan1.tm_mon = 0;
+  jan1.tm_mday = 1;
+  jan1.tm_hour = 0;
+  jan1.tm_min = 0;
+  jan1.tm_sec = 0;
+  jan1.tm_isdst = -1;
+
+  /* Convert from gmtime to seconds since January 1, 1970. */
+  if ( ( jan1_secs = mktime( &jan1 ) ) == -1 ) 
+    {
+      /* Print an error message. */
+      printf( "\n\nError converting gmtime to "
+	      "seconds since 1970 in %s line %d.\n\n",
+	      __FILE__, __LINE__ );
+
+      /* Return failure. */
+      return EXIT_FAILURE;
+    }
+
+  /* Adjust the timezone accordingly. */
+  jan1_secs -= timezone;
+
+  /* Seconds since January 1 of the current year. */
+  diff_secs = t - jan1_secs;
+
+  /* Get the julian time since January 1 of the current year. */
+  *jt = diff_secs / 86400.0;
+
+  /* Return success. */
+  return EXIT_SUCCESS;
+}
+
+/***************************************************
+
+function: epoch2serial
+
+DESCRIPTION:
+    This procedure will convert epoch time into
+    serial time.
+
+    The epoch time is seconds since 00:00:00 UTC
+    January 1, 1970.
+
+    The serial time is the number of days since
+    00 AD.
+
+INPUT VARIABLES:
+       time_t t - Seconds since 00:00:00 UTC January 1, 1970.
+
+OUTPUT VARIABLES:
+      double *s - The number of days since 00 AD.
+
+RETURN VARIABLE:
+    int -
+        EXIT_SUCCESS - If the procedure completed successfully.
+        EXIT_FAILURE - If there was an error in the procedure.
+
+SEE ALSO:
+
+***************************************************/
+int epoch2serial( time_t t, 
+		  double *s )
+{
+  long       y;
+  double     ydays, jdays;
+  struct tm *epoch;
+
+  /* Error check. */
+  if ( s == ( double * )NULL ) 
+    {
+      /* Return failure. */
+      return EXIT_FAILURE;
+    }
+
+  /* Convert seconds since 1970 to time struct. */
+  epoch = gmtime( &t );
+
+  /* Count days in whole years from 00 AD. */
+  ydays = ( epoch->tm_year + 1900 ) * 365.0;
+  for ( y = 0; y < ( epoch->tm_year + 1900 ); y += 4 ) 
+    {
+      if ( isleapyear( y ) ) 
+	{
+	  ydays += 1.0;
+	}
+    }
+
+  /* Find julian day. */
+  if ( epoch2julian( t, &jdays ) == EXIT_FAILURE ) 
+    {
+      /* Return failure. */
+      return EXIT_FAILURE;
+    }
+
+  /* Add julian day to year days from 00 AD.
+   *
+   * NOTE - 1 is added for Matlab serial time format
+   *        compatibility.  Jan 1 = Day 1 vs. Day 0.
+   */
+  *s = ( double )( jdays + ydays + 1.0 );
+
+  /* Return success. */
+  return EXIT_SUCCESS;
+}
+
+void AlexandrovCloudScreen(int    nsamples,          /* input, sample count */
+			   double time[],            /* input, sample times */
+			   double tau[],             /* input, total OD, filter 5 */
+			   double eps[],             /* output, epsilon (value used to 
+						      * determine if cloud is present
+						      */
+			   int    cloud_detected[])  /* output: flag indicating 
+						      * cloud detection for each sample 
+						      */
+{
+  int    i,j,c,p,t,count;
+    
+  int    bar[4318],
+         pos_tau[4318],
+         pos_tau_prime[4318];
+  double abstime,meantau,meanlogtau;
+  double tau_bar[4318],
+         tau_prime[4318],
+         tau_prime_bar[4318],
+         bar_log_tau_prime[4318];
+
+  double gAOT       = 0.2;
+  double gThreshold = 0.0001;
+
+  /*    aod = gAOT;  value is 0.2 (default) */
+
+  /* 
+   * =======================================================
+   * PASS 1: Compute tau prime, a renormalized tau.
+   *         
+   *         Note: examining filter 5 data only
+   *
+   * =======================================================
+   */
+
+  c                  = 0;
+  for (i=0; i<nsamples; i++) 
+    {
+      tau_bar[i]     = 0.0;
+      tau_prime[i]   = 0.0;
+      if (tau[i] > 0.0 && tau[i] < 2.0) 
+	{
+	  pos_tau[c++] = i;
+	}
+    }
+
+  count              = c;
+
+  for (i=0; i<count; i++) 
+    {
+      t              = pos_tau[i];
+
+      c=0;
+      for (j=0; j<count; j++) 
+	{
+	  abstime    = time[t] - time[pos_tau[j]];
+	  if (abstime < 0.0)
+	    {
+	      abstime = -1.0 * abstime;
+	    }
+	  if (abstime <= (double)(5.0/(24.0*60.0))) 
+	    {
+	      bar[c++] = pos_tau[j];
+	    }
+	}
+
+      if (c > 5) 
+	{
+	  meantau      = 0.0;
+	  for (j=0; j<c; j++) 
+	    {
+	      p        = bar[j];
+	      meantau  = meantau + tau[p];
+	    }
+	  meantau      = meantau /(double)c;
+
+	  tau_bar[t]   = meantau;
+	  tau_prime[t] = (tau[t] - tau_bar[t]) + gAOT;
+        }
+
+    } /* End PASS 1 */
+
+  /*
+   * =========================================================
+   * PASS 2: Compute tau_prime_bar, the average of tau_prime.
+   * =========================================================
+   */
+
+  c                     = 0;
+  for (i=0; i<nsamples; i++) 
+    {
+      tau_prime_bar[i]  = 1.0;
+
+      cloud_detected[i] = 0;       /* 0 = no cloud detected, 
+				    * 1 = cloud has been detected 
+				    */
+
+      eps[i]            = 1.0;     
+				   
+      if (tau_prime[i] > 0.0) 
+	{
+	  pos_tau_prime[c++] = i;   /* index to a positive tau_prime */
+	}
+    }
+  count = c;                        /* number of positive tau primes */
+
+  for (i=0; i<count; i++)
+    {
+      t               = pos_tau_prime[i];
+      
+      c               = 0;
+      for (j=0; j<count; j++) 
+	{
+	  abstime     = time[t] - time[pos_tau_prime[j]];
+	  if (abstime < 0.0) 
+	    {
+	      abstime = -1.0 * abstime;
+	    }
+	  if (abstime <= (5.0/(24.0*60.0))) 
+	    {
+	      bar[c++]  = pos_tau_prime[j];
+	    }
+	}
+
+      if (c > 5) 
+	{
+	  meantau                    = 0.0;
+	  meanlogtau                 = 0.0;
+	  for (j=0; j<c; j++) 
+	    {
+	      p                      = bar[j];
+	      meantau                = meantau + tau_prime[p];
+	      meanlogtau             = meanlogtau + log(tau_prime[p]);
+	    }
+	  meantau                    = meantau / (double)c;
+	  meanlogtau                 = meanlogtau / (double)c;
+
+	  tau_prime_bar[t]           = meantau;
+	  bar_log_tau_prime[t]       = meanlogtau;
+
+	  eps[t]                     = 1.0 - (exp(meanlogtau)/meantau);
+
+	  /* if eps  > 0.0001, a cloud has been detected. Set cloud_detected[] to 1 */
+	  if (eps[t] < gThreshold) 
+	    {
+	      cloud_detected[t]      = 0;
+	    }
+	  else
+	    {
+	      cloud_detected[t]      = 1;
+	    }
+	}
+    } /* End PASS 2 */
+
+} /* End of AlexandrovScreen() */
+
+

src/airmass.c

+/*
+ * Atmospheric Sciences Research Center
+ * 100 Fuller Rd.
+ * Albany,  NY 12205
+ */
+
+/*
+ * Version History
+ * ---------------
+ *
+ * Early 1994	First crack.	v 1.0.0
+ *
+ * 5/18/94	Modified for Splus .C interface.  Compile with the Splus 
+ *		COMPILE option (i.e., 'Splus COMPILE airmass.c').  Note that
+ *		you'll have to define _SPLUS_ at compile time to get this
+ *		to work properly.  I did this by changeing the S_makefile in
+ *		$SHOME/newfun/lib so that -D_SPLUS_ was a C compiler flag.
+ *		I also changed CC=gcc so that I could keep things ANSI.
+ *		v 1.0.1		jim
+ * 5/29/97	Added command line date option, as in "airmass 34454.75"
+ *		v 1.0.2		jim
+ * 
+ * Tue Jul  1 13:12:23 EDT 1997
+ *		Added a bunch of other flags including -n, -s, -d.
+ *		v 1.1.0		jim
+ *
+ * 2/20/98	Splus interface is less needed now, but I'll leave it.  Added
+ *		ifdef's to just make the sunae() function visible so that 
+ *		all my programs are using the same source for airmass calcs.
+ *		v 1.1.1		jim
+ *
+ * 2/16/00	fixed y2k leap year bug (rob).  Was missing 2000 as leap year.
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <string.h>
+
+/****** Zebra Includes ******/
+#include "config.h"
+#include "defs.h"
+#include "message.h"
+#include "timer.h"
+#include "DataStore.h"
+
+/****** BW Includes ******/
+#define BW_CODE
+#include "bw_main.h"
+
+typedef struct 
+{
+  int		year, doy;
+  double	hour, lat, lon, az, el, ha, dec, zen, soldst, am;
+} ae_pack;
+
+#ifndef TRUE
+#define TRUE 1
+#define FALSE 0
+#endif
+
+#ifndef M_DTOR
+#define	M_DTOR	0.0174532925199433
+#endif
+#ifndef M_RTOD
+#define	M_RTOD	57.2957795130823230
+#endif
+#ifndef M_HTOR
+#define	M_HTOR	0.2617993877991494
+#endif
+#ifndef M_RTOH
+#define	M_RTOH	3.8197186342054881
+#endif
+#ifndef M_HTOD
+#define	M_HTOD	15.0
+#endif
+#ifndef M_DTOH
+#define	M_DTOH	0.0666666666666667
+#endif
+#ifndef M_2PI
+#define M_2PI      6.2831853071795862320E0 /*Hex  2^ 2 * 1.921FB54442D18 */
+#endif
+
+#define DAYS_1900_1970  25568L   /* days from 1/1/00 to 1/1/1970  */
+
+double airmass(double), Start, End, Increment;
+void   sunae(ae_pack *);
+const char *Version = "1.1.2";
+
+char *CmdLineDate, NoNight, Degrees;
+double Lat=42.698364, Lng=-73.828125;  /* Albany */
+
+/*
+ * if we're compiling for the Splus .C interface, we don't need all
+ * the command line deciphering, usage, and main() crappola.
+ */
+
+char *Progname, Quiet=FALSE;
+char Filename[128];
+FILE *Fp;
+
+int getopt();
+extern char *optarg;
+extern int optind;
+
+
+/* =========================================
+ *
+ * The airmass.c file contains the following
+ * functions:
+ *
+ *    airmass()
+ *    Cairmass()
+ *    sunae() 
+ *
+ * =========================================
+ */
+
+
+/*****************************************************************************
+ *
+ * airmass() function
+ *
+ *****************************************************************************
+ */
+/*
+ *   this function calculates air mass using kasten's
+ *   approximation to bemporad's tables
+ */
+
+double airmass(double el)    /* elevation given in degrees */
+{
+  double  zr;
+ 
+  if(el < 0.0)  /* sun below the horizon */
+    {
+      return(-1.0);
+    }
+ 
+  zr = M_DTOR * (90.0 - el);
+  return (1.0/(cos(zr) + 0.50572*pow(6.07995 + el, -1.6364)));
+ 
+} /* End of airmass() function */
+
+
+/*****************************************************************************
+ *
+ * Cairmass() function
+ *
+ *    uses: sunae() function
+ *
+ *****************************************************************************
+ */
+/* Cairmass - Splus .C function
+ *
+ * times - vector of times in the form 33445.2356 ("double")
+ * am	 - vector of return values for airmass ("double")
+ * az	 - vector of return values for azimuth ("double")
+ * el	 - vector of return values for elevation ("double")
+ * sd	 - vector of return values for solar distance ("double")
+ * lat   - latitude [-90,90] ("double")
+ * lng   - [-180, 180], west negative ("double")
+ * n     - length of times vector ("integer")
+ *
+ */
+ 
+void Cairmass(double *times, 
+	      double *am, 
+	      double *sd, 
+	      double *lat, 
+	      double *lng, 
+	      int    *n)
+{
+  long        i;
+  ae_pack     aep;
+  double      jd, jdays;
+  time_t      secs;
+  struct tm  *tm;
+
+  int         solar_noon_sample;
+  float       solar_noon_hour;
+
+  /* 
+   * first check to see if proper arguments were supplied.  If we proceed
+   * without this check, Splus will probably crash on a Bus Error.
+   */
+
+  if(*lng < -180.0 || *lng > 180.0) 
+    {
+      msg_ELog(EF_INFO, "Cairmass(): longitude is outside of (-180,180): %.4f\n",
+	       *lng);
+      return;
+    }
+
+  if(*lat < -90.0 || *lat > 90.0) 
+    {
+      msg_ELog(EF_INFO, "Cairmass(): latitude is outside of (-90,90): %.4f\n",
+	       *lat);
+      return;
+    }
+
+  aep.lat = *lat;
+  aep.lon = *lng;
+
+  if (*n > 1)
+    {
+      solar_noon_sample = (-1);
+    }
+
+  /*  for(i=0;i<*n;i++)   tracking down "last sample" problem Nov. 11, 2010 */
+  for(i=0;i < *n;i++)   /* tracking down "last sample" problem Nov. 11, 2010 */
+    {
+      jd    = times[i];
+
+      jdays = jd - (double) DAYS_1900_1970;
+
+      if(jdays < 0.0) 
+	{
+	  if (i == *n)
+	    {
+	    }
+	  else
+	    {
+	      msg_ELog(EF_INFO, "Cairmass(): Illegal date in times[%ld]: %.4f (must be >= %ld)\n", 
+		       i, jd, DAYS_1900_1970);
+	    }
+	  return;
+	}
+
+      secs = (long) (jdays * 3600.0 * 24.0);
+
+      tm       = gmtime(&secs);
+      aep.year = tm->tm_year+1900;
+      aep.doy  = tm->tm_yday+1;
+      aep.hour = (jdays - ((int) jdays)) * 24.0;
+
+      sunae(&aep);
+
+      if (aep.az >= 6.28)
+	{
+	  aep.az = aep.az * 0.5;
+	}
+
+      if (aep.az > 3.14 && aep.az < 3.142 && i > 1)
+	{
+	  solar_noon_sample = i;  /* Want a sample for which azimuth is as close to 180 as we will get */
+	  solar_noon_hour   = (float) aep.hour;
+	}
+
+      /* should change parameter am so that it is a matrix - so we can return 
+       * all the other calculations
+       */
+
+      am[i] = aep.am;
+      if(sd != NULL)
+	{
+	  sd[i] = aep.soldst;
+	}
+
+    }
+
+
+} /* End of Cairmass() function */
+
+
+/*****************************************************************************
+ *
+ * sunae() function
+ *
+ *   uses airmass() function
+ *
+ *****************************************************************************
+ */
+
+/*
+ * sunae()
+ *  expects:
+ *
+ *	o GMT
+ *	o latitude 
+ *	o longitude west negative
+ *	o aep->year to be year with century, e.g. 1992
+ *	o aep->doy starting at 1
+ *	o aep->hours to be hours.frac
+ *
+ *  returns:
+ *
+ *	o aep->{az,el,ha,dec,zen} in radians
+ *	o aep->soldst in AU
+ *	o aep->am
+ */
+
+
+void sunae(ae_pack *aep)
+{
+  double	delta, leap, jd, tm, mnlong, mnanom, eclong, oblqec,
+                num, den, ra, lat, sd, cd, sl, cl, gmst, lmst, refrac;
+  double        el_deg, el_rad; 
+
+  delta = aep->year - 1949.0;
+  leap  = (int) (0.25 * delta); 
+  jd    = 32916.5 + delta*365.0 + leap + aep->doy + aep->hour/24.0;
+
+  /*
+   * "the last year of a century is not a leap year therefore"
+   * except for those divisible by 400
+   * 
+   * Actually, if we just did a year % 4 we'd be set for a long time
+   */
+  
+  if (aep->year % 100 == 0 && aep->year % 400 != 0)
+    {
+      jd -= 1.0;
+    }
+
+  tm = jd - 51545.0;
+
+#ifdef DEBUG
+  printf("sunae(): year=%d doy=%d hour=%f jd=%f, tm=%f\n", aep->year, aep->doy, aep->hour, jd, tm); 
+#endif
+
+  mnlong = fmod(280.460 + 0.9856474*tm, 360.0);
+  if (mnlong < 0.0)
+    {
+      mnlong += 360.0;
+    }
+
+  mnanom = fmod(357.528 + 0.9856003*tm, 360.0);
+  mnanom *= M_DTOR;
+
+  eclong = fmod(mnlong + 1.915*sin(mnanom) + 0.020*sin(2.0*mnanom), 360.0);
+  eclong *= M_DTOR;
+
+  oblqec = 23.439 - 0.0000004*tm;
+  oblqec *= M_DTOR;
+
+  num = cos(oblqec)*sin(eclong);
+  den = cos(eclong);
+  ra  = atan(num/den);
+
+  /* force ra between 0 and 2*pi */
+  if (den < 0.0)
+    {
+      ra += M_PI;
+    }
+  else if (num < 0.0)
+    {
+      ra += M_2PI;
+    }
+
+#ifdef DEBUG
+  printf("sunae(): intermediate stuff: tm=%f mnanom=%f eclong=%f oblqec=%f\n", 
+	 tm, mnanom, eclong, oblqec);
+#endif
+
+  /* dec in radians */
+  aep->dec = asin(sin(oblqec)*sin(eclong));
+
+  /* calculate Greenwich mean sidereal time in hours */
+  gmst = fmod(6.697375 + 0.0657098242*tm + aep->hour, 24.0);
+
+  /* calculate local mean sidereal time in radians  */
+  lmst = fmod(gmst + aep->lon*M_DTOH, 24.0);
+  lmst = lmst*M_HTOR;
+
+  /* calculate hour angle in radians between -pi and pi */
+  aep->ha = lmst - ra;
+  if (aep->ha < -M_PI)
+    {
+      aep->ha += M_2PI;
+    }
+  if (aep->ha > M_PI)
+    {
+      aep->ha -= M_2PI;
+    }
+
+  /* change latitude to radians */
+  lat = aep->lat*M_DTOR;
+  sd  = sin(aep->dec);
+  cd  = cos(aep->dec);
+  sl  = sin(lat);
+  cl  = cos(lat);
+
+#ifdef DEBUG
+  printf("sunae(): intermediate stuff: lat=%f sd=%f cd=%f sl=%f cl=%f\n", 
+	 lat,sd,cd,sl,cl);
+#endif
+
+  /* calculate azimuth and elevation */
+  el_rad  = asin(sd*sl + cd*cl*cos(aep->ha));
+  aep->az = asin(-cd*sin(aep->ha)/cos(el_rad));
+
+  /* this puts azimuth between 0 and 2*pi radians */
+  if (sd - sin(el_rad)*sl >= 0.0) 
+    {
+      if(sin(aep->az) < 0.0)
+	{
+	  aep->az += M_2PI;
+	}
+    } 
+  else
+    {
+      aep->az = M_PI - aep->az;
+    }
+
+  /* calculate refraction correction for US stan. atmosphere */
+
+  /* Refraction correction and airmass updated on 11/19/93 */
+
+  /* need to have el in degs before calculating correction */
+  el_deg = el_rad * M_RTOD;
+
+  if (el_deg >= 19.225)
+    {
+      refrac = 0.00452 * 3.51823/tan(el_rad);
+    }
+
+  else if (el_deg > -0.766 && el_deg < 19.225)
+    {
+      refrac = 3.51823*(0.1594 + 0.0196*el_deg + 0.00002 * el_deg*el_deg) /
+	       (1.0 + 0.505 * el_deg + 0.0845 * el_deg*el_deg);
+    }
+  else
+    {
+      refrac = 0.0;
+    }
+
+  /*
+   * old refraction equations
+   *   
+   * if (aep->el > -0.762)
+   *   refrac = 3.51561*(0.1594 + 0.0196*aep->el + 0.00002*aep->el*aep->el) /
+   *     (1.0 + 0.505*aep->el + 0.0845*aep->el*aep->el);
+   * else
+   *   refrac = 0.762;
+   */
+ 
+  /* note that 3.51561=1013.2 mb/288.2 C */
+ 
+  el_deg += refrac;
+
+  /* calculate distance to sun in A.U. & diameter in degs */
+  aep->soldst = 1.00014 - 0.01671*cos(mnanom) - 0.00014*cos(mnanom + mnanom);
+ 
+  /* see if airmass needs doing */
+ 
+  aep->am = airmass(el_deg);   /* use elevation in degrees */
+ 
+  /* convert back to radians - return EVERYTHING in radians */
+
+  aep->el = el_deg * M_DTOR;
+
+  /* do zenith angle */
+  aep->zen = M_PI/2.0 - aep->el;
+ 
+  /* 
+     printf("sunae(): am=%f el=%f sd=%f\n", aep->am, aep->el, aep->soldst); 
+  */
+
+    return; /* 12.0 + aep->ha*M_DTOH; */
+} /* End of sunae() function */
+
+
+
+

src/aod_nimfrsr_batch.pro

+; $Id: mfrod1barnmich_batch.pro,v 1.4 2008-07-24 20:00:29 koontz Exp $
+;--------------------------------------------------------------------------------
+;+
+; Abstract:
+;
+;    This is the top-level IDL commands needed to run quicklook plot
+;    scripts for the MFRSROD1MICHALSKY VAP. 
+;
+;    Environment variables, set by the C-code portion of the VAP, 
+;    are how the arguments are transferred to IDL.
+;
+; Author:
+;    Annette Koontz, PNNL
+;
+; Date Created:
+;    September, 2001
+;
+; Date Last Modified:
+;    $Date: 2008-07-24 20:00:29 $
+;
+;-
+
+pro aod_nimfrsr_batch
+
+;  date_str     = 'date='      + '"' + getenv("DATE")       + '"'
+;  vapname_str  = 'vapname='   + '"' + getenv("VAPNAME")    + '"'
+;  site_str     = 'site='      + '"' + getenv("SITE")       + '"'
+;  facility_str = 'facility='  + '"' + getenv("FACILITY")   + '"'
+;  which_mfr_str = 'which_mfr='    + '"' + getenv ("WHICH_MFR") + '"'
+  date        = long(getenv("DATE"))
+  ;print, 'in batch: date=' + date
+  vapname     = getenv("VAPNAME")
+  ;print, 'in batch: vapname=' + vapname
+  site        = getenv("SITE")
+  ;print, ' in batch: site=' + site
+  facility    = getenv("FACILITY")
+  ;print, ' in batch: facility=' + facility
+  which_mfr   = getenv("WHICH_MFR")
+  ;print, ' in batch: which_mfr=' + which_mfr
+  gNoForgan   = getenv("GNOFORGAN");
+  ;print, ' in batch: gNoForgan=' + gNoForgan
+
+  OpenW, 1, '/tmp/IDL_odvap_batch.tmp'
+; printf, 1, date
+; printf, 1, vapname
+; printf, 1, site
+; printf, 1, facility	
+; printf, 1, which_mfr
+; printf, 1, gNoForgan
+; printf, 1, 'calling plots NOW ... '
+
+  aod_nimfrsr_plots, date, vapname, site, facility, which_mfr, gNoForgan
+
+;  printf, 1, 'back from plots NOW ...'
+
+  return
+
+end

src/aod_nimfrsr_plots.pro

+pro aod_nimfrsr_plots, date, vapname, site, facility, which_mfr, gNoForgan
+
+set_plot, 'Z'
+
+; Checking the keywords
+;if ( n_elements(date) eq 0 ) then date = ''
+;if ( n_elements(facility) eq 0 ) then facility = 'C1'
+;if ( n_elements(site) eq 0 ) then site = 'sgp'
+;if ( n_elements(vapname) eq 0 ) then vapname = 'mfrod1barnmich'
+;if ( n_elements(which_mfr) eq 0) then which_mfr = 'mfrsr'
+
+which_output = 'c1'
+
+dataHome = getenv("DATASTREAM_DATA")
+if(n_elements(dataHome) eq 0) then begin
+     message, ' *** DATASTREAM_DATA is not defined ***'
+endif
+
+; IDL version to print in the left bottom corner
+l_idlVersion = 'IDL Version 8.2'
+l_idlDate    = '$Date: 2012-11-02 $'
+
+l_path       = dataHome + '/' + site + '/' + site + which_mfr + 'aod1mich' +$
+               facility + '.' + which_output
+
+qlookHome = getenv("QUICKLOOK_DATA")
+;print, ' in mfrod1barnmich_plots, qlookHome=' + qlookHome
+if(n_elements(qlookHome) eq 0) then begin
+    print, '        *** QUICKLOOK_DATA is not defined ***'
+endif
+
+; Read in the data 
+
+str_date = string(format='(I0)',date)
+str_site = string(site)
+str_fac  = string(facility)
+str_mfr  = string(which_mfr)
+l_file   = str_site + str_mfr + 'aod1mich' + str_fac + '.' + which_output + '.' + str_date + '*cdf'
+
+look_for = l_path + '/' + l_file
+
+files    = findfile(look_for, count=count)
+
+if (count eq 0) then begin
+   print, ' No file found with name ' + l_file
+   exit, /NO_CONFIRM, status=-1
+endif
+
+l_fid = ncdf_open(files[0])
+
+ncdf_varget, l_fid, 'base_time', l_baseTime
+ncdf_varget, l_fid, 'time_offset', l_timeOffset
+
+ncdf_attget, l_fid, 'process_version', /GLOBAL,  create_version
+create_ver_str = string(create_version)
+ncdf_attget, l_fid, 'history', /GLOBAL,  create_date
+create_date_str = string(create_date)
+;ncdf_attget, l_fid, 'input_data_used', /GLOBAL,  data_used
+data_used_str = 'not used'
+
+ncdf_varget, l_fid, 'lat', l_lat
+ncdf_varget, l_fid, 'lon', l_lon
+
+ncdf_varget, l_fid, 'aerosol_optical_depth_filter1', l_aod_f1
+ncdf_varget, l_fid, 'qc_aerosol_optical_depth_filter1', l_qc_aod_f1
+ncdf_varget, l_fid, 'aerosol_optical_depth_filter2', l_aod_f2
+ncdf_varget, l_fid, 'qc_aerosol_optical_depth_filter2', l_qc_aod_f2
+ncdf_varget, l_fid, 'aerosol_optical_depth_filter3', l_aod_f3
+ncdf_varget, l_fid, 'qc_aerosol_optical_depth_filter3', l_qc_aod_f3
+ncdf_varget, l_fid, 'aerosol_optical_depth_filter4', l_aod_f4
+ncdf_varget, l_fid, 'qc_aerosol_optical_depth_filter4', l_qc_aod_f4
+ncdf_varget, l_fid, 'aerosol_optical_depth_filter5', l_aod_f5
+ncdf_varget, l_fid, 'qc_aerosol_optical_depth_filter5', l_qc_aod_f5
+
+ncdf_varget, l_fid, 'direct_normal_narrowband_filter1', l_dirnor_f1
+ncdf_varget, l_fid, 'qc_direct_normal_narrowband_filter1', l_qc_dirnor_f1
+ncdf_varget, l_fid, 'direct_normal_narrowband_filter2', l_dirnor_f2
+ncdf_varget, l_fid, 'qc_direct_normal_narrowband_filter2', l_qc_dirnor_f2
+ncdf_varget, l_fid, 'direct_normal_narrowband_filter3', l_dirnor_f3
+ncdf_varget, l_fid, 'qc_direct_normal_narrowband_filter3', l_qc_dirnor_f3
+ncdf_varget, l_fid, 'direct_normal_narrowband_filter4', l_dirnor_f4
+ncdf_varget, l_fid, 'qc_direct_normal_narrowband_filter4', l_qc_dirnor_f4
+ncdf_varget, l_fid, 'direct_normal_narrowband_filter5', l_dirnor_f5
+ncdf_varget, l_fid, 'qc_direct_normal_narrowband_filter5', l_qc_dirnor_f5
+
+if (which_mfr eq 'mfrsr') then begin
+  ncdf_varget, l_fid, 'diffuse_hemisp_narrowband_filter1', l_diffuse_f1
+  ncdf_varget, l_fid, 'qc_diffuse_hemisp_narrowband_filter1', l_qc_diffuse_f1
+  ncdf_varget, l_fid, 'diffuse_hemisp_narrowband_filter2', l_diffuse_f2
+  ncdf_varget, l_fid, 'qc_diffuse_hemisp_narrowband_filter2', l_qc_diffuse_f2
+  ncdf_varget, l_fid, 'diffuse_hemisp_narrowband_filter3', l_diffuse_f3
+  ncdf_varget, l_fid, 'qc_diffuse_hemisp_narrowband_filter3', l_qc_diffuse_f3
+  ncdf_varget, l_fid, 'diffuse_hemisp_narrowband_filter4', l_diffuse_f4
+  ncdf_varget, l_fid, 'qc_diffuse_hemisp_narrowband_filter4', l_qc_diffuse_f4
+  ncdf_varget, l_fid, 'diffuse_hemisp_narrowband_filter5', l_diffuse_f5
+  ncdf_varget, l_fid, 'qc_diffuse_hemisp_narrowband_filter5', l_qc_diffuse_f5
+endif else begin
+  l_diffuse_f1 = -9999;
+  l_diffuse_f2 = -9999;
+  l_diffuse_f3 = -9999;
+  l_diffuse_f4 = -9999;
+  l_diffuse_f5 = -9999;
+  l_qc_diffuse_f1 = -9999;
+  l_qc_diffuse_f2 = -9999;
+  l_qc_diffuse_f3 = -9999;
+  l_qc_diffuse_f4 = -9999;
+  l_qc_diffuse_f5 = -9999;
+endelse
+
+ncdf_varget, l_fid, 'angstrom_exponent', l_angstrom
+ncdf_varget, l_fid, 'qc_angstrom_exponent', l_qc_angstrom
+
+; get Io times and Io (interquartile and gauss) values from the netcdf
+; filen
+
+ncdf_varget, l_fid, 'Io_interquartile_time',   l_interq_times
+ncdf_varget, l_fid, 'Io_interquartile_values', l_interq_values
+ncdf_varget, l_fid, 'Io_gauss_time',           l_gauss_times
+ncdf_varget, l_fid, 'Io_gauss_values',         l_gauss_values
+ncdf_varget, l_fid, 'Io_filter1',              l_Io_filter1
+ncdf_varget, l_fid, 'Io_filter2',              l_Io_filter2
+ncdf_varget, l_fid, 'Io_filter3',              l_Io_filter3
+ncdf_varget, l_fid, 'Io_filter4',              l_Io_filter4
+ncdf_varget, l_fid, 'Io_filter5',              l_Io_filter5
+
+; Get Io data arranged into individual filter arrays
+;
+n_Io_pts           = n_elements(l_interq_times) ; we will make this smaller as we weed out unset values
+n_Io_interquartile = n_Io_pts
+n_gauss_pts        = n_elements(l_gauss_times)
+n_Io_gauss         = n_Io_pts
+
+; First do the interquartile
+
+l_Io_iq_time       = fltarr(n_Io_pts)
+l_Io_iq_f1         = fltarr(n_Io_pts)
+l_Io_iq_f2         = fltarr(n_Io_pts)
+l_Io_iq_f3         = fltarr(n_Io_pts)
+l_Io_iq_f4         = fltarr(n_Io_pts)
+l_Io_iq_f5         = fltarr(n_Io_pts)
+
+n_iq               = 0
+
+l_Io_iq_time[0]    = l_interq_times[0]
+
+; Send only the data we intend to plot to the do_io_plot routine
+;
+; when we are running the VAP using an ascii Io file (the -Z option)
+; it would be nice to limit the plot to +/- 20 days around the rundate
+;
+
+;print, ' In aod_nimfrsr_plots'
+;print, ' ===================='
+
+;print, ' gNoForgan = ', gNoForgan
+;print, ' ----------------------'
+
+if (gNoForgan EQ 1) then begin
+  ; use a different function for the Io plot (gauss only)
+  ; aaa
+
+  ; we went to send only the times +/- 20 days around the gauss_time
+  ; that is closes to l_basetime
+
+  found_it = (-1)
+
+;  print, ' l_gauss_times has ', n_elements(l_gauss_times), ' elements' 
+
+;  for n = 0, (n_elements(l_gauss_times)-1) do begin
+;    print, 'n, l_gauss_time[n], l_basetime: ', n, long(l_gauss_times[n]), long(l_basetime)
+;  endfor
+
+;  for n = 0, (n_elements(l_gauss_times)-2) do begin
+;    print, 'n, l_gauss_time[n], l_basetime: ', n, long(l_gauss_times[n]), long(l_basetime)
+;
+;    if (l_gauss_times[n] LE l_basetime) then begin
+;       print, ' n, l_gauss_times[n]=', n, long(l_gauss_times[n])
+;       print, ' n+1, l_gauss_times[n+1], l_base_time=', (n+1), long(l_gauss_times[n+1]), long(l_basetime)
+;      if (l_gauss_times[n+1] GT l_basetime) then begin
+;        found_it = n
+;        break
+;      endif
+;    endif
+;    if (found_it GE 0) then break;
+;  endfor
+
+;  print, ' will use gauss time: ', found_it
+
+  ; We want to plot (if possible) +/- 30 days around the rundate
+
+;  if (found_it GT 31) then begin
+;    ifirst = found_it - 30
+;  endif else begin
+;    ifirst = 0
+;  endelse
+
+;  if ( (n_gauss_pts - found_it) LE 30) then begin
+;     ilast = n_gauss_pts
+;  endif else begin
+;     ilast = found_it + 30
+;  endelse
+
+;  print, ' setting up gauss: ifirst, ilast=', ifirst, ilast
+
+  ;print, 'SCALE: l_Io_filter1, etc.= ', l_Io_filter1,l_Io_filter2,l_Io_filter3,l_Io_filter4,l_Io_filter5
+
+  l_Io_gauss_f1 = fltarr(n_elements(l_gauss_times))
+  l_Io_gauss_f2 = fltarr(n_elements(l_gauss_times))
+  l_Io_gauss_f3 = fltarr(n_elements(l_gauss_times))
+  l_Io_gauss_f4 = fltarr(n_elements(l_gauss_times))
+  l_Io_gauss_f5 = fltarr(n_elements(l_gauss_times))
+
+  for n = 0, (n_elements(l_gauss_times)-1) do begin
+    l_Io_gauss_f1[n]     = l_gauss_values[0,n] / l_Io_filter1
+    l_Io_gauss_f2[n]     = l_gauss_values[1,n] / l_Io_filter2
+    l_Io_gauss_f3[n]     = l_gauss_values[2,n] / l_Io_filter3
+    l_Io_gauss_f4[n]     = l_gauss_values[3,n] / l_Io_filter4
+    l_Io_gauss_f5[n]     = l_gauss_values[4,n] / l_Io_filter5
+  endfor
+
+  ;print, 'l_baseTime = ', l_baseTime
+
+;  print, ' CALLING DO_IO_PLOT2'
+  device,set_resolution=[1024,768]
+  erase
+  desire_png = 1
+
+  do_io_plot2, which_mfr, date, create_ver_str, create_date_str, l_idlVersion,$
+    l_baseTime, l_gauss_times, $
+    l_Io_gauss_f1, l_Io_gauss_f2, l_Io_gauss_f3, l_Io_gauss_f4, l_Io_gauss_f5
+
+    vap_lib = getenv("VAP_HOME") + '/lib'
+    tmp_vapname = str_mfr + 'aod1mich' 
+
+    get_quicklook_dir, date, tmp_vapname, site, facility, which_output, $
+                   /create,  quicklook_dir = pngpath, $
+                   errornum = enum
+
+    pngname = str_site + str_mfr + 'aod1mich' + str_fac + which_output + '.' + string(format='(I0)',date) + '.000000.' + 'Io.png'
+
+    desire_png = 1
+    if(desire_png eq 1) then begin
+      name = pngpath + '/' + pngname
+      tvlct, r,g,b, /GET
+      write_png, name, tvrd(), r,g,b
+  endif
+
+endif else begin
+
+;  print, 'gNoForgan is: ', gNoForgan
+;  print, '====================='
+
+;  print, 'n_Io_pts=', n_Io_pts
+
+  for n = 0, n_Io_pts-1 do begin
+      l_Io_iq_time[n] = l_interq_times[n] ; time is relative to first interquartile time
+;      print, '   ... l_Io_iq_time[n] = ', l_Io_iq_time[n], ' ... '
+
+      l_Io_iq_f1[n]   = l_interq_values[0, n] / l_Io_filter1
+      l_Io_iq_f2[n]   = l_interq_values[1, n] / l_Io_filter2
+      l_Io_iq_f3[n]   = l_interq_values[2, n] / l_Io_filter3
+      l_Io_iq_f4[n]   = l_interq_values[3, n] / l_Io_filter4
+      l_Io_iq_f5[n]   = l_interq_values[4, n] / l_Io_filter5
+      n_iq            = n_iq + 1
+  endfor
+
+;  print, 'n_gauss_pts=', n_gauss_pts
+
+  l_Io_gauss_time      = fltarr(n_gauss_pts)
+  l_Io_gauss_f1        = fltarr(n_gauss_pts)
+  l_Io_gauss_f2        = fltarr(n_gauss_pts)
+  l_Io_gauss_f3        = fltarr(n_gauss_pts)
+  l_Io_gauss_f4        = fltarr(n_gauss_pts)
+  l_Io_gauss_f5        = fltarr(n_gauss_pts)
+  n_gauss              = 0
+  l_Io_gauss_time[0]   = l_gauss_times[0]
+  for n = 0, n_gauss_pts-1 do begin
+    if (l_gauss_values[0,n] GT 0) then begin
+      l_Io_gauss_time[n] = l_gauss_times[n]  ; time is relative to first interquartile time
+;      print, ' n=', n, ' gauss_time = ', l_Io_gauss_time[n]
+
+      l_Io_gauss_f1[n]   = l_gauss_values[0, n] / l_Io_filter1
+      l_Io_gauss_f2[n]   = l_gauss_values[1, n] / l_Io_filter2
+      l_Io_gauss_f3[n]   = l_gauss_values[2, n] / l_Io_filter3
+      l_Io_gauss_f4[n]   = l_gauss_values[3, n] / l_Io_filter4
+      l_Io_gauss_f5[n]   = l_gauss_values[4, n] / l_Io_filter5
+      n_gauss            = n_gauss + 1
+    endif
+  endfor
+
+  ;print, 'l_baseTime = ', l_baseTime
+
+  device,set_resolution=[1024,768]
+  erase
+  desire_png = 1
+
+  do_io_plot, which_mfr, date, create_ver_str, create_date_str, l_idlVersion,$
+    l_baseTime, $
+    n_iq,    l_Io_iq_time,    l_Io_iq_f1,    l_Io_iq_f2,    l_Io_iq_f3,    l_Io_iq_f4,    l_Io_iq_f5, $
+    n_gauss, l_Io_gauss_time, l_Io_gauss_f1, l_Io_gauss_f2, l_Io_gauss_f3, l_Io_gauss_f4, l_Io_gauss_f5
+
+    vap_lib = getenv("VAP_HOME") + '/lib'
+    tmp_vapname = str_mfr + 'aod1mich' 
+
+    get_quicklook_dir, date, tmp_vapname, site, facility, which_output, $
+                   /create,  quicklook_dir = pngpath, $
+                   errornum = enum
+
+    pngname = str_site + str_mfr + 'aod1mich' + str_fac + which_output + '.' + string(format='(I0)',date) + '.000000.' + 'Io.png'
+
+    desire_png = 1
+    if(desire_png eq 1) then begin
+      name = pngpath + '/' + pngname
+      tvlct, r,g,b, /GET
+      write_png, name, tvrd(), r,g,b
+  endif
+ endelse
+
+ncdf_attget, l_fid, 'command_line', /GLOBAL, Command_line
+Comm_str     = string(command_line)
+
+comm_parts   = strsplit(Comm_str, ' ', /EXTRACT)
+num_pts      = n_elements(comm_parts)
+
+match        = 0
+for i=0,num_pts-1 do begin
+   result = strcmp(comm_parts[i], '-c')
+   if (result eq 1) then begin
+       match = i+1
+       break
+   endif
+endfor
+
+if (match eq 0) then begin
+   channels = '1,5'
+endif else begin
+   channels      = strsplit (comm_parts[match], ',', /EXTRACT)
+endelse
+
+reads, channels, channel1, format='(i5)'
+
+; which_langley = 1 is michalsky Langleys
+;               = 2 is barnard Langleys
+;               = 0 is both michalsky and barnard langleys
+
+which_langley = 1
+
+for i=0,num_pts-1 do begin
+   result = strcmp(comm_parts[i], '-b')
+   if (result eq 1) then begin
+       which_one = comm_parts[i+1]
+       result = strcmp(which_one, 'barnmich')
+       if (result eq 1) then begin
+          which_langley = 0
+       endif else begin
+          result = strcmp(which_one, 'barn')
+          if (result eq 1) then begin
+             which_langley = 2
+          endif
+       endelse
+   endif
+endfor
+
+;if (ncdf_varid(l_fid, 'direct_normal_narrowband_filter2') ne -1) then begin
+;    ncdf_varget, l_fid, 'direct_normal_narrowband_filter2', l_direct_beam
+;endif
+
+if (ncdf_varid(l_fid, 'lat') ne -1) then begin
+    ncdf_varget, l_fid, 'lat', l_lat
+endif
+if (ncdf_varid(l_fid, 'lon') ne -1) then begin
+    ncdf_varget, l_fid, 'lon', l_lon
+endif
+
+; close the netCDF file
+ncdf_close, l_fid
+
+num_samples  = n_elements(l_aod_f2)
+
+; Store the yy, mm, dd (year/month/day) and hh, mm, ss (hh:mm:ss)
+; stuff in arrays.
+
+systime2ymdhms, l_baseTime + l_timeOffset, year, mon, day, hh, mm, ss
+time_samples  = make_array(num_samples, /FLOAT, value=0.0)
+
+flthr  = float(hh)
+fltmin = float(mm)
+fltsec = float(ss)
+
+; The DQ office requests that we plot time axes using hh:mm format
+; It turns out that IDL specifies a julian day of 0.0 is really
+; noon on that day, so if we add 0.5, we should get the
+; time axis with hours 00:00, 04:00, ...)
+;
+; This "seems" to work for Davos, Switzerland data. I will have
+; to see what happens with TWP data
+;
+
+my_juljan01 = julday(1, 1, year[0]) + 0.5
+my_julday   = julday(mon[0], day[0], year[0]) - my_juljan01 + 1
+
+for i=0,num_samples-1 do begin
+; want time in terms of  days (day.frac_day), as a float
+
+  dayfrac         = float(((hh[i]*3600.0) + float(mm[i] * 60.0) + float(ss[i]) )/ 86400.0)
+  time_samples[i] = float(my_julday) + dayfrac
+  if (i gt 0) then begin
+    if (time_samples[i] lt time_samples[i-1]) then begin
+       time_samples[i] = time_samples[i] + 1
+    endif
+  endif
+
+endfor
+
+; Use a C function to compute sunrise and sunset so that we 
+; can put all the daylight hours in the netCDF file, and so
+; we can center the plots properly.
+;
+
+;
+;   C code to call the function is shown below.
+;   
+;      year, month and day are "int"
+;      lon (+ is east) and lat (+ is north) are double
+;      rise and set are double
+;
+;            rs   = __sunriset__         ( year, month, day, 
+;		                           lon, lat,
+;                                          (-35./60.), 1,     /* So we get back sunrise and sunset, not some other stuff */
+;                                          &rise, &set );
+;
+
+year      = double(year[0])
+month     = double(mon[0])
+day       = double(day[0])
+longitude = double(l_lon)
+latitude  = double(l_lat)
+sunrise   = double(0.0)
+sunset    = double(0.0)
+
+
+; I do not think we use sunrise/sunset in these plots anymore
+; and in fact, we seem to have lost the source code for the
+; library, so I can't rebuild the library ...
+;
+; A. Koontz, PNNL, March, 2011
+
+;result = call_external (vap_lib + '/libsunriset.so', $
+;                        'idl_sun_rise_set', $
+;                        year, month, day, $
+;                        longitude, latitude, $
+;                        sunrise, sunset)
+
+; Verify the directory exists, and if not, create it
+command = 'if (! -d ' + pngpath + ') mkdir -m 775 -p ' + pngpath
+spawn, command
+vap = 'mfrod1barnmich'
+
+device,set_resolution=[1024,768]
+erase
+desire_png = 1
+
+do_od_plot1, which_mfr, $
+     date, time_samples, create_ver_str, create_date_str, l_idlVersion, $
+     Langley_used_str, data_used_str, $
+     l_aod_f2,  $
+     l_qc_aod_f2,  $
+     l_aod_f5,  $
+     l_qc_aod_f5, $
+;     sunrise, sunset, $
+     l_angstrom, l_qc_angstrom, $
+     l_lat,      l_lon             ; need lat/lon for DQ-style plots
+
+pngname = str_site + str_mfr + 'aod1mich' + str_fac + which_output + '.' + string(format='(I0)',date) + '.000000.' + 'aod_f2_f5.png'
+
+desire_png = 1
+if(desire_png eq 1) then begin
+   name = pngpath + '/' + pngname
+   tvlct, r,g,b, /GET
+   write_png, name, tvrd(), r,g,b
+endif
+
+device,set_resolution=[1024,768]
+erase
+desire_png = 1
+
+do_od_plot2, which_mfr, $
+     date, time_samples, create_ver_str, create_date_str, l_idlVersion, $
+     Langley_used_str, data_used_str, $
+     l_aod_f1,       l_aod_f2,       l_aod_f3,       l_aod_f4,       l_aod_f5, $
+     l_qc_aod_f1,    l_qc_aod_f2,    l_qc_aod_f3,    l_qc_aod_f4,    l_qc_aod_f5, $
+;     sunrise, sunset, $
+     l_angstrom, l_qc_angstrom, $
+     l_lat,      l_lon     ; lat/lon needed for DQ style plots
+
+pngname = str_site + str_mfr + 'aod1mich' + str_fac + which_output + '.' + string(format='(I0)',date) + '.000000.' + 'aodall.png'
+
+desire_png = 1
+if(desire_png eq 1) then begin
+   name = pngpath + '/' + pngname
+   tvlct, r,g,b, /GET
+   write_png, name, tvrd(), r,g,b
+endif
+
+device,set_resolution=[1024,768]
+erase
+desire_png = 1
+
+do_od_plot3, which_mfr, $
+     date, time_samples, create_ver_str, create_date_str, l_idlVersion, $
+     Langley_used_str, data_used_str, $
+     l_dirnor_f1,     l_dirnor_f2,     l_dirnor_f3,     l_dirnor_f4,     l_dirnor_f5, $
+     l_qc_dirnor_f1,  l_qc_dirnor_f2,  l_qc_dirnor_f3,  l_qc_dirnor_f4,  l_qc_dirnor_f5, $
+     l_diffuse_f1,    l_diffuse_f2,    l_diffuse_f3,    l_diffuse_f4,    l_diffuse_f5, $
+     l_qc_diffuse_f1, l_qc_diffuse_f2, l_qc_diffuse_f3, l_qc_diffuse_f4, l_qc_diffuse_f5, $
+;     sunrise, sunset, $
+     l_angstrom, l_qc_angstrom, $
+     l_lat, l_lon           ; lat/lon needed for DQ style plots
+pngname = str_site + str_mfr + 'aod1mich' + str_fac + which_output + '.' + string(format='(I0)',date) + '.000000.' + 'irrad.png'
+
+desire_png = 1
+if(desire_png eq 1) then begin
+   name = pngpath + '/' + pngname
+   tvlct, r,g,b, /GET
+   write_png, name, tvrd(), r,g,b
+endif
+
+
+
+
+desire_PS = 1
+
+
+
+
+end

src/colors.include

+; $Id: colors.include,v 1.1 2010-11-09 16:27:53 koontz Exp $
+;
+;  This file was meant to be included into IDL source (using the "@" sign
+;       if you want these very discrete colors, which are added to the current
+;	color table.  The version that is executable is color_syms.pro
+
+	; Set up our own colors in the current color table
+  dred=  [0,255,0,255,0,255,0,255,0,255,255,112,219,127,0,255,255,0,155,255,240,185,0]
+  dgreen=[0,0,255,255,255,0,0,255,0,187,127,219,112,127,163,171,153,204,0,255,240,110,150]
+  dblue= [0,255,255,0,0,0,255,255,115,0,127,147,219,127,255,127,0,0,200,230,240,70,0]
+
+  tvlct,dred,dgreen,dblue
+
+	; For the user 
+black = 0
+magenta = 1
+cyan = 2
+yellow = 3
+green = 4
+red = 5
+blue = 6
+white = 7
+navy = 8
+gold = 9
+pink = 10
+aquamarine = 11
+orchid = 12
+gray = 13
+sky = 14
+beige = 15
+orange=16
+darkgreen=17
+purple=18
+daylight=19
+night=20
+brown=21
+pinegreen=22
+  if(!d.name eq 'PS') then pmain = black else pmain = white
+
+  	; Set some personal plotting symbols
+  	;   You use them by using the command      usersym,___,/fill
+  	;   Then, use the command 		     plot,a,b,psym=8
+  d_triangle = [[0,1.5],[1.3,-0.75],[-1.3,-0.75]]
+  d_inv_triangle = [[0,-1.5],[1.3,0.75],[-1.3,0.75]]
+  d_square = [[1,1],[1,-1],[-1,-1],[-1,1]]
+  d_diamond = [[1,0],[0,-1],[-1,0],[0,1]]
+  foo = findgen(17) * (!PI*2/16.0)
+  d_circle = [transpose(cos(foo)),transpose(sin(foo))]

src/day_night_background.pro

+; docformat = 'rst'
+
+;+
+; This procedure will add day and night background color to an existing
+; plot.  It will use POLYFILL IDL procedure, which will plot over
+; any existing plotted information.  Use after plot has been created
+; and before plotting data with OPLOT.
+;
+; :Post: 
+;    This will only work with the native IDL time of Julian days.
+;
+; :Params:
+;   lattitude : in, required, type="scalar"
+;     Lattitude of instrument.
+;   longitude : in, required, type="scalar"
+;     Longitude of instrument.
+;
+; :Keywords:
+;    NOON : in, optional, type=flag
+;      If set will add vertical dashed line at local solar noon.
+;
+; :Uses:
+;    zensun.pro
+;
+; :Examples:
+;     day_night_background, lat, lon, /NOON
+; 
+; :Author: Ken Kehoe, ARM Data Quality Office, University of Oklahoma
+; :History: Created on January 20, 2010
+; :Version: $Id: day_night_background.pro,v 1.1 2010-11-09 16:28:49 koontz Exp $:
+;
+;-
+PRO day_night_background, lattitude, longitude, NOON=show_noon
+
+lon=longitude
+lat=lattitude
+@colors.include
+POLYFILL, [!X.CRANGE[0],!X.CRANGE[0],!X.CRANGE[1],!X.CRANGE[1]], $
+	[!Y.CRANGE[0],!Y.CRANGE[1],!Y.CRANGE[1],!Y.CRANGE[0]],color=night
+
+FOR multi_date=0, !X.CRANGE[1]-!X.CRANGE[0]-1 DO BEGIN
+ FOR jj=-1, 1 DO BEGIN
+  p_date=!X.CRANGE[0]+jj+multi_date
+  CALDAT, p_date, p_month,p_day,p_year
+  p_dayOfYear=FLOAT(p_date - JULDAY(1,0,p_year,0,0,0))
+  zensun, p_dayOfYear,12.0, lat, lon, zenith, azimuth, solfac,	$
+ 	sunrise=sunrise,sunset=sunset,noon=noon
+  noon = JULDAY(1,p_dayOfYear,p_year,0,0,noon*60.*60.)
+  sunrise = JULDAY(1,p_dayOfYear,p_year,0,0,sunrise*60.*60.)
+  sunset= JULDAY(1,p_dayOfYear,p_year,0,0,sunset*60.*60.)
+  index=WHERE(sunrise LT !X.CRANGE[0],indexCt)
+  IF(indexCt GT 0) THEN sunrise[index]=!X.CRANGE[0]
+  index=WHERE(sunrise GT !X.CRANGE[1],indexCt)
+  IF(indexCt GT 0) THEN sunrise[index]=!X.CRANGE[1]
+  index=WHERE(sunset LT !X.CRANGE[0],indexCt)
+  IF(indexCt GT 0) THEN sunset[index]=!X.CRANGE[0]
+  index=WHERE(sunset GT !X.CRANGE[1],indexCt)
+  IF(indexCt GT 0) THEN sunset[index]=!X.CRANGE[1]
+  POLYFILL, [sunrise,sunrise,sunset,sunset],$
+ 	[!Y.CRANGE[0],!Y.CRANGE[1],!Y.CRANGE[1],!Y.CRANGE[0]],COLOR=daylight
+  IF(N_ELEMENTS(show_noon) GT 0) THEN OPLOT, [noon,noon],[!Y.CRANGE[0],!Y.CRANGE[1]], $
+ 	COLOR=yellow,THICK=1,LINESTYLE=5 
+ ENDFOR ; jj FOR
+ENDFOR ; multi_date FOR
+
+END ; Procedure End

src/do_io_plot.pro

+pro do_io_plot, which_mfr, date, create_ver_str, create_date_str, l_idlVersion, $
+  l_base_time, $
+  n_iq,    l_Io_iq_time,    l_Io_iq_f1,    l_Io_iq_f2,    l_Io_iq_f3,    l_Io_iq_f4,    l_Io_iq_f5, $
+  n_gauss, l_Io_gauss_time, l_Io_gauss_f1, l_Io_gauss_f2, l_Io_gauss_f3, l_Io_gauss_f4, l_Io_gauss_f5
+    
+@colors.include           ; DQ office color table
+
+usersym, d_circle, /fill
+
+Compile_Opt DEFINT32
+
+str_tmp         = string(format='(I0)',date) 
+str_year        = strmid(str_tmp, 0, 4)
+str_month       = strmid(str_tmp, 4, 2)
+str_day         = strmid(str_tmp, 6, 2)
+
+str_date2       = str_year+'/'+str_month+'/'+str_day
+
+; This plot region is only for getting a legible title, since we set the
+; x.charsize in the first real plot to a very small value (so as to not get
+; label on our x-axis ticks.
+
+title_str       = strupcase(which_mfr) + ' Interquartile and Gaussian Io '+str_date2
+
+;!p.Multi        = [0, 1, 3, 0, 1]
+!p.region       = [0.05, 0.1, 1, 1.]
+!y.charsize     = 3.0
+!x.ticklen      = -0.02
+!y.ticklen      = -0.005
+!y.tickformat   = '(f6.2)'
+
+; First real plot (Io_interquartile)
+
+;print, 'IN DO_IO_PLOT...'
+;print, 'n_elements(l_Io_iq_time)   ', n_elements(l_Io_iq_time)
+;print, 'n_elements(l_Io_gauss_time)', n_elements(l_Io_gauss_time)
+xmin1 = 2500000000000
+xmax1 = 0
+for n = 0, (n_elements(l_Io_iq_time)-1) do begin
+
+;   print, 'n=', n, ' io_iq_time = ', long(l_Io_iq_time[n])
+
+   if (l_Io_iq_time[n] LE 0.0) then break
+   if (l_Io_iq_time[n] GT 0.0) then begin
+      if (l_Io_iq_time[n] LT xmin1) then begin
+          xmin1 = l_Io_iq_time[n]
+      endif
+   endif
+   if (l_Io_iq_time[n] GT xmax1) then begin
+      xmax1    = l_Io_iq_time[n]
+   endif
+endfor
+
+;print, ' l_io_iq_time min is ', long(xmin1)
+;print, ' l_io_iq_time max is ', long(xmax1)
+
+xmin2 = 2500000000000
+xmax2 = 0
+
+;print, ' n_elements(l_Io_gauss_time) = ', n_elements(l_Io_gauss_time)
+
+for n = 0, (n_elements(l_Io_gauss_time)-1) do begin
+;   print, 'n=', n, ' gauss_time=', long(l_Io_gauss_time[n])
+   if (l_Io_gauss_time[n] LE 0.0) then break
+   if (l_Io_gauss_time[n] GT 0.0) then begin
+      if (l_Io_gauss_time[n] LT xmin2) then begin
+          xmin2 = l_Io_gauss_time[n]
+      endif
+   endif
+   if (l_Io_gauss_time[n] GT xmax2) then begin
+      xmax2    = l_Io_gauss_time[n]
+   endif
+endfor
+
+;print, ' l_io_gauss_time min is ', long(xmin2)
+;print, ' l_io_gauss_time max is ', long(xmax2)
+
+xmin = xmin2
+if (xmin1 LT xmin2) then begin
+  xmin = xmin1
+endif
+
+xmax = xmax2
+if (xmax1 GT xmax2) then begin
+  xmax = xmax1
+endif
+
+;print, 'OVERALL: xmin, xmax=', long(xmin), long(xmax)
+
+; scale the times 
+
+npts_Iq = (-1)
+tmp     = 0L
+for n = 0, (n_elements(l_Io_iq_time)-1) do begin
+   if (l_io_iq_time[n] GT 0) then begin
+
+      npts_Iq = npts_Iq + 1
+
+      tmp = ((l_Io_iq_time[n] - xmin) / 86400)
+;      print, ' l_Io_iq_time[n], xmin, tmp= ', long(l_Io_iq_time[n]), long(xmin), tmp
+      l_Io_iq_time[n] = (tmp - 31)
+   endif
+endfor
+
+; Find the "closest" l_io_gauss_time to l_base_time
+
+found_it = (-1)
+
+for n = 0, (n_elements(l_Io_gauss_time)-2) do begin
+;  print, 'n, l_Io_gauss_time[n], l_base_time: ', n,  long(l_Io_gauss_time[n]), l_base_time
+
+  if (l_Io_gauss_time[n] LE l_base_time) then begin
+     if (l_Io_gauss_time[n+1] GE l_base_time) then begin
+        found_it = n
+        break
+     endif
+  endif
+  if (found_it GE 0) then break;
+endfor
+
+;print, ' will use gauss time: ', found_it
+
+npts_gauss = (-1)
+for n = 0, (n_elements(l_Io_gauss_time)-1) do begin
+   if (l_Io_gauss_time[n] GT 0) then begin
+      npts_gauss = npts_gauss + 1
+
+      tmp = ((l_Io_gauss_time[n] -  xmin) / 86400)
+;      print, ' n=', n, ' l_Io_gauss_time[n], xmin, tmp, (tmp-31) = ', long(l_Io_gauss_time[n]), xmin, tmp, (tmp-31)
+      l_Io_gauss_time[n] = (tmp - 31)
+   endif
+endfor
+
+;print, 'xmin, xmax=', xmin, xmax
+
+ymin1 = min(l_Io_iq_f2)
+ymax1 = max(l_Io_iq_f5)
+ymin2 = min(l_Io_gauss_f2)
+ymax2 = max(l_Io_gauss_f5)
+
+;print, 'ymin1, ymax1=', ymin1, ymax1
+;print, 'ymin2, ymax2=', ymin2, ymax2
+
+if (ymin1 LE 0 AND ymax1 LE 0) then begin
+   ymin = ymin2
+   ymax = ymax2
+endif else begin
+   ymin  = min(ymin1, ymin2)
+   ymax  = max(ymax1, ymax2)
+endelse
+
+;print, 'ymin, ymax=', ymin, ymax
+
+; data is normalized to the Io_filter* values, so a scale of 0.95 to
+; 1.05 should be sufficient.
+
+ymin  = 0.95
+ymax  = 1.05
+
+!x.charsize     = 1.5
+!y.charsize     = 1.5
+
+plot,l_Io_iq_time, $
+     title      = title_str, $
+     ystyle     = 1, $ 
+     xtitle     = 'Time, days', $
+     ytitle     = 'Normalized Io Values', $
+     background = white, $
+     color      = black, $
+     yrange     = [ymin, ymax], $
+     xstyle     = 1, $
+     xrange     = [-31, 31.0], $
+     /NODATA
+
+  oplot, l_Io_iq_time[0:npts_iq], l_Io_iq_f1[0:npts_iq], color=blue,   psym=1, symsize=1  
+  oplot, l_Io_iq_time[0:npts_iq], l_Io_iq_f2[0:npts_iq], color=green,  psym=1, symsize=1  
+  oplot, l_Io_iq_time[0:npts_iq], l_Io_iq_f3[0:npts_iq], color=yellow, psym=1, symsize=1  
+  oplot, l_Io_iq_time[0:npts_iq], l_Io_iq_f4[0:npts_iq], color=orange, psym=1, symsize=1  
+  oplot, l_Io_iq_time[0:npts_iq], l_Io_iq_f5[0:npts_iq], color=red,    psym=1, symsize=1  
+
+  oplot, l_Io_gauss_time[0:npts_gauss], l_Io_gauss_f1[0:npts_gauss], color = blue, linestyle=0 
+  oplot, l_Io_gauss_time[0:npts_gauss], l_Io_gauss_f2[0:npts_gauss], color = green, linestyle=0
+  oplot, l_Io_gauss_time[0:npts_gauss], l_Io_gauss_f3[0:npts_gauss], color = yellow, linestyle=0
+  oplot, l_Io_gauss_time[0:npts_gauss], l_Io_gauss_f4[0:npts_gauss], color = orange, linestyle=0
+  oplot, l_Io_gauss_time[0:npts_gauss], l_Io_gauss_f5[0:npts_gauss], color = red, linestyle=0
+
+  aline    = fltarr(2)
+  bline    = fltarr(2)
+  
+; define a vertical line where the "base_time" is, and draw a dashed
+; line there
+
+  if (found_it GT 0) then begin
+
+;     print, ' location of vertical line = ', l_Io_gauss_time[found_it] 
+
+    aline[0] = l_Io_gauss_time[found_it]
+    aline[1] = l_Io_gauss_time[found_it]
+    bline[0] = ymin
+    bline[1] = ymax
+
+    oplot, aline[0:1], bline[0:1], color = black, linestyle=1
+
+    xyouts, aline[0], (ymin + 0.01), str_date2, orientation=90., color=black, charsize=0.8
+
+ endif
+
+ ; generate a horizontal line  at the Io=1.0 location
+
+    aline[0] = l_Io_gauss_time[0]
+    aline[1] = l_Io_gauss_time[npts_gauss]
+    bline[0] = 1.0
+    bline[1] = 1.0
+
+    oplot, aline[0:1], bline[0:1], color = black, linestyle=1
+
+;    xyouts, aline[0], (ymin + 0.01), str_date2, orientation=90., color=black, charsize=0.8
+
+
+; legend info copied from do_od_plot2.pro
+;=============================================================
+
+; annotate all plots in the plot at the bottom
+
+;!p.region = [0., 0., 1., .1]
+!x.charsize     = .5
+!y.charsize     = .5
+
+
+xline    = fltarr(2)
+yline    = fltarr(2)
+
+xyouts, .20,  0.345,'Interquartile Io', color=black, /normal, charsize=1.0
+xyouts, .70,  0.345,'Gaussian Io', color=black, /normal, charsize=1.0
+
+xyouts, .25, 0.32, '+' , color=blue, /normal, charsize=1.0
+xyouts, .30, 0.32, ' 415 nm (filter 1)', color=black, /normal, charsize=0.8
+xline[0] = 0.75
+yline[0] = 0.32
+xline[1] = 0.80
+yline[1] = 0.32
+plots, xline[0], yline[0], /normal, color=blue
+plots, xline[1], yline[1], /normal, color=blue, /continue
+
+xyouts, .25, 0.3, '+' ,  color=green, /normal, charsize=1.0
+xyouts, .30, 0.3, ' 500 nm (filter 2)', color=black, /normal, charsize=0.8
+xline[0] = 0.75
+yline[0] = 0.3
+xline[1] = 0.80
+yline[1] = 0.3
+plots, xline[0], yline[0], /normal, color=green
+plots, xline[1], yline[1], /normal, color=green, /continue
+
+xyouts, .25, 0.28, '+' ,  color=yellow, /normal, charsize=1.0
+xyouts, .30, 0.28, ' 615 nm (filter 3)', color=black, /normal, charsize=0.8
+xline[0] = 0.75
+yline[0] = 0.28
+xline[1] = 0.80
+yline[1] = 0.28
+plots, xline[0], yline[0], /normal, color=yellow
+plots, xline[1], yline[1], /normal, color=yellow, /continue
+
+xyouts, .25, 0.26, '+' , color=orange, /normal, charsize=1.0
+xyouts, .30, 0.26, ' 673 nm (filter 4)', color=black, /normal, charsize=0.8
+xline[0] = 0.75
+yline[0] = 0.26
+xline[1] = 0.80
+yline[1] = 0.26
+plots, xline[0], yline[0], /normal, color=orange
+plots, xline[1], yline[1], /normal, color=orange, /continue
+
+xyouts, .25, 0.24, '+' , color=red, /normal, charsize=1.0
+xyouts, .30, 0.24, ' 870 nm (filter 5)', color=black, /normal, charsize=0.8
+xline[0] = 0.75
+yline[0] = 0.24
+xline[1] = 0.80
+yline[1] = 0.24
+plots, xline[0], yline[0], /normal, color=red
+plots, xline[1], yline[1], /normal, color=red, /continue
+
+; move "created on" to bottom
+xyouts, 0.1, 0.05, 'Created on: ' + create_date_str, $
+      /normal, charsize=0.7, color=black
+xyouts, 0.1, 0.03,  'Quicklook Version: ' + l_idlVersion, $ 
+      /normal, charsize=0.7, color=black
+xyouts, 0.1, 0.01, 'mfrod1barnmich Version: ' + create_ver_str, $
+      /normal, charsize=0.7, color=black
+
+;=============================================================
+
+end

src/do_io_plot2.pro

+pro do_io_plot2, which_mfr, date, create_ver_str, create_date_str, l_idlVersion, $
+                 l_base_time, l_gauss_times, $
+                 l_Io_gauss_f1, l_Io_gauss_f2, l_Io_gauss_f3, l_Io_gauss_f4, l_Io_gauss_f5
+    
+@colors.include           ; DQ office color table
+
+usersym, d_circle, /fill
+
+Compile_Opt DEFINT32
+
+str_tmp         = string(format='(I0)',date) 
+str_year        = strmid(str_tmp, 0, 4)
+str_month       = strmid(str_tmp, 4, 2)
+str_day         = strmid(str_tmp, 6, 2)
+
+str_date2       = str_year+'/'+str_month+'/'+str_day
+
+; This plot region is only for getting a legible title, since we set the
+; x.charsize in the first real plot to a very small value (so as to not get
+; label on our x-axis ticks.
+
+title_str       = strupcase(which_mfr) + ' Io '+str_date2
+
+;!p.Multi        = [0, 1, 3, 0, 1]
+!p.region       = [0.05, 0.1, 1, 1.]
+!y.charsize     = 3.0
+!x.ticklen      = -0.02
+!y.ticklen      = -0.005
+!y.tickformat   = '(f6.2)'
+
+; We are only showing the gaussian data because we ran the AOD
+; vap using ascii Ios. The Ios we obtain from the ascii file
+; are in the gauss arrays
+;
+
+;print, ' IN DO_IO_PLOT2 ... '
+;print, ' n_elements(l_gauss_time) = ', n_elements(l_gauss_time)
+;print, ' base_time = ', long(l_base_time)
+
+;print, ' n_elements(l_gauss_time) = ', n_elements(l_gauss_time)
+
+xmin = 2500000000000
+xmax = -1
+ymin = 0.95
+ymax = 1.05
+
+found_it = (-1)
+
+;print, ' l_gauss_times has ', n_elements(l_gauss_times), ' elements'
+
+  for n = 0, (n_elements(l_gauss_times)-1) do begin
+    if (l_gauss_times[n] LT xmin) then begin
+      xmin = l_gauss_times[n]
+    endif
+    if (l_gauss_times[n] GT xmax) then begin
+      xmax = l_gauss_times[n]
+    endif 
+;    print, 'n, l_gauss_time[n], l_base_time: ', n, long(l_gauss_times[n]), long(l_base_time)
+ endfor
+
+;print, 'xmin, xmax=', xmin, xmax
+
+  for n = 0, (n_elements(l_gauss_times)-2) do begin
+;    print, 'n, l_gauss_time[n], l_base_time: ', n, long(l_gauss_times[n]), long(l_base_time)
+
+    if (l_gauss_times[n] LE l_base_time) then begin
+;       print, ' n, l_gauss_times[n]=', n, long(l_gauss_times[n])
+;       print, ' n+1, l_gauss_times[n+1], l_base_time=', (n+1), long(l_gauss_times[n+1]), long(l_base_time)
+      if (l_gauss_times[n+1] GT l_base_time) then begin
+        found_it = n
+        break
+      endif
+    endif
+    if (found_it GE 0) then break;
+  endfor
+
+;  print, ' will use gauss time: ', found_it
+
+  for n = 0, (n_elements(l_gauss_times)-1) do begin
+     l_gauss_times[n] = (l_gauss_times[n] - xmin) / 86400 ; converting to days, and normalizing
+;     print, ' NOW: n, l_gauss_times[n] = ', n, l_gauss_times[n]
+  endfor
+  xmax = (xmax - xmin) / 86400
+  xmin = 0
+ 
+;print, 'xmin, xmax=', xmin, xmax
+
+!x.charsize     = 1.5
+!y.charsize     = 1.5
+
+plot,l_gauss_times,             $
+     title      = title_str,    $
+     ystyle     = 1,            $ 
+     xtitle     = 'Time, days', $
+     ytitle     = 'Normalized Io Values', $
+     background = white,        $
+     color      = black,        $
+     xrange     = [xmin, xmax], $
+     yrange     = [ymin, ymax], $
+     xstyle     = 1,            $
+     /NODATA
+
+  oplot, l_gauss_times, l_Io_gauss_f1, color = blue,   linestyle=0 
+  oplot, l_gauss_times, l_Io_gauss_f2, color = green,  linestyle=0
+  oplot, l_gauss_times, l_Io_gauss_f3, color = yellow, linestyle=0
+  oplot, l_gauss_times, l_Io_gauss_f4, color = orange, linestyle=0
+  oplot, l_gauss_times, l_Io_gauss_f5, color = red,    linestyle=0
+
+  aline    = fltarr(2)
+  bline    = fltarr(2)
+  
+; define a vertical line where the "base_time" is, and draw a dashed
+; line there
+
+  aline[0] = l_gauss_times[found_it]
+  aline[1] = l_gauss_times[found_it]
+
+;  print, 'vert line at ', aline[0]
+
+  bline[0] = ymin
+  bline[1] = ymax
+
+  oplot, aline[0:1], bline[0:1], color = black, linestyle=1
+  xyouts, aline[0], (ymin + 0.01), str_date2, orientation=90., color=black, charsize=0.8
+
+ ; generate a horizontal line  at the Io=1.0 location
+
+  aline[0] = l_gauss_times[found_it]
+  aline[1] = l_gauss_times[found_it]
+  bline[0] = ymin
+  bline[1] = ymax
+
+  oplot, aline[0:1], bline[0:1], color = black, linestyle=1
+
+; legend info copied from do_od_plot2.pro
+;=============================================================
+
+; annotate all plots in the plot at the bottom
+
+;!p.region = [0., 0., 1., .1]
+!x.charsize     = .5
+!y.charsize     = .5
+
+xline    = fltarr(2)
+yline    = fltarr(2)
+
+xyouts, .32,  0.345,'Io', color=black, /normal, charsize=1.0
+
+xyouts, .30, 0.32, '  415 nm (filter 1)', color=black, /normal, charsize=0.8
+xline[0] = 0.25
+yline[0] = 0.32
+xline[1] = 0.30
+yline[1] = 0.32
+plots, xline[0], yline[0], /normal, color=blue
+plots, xline[1], yline[1], /normal, color=blue, /continue
+
+xyouts, .30, 0.3, '  500 nm (filter 2)', color=black, /normal, charsize=0.8
+xline[0] = 0.25
+yline[0] = 0.3
+xline[1] = 0.30
+yline[1] = 0.3
+plots, xline[0], yline[0], /normal, color=green
+plots, xline[1], yline[1], /normal, color=green, /continue
+
+xyouts, .30, 0.28, '  615 nm (filter 3)', color=black, /normal, charsize=0.8
+xline[0] = 0.25
+yline[0] = 0.28
+xline[1] = 0.30
+yline[1] = 0.28
+plots, xline[0], yline[0], /normal, color=yellow
+plots, xline[1], yline[1], /normal, color=yellow, /continue
+
+xyouts, .30, 0.26, '  673 nm (filter 4)', color=black, /normal, charsize=0.8
+xline[0] = 0.25
+yline[0] = 0.26
+xline[1] = 0.30
+yline[1] = 0.26
+plots, xline[0], yline[0], /normal, color=orange
+plots, xline[1], yline[1], /normal, color=orange, /continue
+
+xyouts, .30, 0.24, '  870 nm (filter 5)', color=black, /normal, charsize=0.8
+xline[0] = 0.25
+yline[0] = 0.24
+xline[1] = 0.30
+yline[1] = 0.24
+plots, xline[0], yline[0], /normal, color=red
+plots, xline[1], yline[1], /normal, color=red, /continue
+
+; move "created on" to bottom
+xyouts, 0.1, 0.05, 'Created on: ' + create_date_str, $
+      /normal, charsize=0.7, color=black
+xyouts, 0.1, 0.03,  'Quicklook Version: ' + l_idlVersion, $ 
+      /normal, charsize=0.7, color=black
+xyouts, 0.1, 0.01, 'mfrod1barnmich Version: ' + create_ver_str, $
+      /normal, charsize=0.7, color=black
+
+;=============================================================
+
+end

src/do_od_plot1.pro

+pro do_od_plot1, which_mfr, $
+                 date, time_samples, create_ver_str, create_date_str, l_idlVersion, $
+                 Langley_used_str, data_used_str, $
+                 l_aod_f2,       $
+                 l_qc_aod_f2,    $
+                 l_aod_f5,       $
+                 l_qc_aod_f5,    $
+;                 sunrise, sunset, $
+	         l_angstrom, l_qc_angstrom, $
+                 l_lat, l_lon
+
+;@color_syms.include   Our color table
+
+@colors.include       ; DQ office color table
+
+usersym, d_circle, /fill
+;dloadct, 42
+
+Compile_Opt DEFINT32
+
+str_tmp         = string(format='(I0)',date)
+str_year        = strmid(str_tmp, 0, 4)
+str_month       = strmid(str_tmp, 4, 2)
+str_day         = strmid(str_tmp, 6, 2)
+
+str_date2       = str_year+'/'+str_month+'/'+str_day
+
+!P.Multi        = [0, 1, 4, 0, 1]  ; erase the page, 1 column of 5 plots, top to bottom
+
+; This plot region is only for getting a legible title, since we set the
+; x.charsize in the first real plot to a very small value (so as to not get
+; label on our x-axis ticks.
+
+!p.region       = [0.0, .9, 1., 1.]
+
+title_str       = strupcase(which_mfr) + ' Aerosol Optical Depths for '+str_date2
+
+!p.region       = [0.1, .5, 1, .95]
+!x.charsize     = .001
+!y.charsize     = 3.0
+!x.ticklen      = -0.02
+!y.ticklen      = -0.005
+!y.tickformat   = '(f6.1)'
+
+; First real plot (aerosol_optical_depth filter 2)
+
+npts = n_elements(l_aod_f2)
+w_qc = intarr(npts)
+w_qc = l_aod_f2
+
+bits    = intarr(10)
+qc_bits = intarr(10)
+qc_bits = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512]
+
+allbad  = 0L
+w_value = 0L
+
+;
+; Connor wants to plot points as yellow if qc is indeterminate
+; and to plot points as green if qc is good
+; not plotting obviously bad points
+;
+; for AOD, QC bit 6 to 10 are an indeterminate bit
+;              
+
+curve_green                = l_aod_f2
+curve_green[0:npts-1]      = 0
+
+sum_od_green_f2            = 0
+sum_od_green_f5            = 0
+
+time_green                 = time_samples
+time_green[0:npts-1]       = 0
+npts_green                 = -1
+
+ymax                       = 0.0
+
+for n  = 1, npts-1 do begin
+;
+  w_value = l_qc_aod_f2[n]
+
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[0:5] are set, the value is indeterminate
+;
+; if ( (bits[5] > 0 OR bits[6] > 0 OR bits[7] > 0 OR bits[8] > 0 OR bits[9] > 0) ) then continue
+  
+  if ( (bits[0] > 0 OR bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0 ) AND (l_aod_f2[n] > (-9000.0)) ) then begin
+  endif else begin
+    npts_green                = npts_green + 1
+    time_green[npts_green]    = time_samples[n]
+    curve_green[npts_green]   = l_aod_f2[n]
+
+    sum_od_green_f2           = sum_od_green_f2 + l_aod_f2[n]
+
+    if (curve_green[npts_green] GT ymax) then begin
+      ymax = curve_green[npts_green]
+    endif
+  endelse
+endfor
+
+xmin = time_samples[0] - .1
+xmax = time_samples[n_elements(l_aod_f2)-1] + .1
+ymax = 1.1*ymax
+
+; xxxxxxxxxxxxxxxxxxxxxxxxx
+;
+;  similar logic for filter 5 (add to top plot)
+;
+
+curve_green5                = l_aod_f5
+npts5                       = n_elements(l_aod_f5)
+curve_green5[0:npts5-1]     = 0
+time_green5                 = time_samples
+time_green5[0:npts5-1]      = 0
+npts_green5                 = -1
+
+for n  = 1, npts5-1 do begin
+;
+  w_value = l_qc_aod_f5[n]
+
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[0:5] are set, the value is indeterminate
+;
+;  if ( (bits[5] > 0 OR bits[6] > 0 OR bits[7] > 0 OR bits[8] > 0 OR bits[9] > 0) ) then continue
+  
+  if ( (bits[0] > 0 OR bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0 ) AND (l_aod_f5[n] > (-9000.0)) ) then begin
+  endif else begin
+    npts_green5                = npts_green5 + 1
+    time_green5[npts_green5]   = time_samples[n]
+    curve_green5[npts_green5]  = l_aod_f5[n]
+
+    sum_od_green_f5            = sum_od_green_f5 + l_aod_f5[n]
+
+  endelse
+endfor
+
+avg_od_green_f2   = 0
+avg_od_green_f5   = 0
+
+if (npts_green GT 0) then begin
+  avg_od_green_f2 = sum_od_green_f2 / npts_green
+endif
+
+if (npts_green5 GT 0) then begin
+  avg_od_green_f5 = sum_od_green_f5 / npts_green5
+endif
+
+averages    = fltarr(2)
+averages[0] = avg_od_green_f2
+averages[1] = avg_od_green_f5
+
+ymax = max(averages) * 1.8
+
+; xxxxxxxxxxxxxxxxxxxxxxxxx
+
+do_the_plot = 0L
+do_the_plot = npts_green + npts_green5
+
+if (do_the_plot LE 0 ) then begin
+;
+; No points to plot, so draw a box, 
+; but no axes
+;
+
+!y.charsize = 0.001
+
+; get time axis in hh:mm form
+
+plot, time_samples, l_aod_f2, $
+      ystyle     = 1,  $
+      ytickformat= '(f5.2)', $
+      ytitle     = 'Aerosol Optical Depth', $
+      background = white, $
+      color      = black, $
+      xrange     = [xmin, xmax], $
+      xstyle     = 1, $
+      XTICKFORMAT='LABEL_DATE', XTICKS=8, XTICKUNITS=['Time', 'Time'], $
+      /NODATA
+
+xyouts, 0.35, 0.8, $
+      'No valid aerosol optical depths for this day.', $
+       /normal, charsize=1.5, color=red
+
+endif else begin
+
+  ; We have something to plot, so draw axes,
+  ; then add the curve
+
+  ; get time axis in hh:mm form
+
+  dummy      = LABEL_DATE(DATE_FORMAT=['%H:%I'])
+
+  time_hours = TIMEGEN(START=time_samples[0], FINAL=(time_samples[0]+1.), UNITS='Hours')
+
+  plot, time_hours, l_aod_f2, $
+      ystyle = 1, $
+      ytickformat= '(f5.2)', $
+      ytitle = 'Aerosol Optical Depth', $
+      background = white, $
+      color      = black, $
+      yrange     = [0, ymax], $
+      xstyle = 1, $
+      XTICKFORMAT='LABEL_DATE', XTICKS=8, XTICKUNITS=['Hours'], xcharsize=0.01, $
+      /NODATA
+
+  ; Try to add the day/night background
+
+  day_night_background, l_lat, l_lon, /NOON
+
+  ; redraw the axes ...
+
+  axis,color=black,charsize=0.01, xticklen=0.0, $
+      XTICKFORMAT='LABEL_DATE', XTICKS=8, XTICKUNITS=['Hours']
+  
+  axis,color=black,yaxis=0,yrange=[0,ymax],ystyle=1,charsize=0.01,xtickunits='Time'
+  axis,color=black,yaxis=1,yrange=[0,ymax],ystyle=1,charsize=0.01,xtickunits='Time'
+
+  if (npts_green GT 0) then begin
+    oplot, time_green[0:npts_green-1], curve_green[0:npts_green-1], color = green, $
+       psym=1, symsize=1.5 
+  endif
+
+  if (npts_green5 GT 0) then begin
+    oplot, time_green5[0:npts_green5-1], curve_green5[0:npts_green5-1], color = red, $
+       psym=1, symsize=1.5
+  endif
+
+
+; next plot (angstrom exponent)
+
+!p.region = [0.1, 0.15, 1.0, 0.45]
+!y.ticks = 2
+!y.charsize = 3
+!x.charsize = 3
+
+curve           = l_angstrom
+num_pts         = n_elements(curve)
+start_line      = 0
+max_time_diff   = 0.0055
+
+curve_green[0:npts-1]      = 0
+time_green                 = time_samples
+time_green[0:npts-1]       = 0
+npts_green                 = -1
+
+for n  = 1, npts-1 do begin
+;
+  w_value = l_qc_angstrom[n]
+
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if bits[0] is set, the value of angstrom is indeterminate
+;
+  if ( (bits[1] > 0 OR bits[2] > 0 OR bits[3] ) ) then continue 
+  
+  if ( (bits[0] > 0 ) AND (l_angstrom[n] > (-9000.0)) ) then begin
+  endif else begin
+    npts_green                = npts_green + 1
+    time_green[npts_green]    = time_samples[n]
+    curve_green[npts_green]   = l_angstrom[n]
+  endelse
+endfor
+
+plot, time_hours, curve, ystyle=1, $
+      ytitle = 'Angstrom Exponent', $
+      xtitle='Time of Day, GMT', $
+      yrange=[-0.5, 2.0], $
+      background=white, color=black, $
+      xstyle = 1, $
+      XTICKFORMAT='LABEL_DATE', XTICKS=8, XTICKUNITS=['Hours'], $
+      /NODATA
+
+; Try to add the day/night background
+
+day_night_background, l_lat, l_lon, /NOON
+
+  ; redraw the axes ...
+
+axis,color=black,charsize=0.01, $
+      XTICKFORMAT='LABEL_DATE', XTICKS=8, XTICKUNITS=['Hours']
+
+axis,color=black,yaxis=0,yrange=[-0.5,2.0],ystyle=1,xtickunits='Time'
+axis,color=black,yaxis=1,yrange=[-0.5,2.0],ystyle=1,charsize=0.01,xtickunits='Time'
+
+if (npts_green GT 0) then begin
+  oplot, time_green, curve_green, color=black, $ 
+       psym=7, symsize=1.5
+endif
+
+endelse
+
+; annotate all plots in the plot at the bottom
+
+!p.region = [0., 0., 1., .1]
+
+;xyouts, 0.5, 0.06, $
+;      'Optical depths calculated from ' + data_used_str, $
+;       /normal, charsize=1., color=black
+xyouts, 0.5, 0.04, $
+      'Created on: ' + create_date_str, $
+      /normal, charsize=1., color=black
+xyouts, 0.5, 0.02, $
+      'Quicklook Version: ' + l_idlVersion, $ 
+      /normal, charsize=1., color=black
+xyouts, 0.5, 0.00, $
+      'mfrod1barnmich Version: ' + create_ver_str, $
+      /normal, charsize=1., color=black
+
+xyouts, 350, 730, title_str, color=black, /device, charsize=1.5, charthick=1.25
+
+xyouts, .08,  0.06, 'Nominal Wavelengths', color=black, /normal, charsize=1.0
+
+xyouts, .08, 0.04, '+' , $
+        color=red, /normal, charsize=1.0
+xyouts, .10, 0.04, ' = 500 nm (filter 2)', $
+        color=black, /normal, charsize=1.0
+ 
+xyouts, .08, 0.02, '+', $
+        color=green, /normal, charsize=1.0
+xyouts, .10, 0.02, ' = 870 nm (filter 5)', $
+        color=black, /normal, charsize=1.0
+ 
+
+
+end

src/do_od_plot2.pro

+pro do_od_plot2, which_mfr, $
+                 date, time_samples, create_ver_str, create_date_str, l_idlVersion, $
+                 Langley_used_str, data_used_str, $
+                 l_aod_f1,       l_aod_f2,       l_aod_f3,       l_aod_f4,       l_aod_f5, $
+                 l_qc_aod_f1,    l_qc_aod_f2,    l_qc_aod_f3,    l_qc_aod_f4,    l_qc_aod_f5, $
+;                 sunrise, sunset, $
+	         l_angstrom, l_qc_angstrom, $
+                 l_lat,    l_lon
+
+; @color_syms.include  our color table
+
+@colors.include     ; DQ office color_table
+
+usersym, d_circle, /fill
+;dloadct, 42
+
+Compile_Opt DEFINT32
+
+str_tmp         = string(format='(I0)',date)
+str_year        = strmid(str_tmp, 0, 4)
+str_month       = strmid(str_tmp, 4, 2)
+str_day         = strmid(str_tmp, 6, 2)
+
+str_date2       = str_year+'/'+str_month+'/'+str_day
+
+!P.Multi        = [0, 1, 3, 0, 1]  ; erase the page, 1 column of 5 plots, top to bottom
+
+; This plot region is only for getting a legible title, since we set the
+; x.charsize in the first real plot to a very small value (so as to not get
+; label on our x-axis ticks.
+
+title_str       = strupcase(which_mfr) + ' Aerosol Optical Depths for '+str_date2
+
+!p.region       = [0.1, .2, 1, .95]
+!x.charsize     = 3.0
+!y.charsize     = 3.0
+!x.ticklen      = -0.02
+!y.ticklen      = -0.005
+!y.tickformat   = '(f6.1)'
+
+; First real plot (aerosol_optical_depth filter 2)
+
+npts = n_elements(l_aod_f2)
+w_qc = intarr(npts)
+w_qc = l_aod_f2
+
+bits    = intarr(10)
+qc_bits = intarr(10)
+qc_bits = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512]
+
+w_value = 0L
+
+;
+; Connor wants to plot points as yellow if qc is indeterminate
+; and to plot points as green if qc is good
+; not plotting obviously bad points
+;
+; for AOD, QC bit 6 to 10 are an indeterminate bit
+;              
+
+curve_f1                   = l_aod_f2
+curve_f2                   = l_aod_f2
+curve_f3                   = l_aod_f2
+curve_f4                   = l_aod_f2
+curve_f5                   = l_aod_f2
+
+time_f1                    = time_samples
+time_f2                    = time_samples
+time_f3                    = time_samples
+time_f4                    = time_samples
+time_f5                    = time_samples
+
+time_f1[0:npts-1]          = 0
+time_f2[0:npts-1]          = 0
+time_f3[0:npts-1]          = 0
+time_f4[0:npts-1]          = 0
+time_f5[0:npts-1]          = 0
+npts_f1                    = -1
+npts_f2                    = -1
+npts_f3                    = -1
+npts_f4                    = -1
+npts_f5                    = -1
+
+sum_od_f1                  = 0
+sum_od_f2                  = 0
+sum_od_f3                  = 0
+sum_od_f4                  = 0
+sum_od_f5                  = 0
+
+ymax = 0.0
+
+for n  = 1, npts-1 do begin
+  w_value = l_qc_aod_f1[n]
+
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+
+;  if ( (bits[5] > 0 OR bits[6] > 0 OR bits[7] > 0 OR bits[8] > 0 OR bits[9] > 0) ) then continue
+  if ( (bits[0] > 0 OR bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0) AND (l_aod_f1[n] > (-9000.0)) ) then begin
+  endif else begin
+;
+;   a bit is set, but the value is > 0, so we will plot it
+;
+    npts_f1           = npts_f1 + 1
+    time_f1[npts_f1]  = time_samples[n]
+    curve_f1[npts_f1] = l_aod_f1[n]
+
+    sum_od_f1         = sum_od_f1 + l_aod_f1[n]
+
+;      print, 'n=', n, ' AOD = ', l_aod_f1[n], ' sum_od_f1=', sum_od_f1
+
+  endelse
+
+endfor
+
+for n  = 1, npts-1 do begin
+  w_value = l_qc_aod_f2[n]
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+
+;  if ( (bits[5] > 0 OR bits[6] > 0 OR bits[7] > 0 OR bits[8] > 0 OR bits[9] > 0) ) then continue
+;  if ( (bits[0] > 0 OR bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR
+;  bits[4] > 0) AND (l_aod_f2[n] > (-9000.0)) ) then begin
+  if ( (bits[0] > 0 OR bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0) AND (l_aod_f2[n] >= (-9000.0)) ) then begin
+  endif else begin
+;
+;   a bit is set, but the value is > 0, so we will plot it
+;
+    npts_f2           = npts_f2 + 1
+    time_f2[npts_f2]  = time_samples[n]
+    curve_f2[npts_f2] = l_aod_f2[n]
+    sum_od_f2         = sum_od_f2 + l_aod_f2[n]
+  endelse
+endfor
+
+for n  = 1, npts-1 do begin
+  w_value = l_qc_aod_f3[n]
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+
+;  if ( (bits[5] > 0 OR bits[6] > 0 OR bits[7] > 0 OR bits[8] > 0 OR bits[9] > 0) ) then continue
+  if ( (bits[0] > 0 OR bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0) AND (l_aod_f3[n] > (-9000.0)) ) then begin
+  endif else begin
+;
+;   a bit is set, but the value is > 0, so we will plot it
+;
+    npts_f3           = npts_f3 + 1
+    time_f3[npts_f3]  = time_samples[n]
+    curve_f3[npts_f3] = l_aod_f3[n]
+    sum_od_f3         = sum_od_f3 + l_aod_f3[n]
+
+  endelse
+endfor
+
+for n  = 1, npts-1 do begin
+  w_value = l_qc_aod_f4[n]
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+
+;  if ( (bits[5] > 0 OR bits[6] > 0 OR bits[7] > 0 OR bits[8] > 0 OR bits[9] > 0) ) then continue
+  if ( (bits[0] > 0 OR bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0) AND (l_aod_f4[n] > (-9000.0)) ) then begin
+  endif else begin
+;
+;   a bit is set, but the value is > 0, so we will plot it
+;
+    npts_f4           = npts_f4 + 1
+    time_f4[npts_f4]  = time_samples[n]
+    curve_f4[npts_f4] = l_aod_f4[n]
+    sum_od_f4         = sum_od_f4 + l_aod_f4[n]
+  endelse
+endfor
+
+for n  = 1, npts-1 do begin
+  w_value = l_qc_aod_f5[n]
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+
+;  if ( (bits[5] > 0 OR bits[6] > 0 OR bits[7] > 0 OR bits[8] > 0 OR bits[9] > 0) ) then continue
+  if ( (bits[0] > 0 OR bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0) AND (l_aod_f5[n] > (-9000.0)) ) then begin
+  endif else begin
+;
+;   a bit is set, but the value is > 0, so we will plot it
+;
+    npts_f5           = npts_f5 + 1
+    time_f5[npts_f5]  = time_samples[n]
+    curve_f5[npts_f5] = l_aod_f5[n]
+    sum_od_f5         = sum_od_f5 + l_aod_f5[n]
+  endelse
+endfor
+
+xmin = time_samples[0] - .1
+xmax = time_samples[n_elements(l_aod_f2)-1] + .1
+
+averages      = fltarr(5)
+averages[0:4] = 0
+
+;  print, ' sum_od_f1 = ', sum_od_f1, ' npts_f1=', npts_f1
+
+if (npts_f1 GT 0) then begin
+  averages[0] = sum_od_f1 / npts_f1
+endif
+if (npts_f2 GT 0) then begin
+  averages[1] = sum_od_f2 / npts_f2
+endif
+if (npts_f3 GT 0) then begin
+  averages[2] = sum_od_f3 / npts_f3
+endif
+if (npts_f4 GT 0) then begin
+  averages[3] = sum_od_f4 / npts_f4
+endif
+if (npts_f5 GT 0) then begin
+  averages[4] = sum_od_f5 / npts_f5
+endif
+
+;print, ' averages=', averages[0], averages[1],averages[2],averages[3],averages[4]
+
+ymax = max(averages)
+
+ymax = ymax*1.5
+
+do_the_plot = 0L
+do_the_plot = npts_f1 + npts_f2 + npts_f3 + npts_f4 + npts_f5 
+
+if (do_the_plot LE 0) then begin
+
+  !y.charsize = 0.001
+
+  plot, time_samples, l_aod_f2, $
+        ystyle = 1, xstyle = 1, $
+        ytickformat='(f5.2)', $
+        ytitle = 'Aerosol Optical Depths', $
+        xtitle = 'Julian Day', $
+        background = white, $
+        color      = black, $
+        yrange     = [0, ymax], $
+        xrange     = [xmin, xmax], $
+        /NODATA
+
+  xyouts, 0.35, 0.8, $
+        'No valid aerosol optical depths for this day.', $
+         /normal, charsize=1.5, color=red
+
+endif else begin
+
+; DQ office prefers time axis in hh:mm form
+
+  dummy = LABEL_DATE(DATE_FORMAT=['%H:%I'])
+
+  time_hours = TIMEGEN(START=time_samples[0], FINAL=(time_samples[0]+1.), UNITS='Hours')
+
+  n_time_hours = n_elements(time_hours)
+
+  plot, time_hours, l_aod_f2, $
+        ystyle = 1, $
+        ytickformat='(f5.2)', $
+        ytitle = 'Aerosol Optical Depths', $
+        xtitle = 'Time of Day, GMT', $
+        background = white, $
+        color      = black, $
+        yrange     = [0, ymax], $
+        xstyle = 1, $
+        XTICKFORMAT='LABEL_DATE', XTICKS=8, XTICKUNITS=['Hours'], $
+;        xrange     = [xmin, xmax], $
+        /NODATA
+
+  ; Try to add the day/night background
+
+  day_night_background, l_lat, l_lon, /NOON
+
+  ; redraw the axes ...
+
+  axis,color=black,charsize=0.01, xticks=1, xticklen=0.0, $
+        XTICKFORMAT='LABEL_DATE', XTICKUNITS=['Hours']
+
+  axis,color=black,yaxis=0,yrange=[0,ymax],ystyle=1,charsize=0.01, $
+        XTICKFORMAT='LABEL_DATE', XTICKS=8, XTICKUNITS=['Hours']
+  axis,color=black,yaxis=1,yrange=[0,ymax],ystyle=1,charsize=0.01,xtickunits='Time'
+
+  if (npts_f1 GT 0) then begin
+    oplot, time_f1[0:npts_f1-1], curve_f1[0:npts_f1-1], color = blue, $
+       psym=1, symsize=1.5
+  endif
+  if (npts_f2 GT 0) then begin
+    oplot, time_f2[0:npts_f2-1], curve_f2[0:npts_f2-1], color = green, $
+       psym=1, symsize=1.5
+  endif
+  if (npts_f3 GT 0) then begin
+    oplot, time_f3[0:npts_f3-1], curve_f3[0:npts_f3-1], color = yellow, $
+       psym=1, symsize=1.5
+  endif
+  if (npts_f4 GT 0) then begin
+    oplot, time_f4[0:npts_f4-1], curve_f4[0:npts_f4-1], color = orange, $
+       psym=1, symsize=1.5
+  endif
+  if (npts_f5 GT 0) then begin
+    oplot, time_f5[0:npts_f5-1], curve_f5[0:npts_f5-1], color = red, $
+       psym=1, symsize=1.5
+  endif
+
+
+ 
+endelse
+
+; annotate all plots in the plot at the bottom
+
+!p.region = [0., 0., 1., .1]
+
+;xyouts, 0.4, 0.075, $
+;      'Optical depths calculated from ' + data_used_str, $
+;       /normal, charsize=1.0, color=black
+xyouts, 0.4, 0.055, $
+      'Created on: ' + create_date_str, $
+      /normal, charsize=1.0, color=black
+xyouts, 0.4, 0.035, $
+      'Quicklook Version: ' + l_idlVersion, $ 
+      /normal, charsize=1.0, color=black
+xyouts, 0.4, 0.015, $
+      'mfrod1barnmich Version: ' + create_ver_str, $
+      /normal, charsize=1.0, color=black
+
+xyouts, 350, 730, title_str, color=black, /device, charsize=1.5, charthick=1.25
+
+
+  xyouts, .10,  0.125,'Nominal Wavelengths', color=black, /normal, charsize=1.0
+
+  xyouts, .08, 0.095, '+' , $
+         color=blue, /normal, charsize=1.0
+  xyouts, .10, 0.095, '= 415 nm (filter 1)', $
+          color=black, /normal, charsize=1.0
+ 
+  xyouts, .08, 0.075, '+' , $
+         color=green, /normal, charsize=1.0
+  xyouts, .10, 0.075, '= 500 nm (filter 2)', $
+          color=black, /normal, charsize=1.0
+ 
+  xyouts, .08, 0.055, '+' , $
+         color=yellow, /normal, charsize=1.0
+  xyouts, .10, 0.055, '= 615 nm (filter 3)', $
+          color=black, /normal, charsize=1.0
+ 
+  xyouts, .08, 0.035, '+' , $
+         color=orange, /normal, charsize=1.0
+  xyouts, .10, 0.035,  '= 673 nm (filter 4)', $
+          color=black, /normal, charsize=1.0
+ 
+  xyouts, .08, 0.015, '+' , $
+         color=red, /normal, charsize=1.0
+  xyouts, .10, 0.015, '= 870 nm (filter 5)', $
+          color=black, /normal, charsize=1.0
+ 
+
+end

src/do_od_plot3.pro

+pro do_od_plot3, which_mfr, $ 
+                 date, time_samples, create_ver_str, create_date_str, l_idlVersion, $
+                 Langley_used_str, data_used_str, $
+                 l_dirnor_f1,     l_dirnor_f2,     l_dirnor_f3,     l_dirnor_f4,     l_dirnor_f5, $
+                 l_qc_dirnor_f1,  l_qc_dirnor_f2,  l_qc_dirnor_f3,  l_qc_dirnor_f4,  l_qc_dirnor_f5, $
+                 l_diffuse_f1,    l_diffuse_f2,    l_diffuse_f3,    l_diffuse_f4,    l_diffuse_f5, $
+                 l_qc_diffuse_f1, l_qc_diffuse_f2, l_qc_diffuse_f3, l_qc_diffuse_f4, l_qc_diffuse_f5, $
+;                 sunrise, sunset, $
+	         l_angstrom, l_qc_angstrom, $
+                 l_lat, l_lon
+; @color_syms.include     ; our color table
+
+@colors.include           ; DQ office color table
+
+usersym, d_circle, /fill
+;dloadct, 42
+
+
+Compile_Opt DEFINT32
+
+str_tmp         = string(format='(I0)',date)
+str_year        = strmid(str_tmp, 0, 4)
+str_month       = strmid(str_tmp, 4, 2)
+str_day         = strmid(str_tmp, 6, 2)
+
+str_date2       = str_year+'/'+str_month+'/'+str_day
+
+!P.Multi        = [0, 1, 4, 0, 1]  ; erase the page, 1 column of 5 plots, top to bottom
+
+; This plot region is only for getting a legible title, since we set the
+; x.charsize in the first real plot to a very small value (so as to not get
+; label on our x-axis ticks.
+
+!p.region       = [0.1, .5, 1, .95]
+
+title_str       = strupcase(which_mfr) + ' Direct Normal and Diffuse Hemispheric '+str_date2
+if (which_mfr EQ 'nimfr') then begin
+   title_str       = strupcase(which_mfr) + ' Direct Normal '+str_date2
+  !p.region       = [0.1, .15, 1., .95]
+endif
+
+!y.charsize     = 3.0
+!x.ticklen      = -0.02
+!y.ticklen      = -0.005
+!y.tickformat   = '(f6.1)'
+
+; First real plot (aerosol_optical_depth filter 2)
+
+npts = n_elements(l_dirnor_f2)
+w_qc = intarr(npts)
+w_qc = l_dirnor_f2
+
+bits    = intarr(10)
+qc_bits = intarr(10)
+qc_bits = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512]
+
+allbad  = 0L
+w_value = 0L
+
+;
+; Connor wants to plot points as yellow if qc is indeterminate
+; and to plot points as green if qc is good
+; not plotting obviously bad points
+;
+; for AOD, QC bit 6 to 10 are an indeterminate bit
+;              
+
+curve_f1               = l_dirnor_f2
+curve_f2               = l_dirnor_f2
+curve_f3               = l_dirnor_f2
+curve_f4               = l_dirnor_f2
+curve_f5               = l_dirnor_f2
+curve_f1[0:npts-1]     = 0
+curve_f2[0:npts-1]     = 0
+curve_f3[0:npts-1]     = 0
+curve_f4[0:npts-1]     = 0
+curve_f5[0:npts-1]     = 0
+time_f1                = time_samples
+time_f2                = time_samples
+time_f3                = time_samples
+time_f4                = time_samples
+time_f5                = time_samples
+time_f1[0:npts-1]      = 0
+time_f2[0:npts-1]      = 0
+time_f3[0:npts-1]      = 0
+time_f4[0:npts-1]      = 0
+time_f5[0:npts-1]      = 0
+npts_f1                = -1
+npts_f2                = -1
+npts_f3                = -1
+npts_f4                = -1
+npts_f5                = -1
+
+sum_f1                 = 0
+sum_f2                 = 0
+sum_f3                 = 0
+sum_f4                 = 0
+sum_f5                 = 0
+
+for n  = 1, npts-1 do begin
+  w_value = l_qc_dirnor_f1[n]
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+
+  if ( (bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0 OR bits[5] > 0 ) ) then continue
+  if ( bits[0] > 0 ) then begin
+    if (l_dirnor_f1[n] >= 0) then begin
+      npts_f1           = npts_f1 + 1
+      time_f1[npts_f1]  = time_samples[n]
+      curve_f1[npts_f1] = l_dirnor_f1[n]
+
+      sum_f1            = sum_f1 + l_dirnor_f1[n]
+
+;      print, 'AT A npts_f1=',npts_f1, ' n=', n, ' l_dirnor_f1[n] = ', l_dirnor_f1[n], ' sum_f1=', sum_f1
+    endif
+
+  endif else begin
+    if (l_dirnor_f1[n] >= 0.0) then begin
+      npts_f1           = npts_f1 + 1
+      time_f1[npts_f1]  = time_samples[n]
+      curve_f1[npts_f1] = l_dirnor_f1[n]
+
+      sum_f1            = sum_f1 + l_dirnor_f1[n]
+  
+;      print, 'AT B npts_f1=',npts_f1, ' n=', n, ' l_dirnor_f1[n] = ', l_dirnor_f1[n], ' sum_f1=', sum_f1
+
+    endif
+  endelse
+endfor
+
+for n  = 1, npts-1 do begin
+  w_value = l_qc_dirnor_f2[n]
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+
+  if ( (bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0 OR bits[5] > 0 ) ) then continue
+  if ( bits[0] > 0 ) then begin
+    if (l_dirnor_f2[n] >= 0) then begin
+      npts_f2           = npts_f2 + 1
+      time_f2[npts_f2]  = time_samples[n]
+      curve_f2[npts_f2] = l_dirnor_f2[n]
+      sum_f2            = sum_f2 + l_dirnor_f2[n]
+    endif
+  endif else begin
+    if (l_dirnor_f2[n] > (-9000.0)) then begin
+      npts_f2           = npts_f2 + 1
+      time_f2[npts_f2]  = time_samples[n]
+      curve_f2[npts_f2] = l_dirnor_f2[n]
+      sum_f2            = sum_f2 + l_dirnor_f2[n]
+    endif
+  endelse
+endfor
+
+for n  = 1, npts-1 do begin
+  w_value = l_qc_dirnor_f3[n]
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+
+  if ( (bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0 OR bits[5] > 0 ) ) then continue
+  if ( bits[0] > 0 ) then begin
+    if (l_dirnor_f3[n] >= 0) then begin
+      npts_f3           = npts_f3 + 1
+      time_f3[npts_f3]  = time_samples[n]
+      curve_f3[npts_f3] = l_dirnor_f3[n]
+      sum_f3            = sum_f3 + l_dirnor_f3[n]
+    endif
+  endif else begin
+    if (l_dirnor_f3[n] >= 0) then begin
+      npts_f3           = npts_f3 + 1
+      time_f3[npts_f3]  = time_samples[n]
+      curve_f3[npts_f3] = l_dirnor_f3[n]
+      sum_f3            = sum_f3 + l_dirnor_f3[n]
+    endif
+  endelse
+endfor
+
+for n  = 1, npts-1 do begin
+  w_value = l_qc_dirnor_f4[n]
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+
+  if ( (bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0 OR bits[5] > 0 ) ) then continue
+  if ( bits[0] > 0 ) then begin
+    if (l_dirnor_f4[n] >= 0.0) then begin
+      npts_f4           = npts_f4 + 1
+      time_f4[npts_f4]  = time_samples[n]
+      curve_f4[npts_f4] = l_dirnor_f4[n]
+      sum_f4            = sum_f4 + l_dirnor_f4[n]
+    endif
+  endif else begin
+    if (l_dirnor_f4[n] >= 0.0) then begin
+      npts_f4           = npts_f4 + 1
+      time_f4[npts_f4]  = time_samples[n]
+      curve_f4[npts_f4] = l_dirnor_f4[n]
+      sum_f4            = sum_f4 + l_dirnor_f4[n]
+    endif
+  endelse
+endfor
+
+for n  = 1, npts-1 do begin
+  w_value = l_qc_dirnor_f5[n]
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+
+  if ( (bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0 OR bits[5] > 0 ) ) then continue
+  if ( bits[0] > 0 ) then begin
+    if (l_dirnor_f4[n] >= 0.0) then begin
+      npts_f5           = npts_f5 + 1
+      time_f5[npts_f5]  = time_samples[n]
+      curve_f5[npts_f5] = l_dirnor_f5[n]
+      sum_f5            = sum_f5 + l_dirnor_f5[n]
+    endif
+  endif else begin
+    if (l_dirnor_f4[n] >= 0.0) then begin
+      npts_f5           = npts_f5 + 1
+      time_f5[npts_f5]  = time_samples[n]
+      curve_f5[npts_f5] = l_dirnor_f5[n]
+      sum_f5            = sum_f5 + l_dirnor_f5[n]
+    endif
+  endelse
+endfor
+
+;max_y           = 2.0
+
+averages         = fltarr(5)
+averages[0:4]    = 0
+
+if (npts_f1 GT 0) then begin
+  averages[0] = sum_f1 / npts_f1
+;  print, 'sum_f1=', sum_f1, ' npts_f1=', npts_f1
+endif
+if (npts_f2 GT 0) then begin
+  averages[1] = sum_f2 / npts_f2
+endif
+if (npts_f3 GT 0) then begin
+  averages[2] = sum_f3 / npts_f3
+endif
+if (npts_f4 GT 0) then begin
+  averages[3] = sum_f4 / npts_f4
+endif
+if (npts_f5 GT 0) then begin
+  averages[4] = sum_f5 / npts_f5
+endif
+
+;print, 'averages[0]=', averages[0]
+
+ymax = max(averages)
+
+;print, 'IRRAD plot: ymax=', ymax
+
+;xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+if (ymax LE 0.0) then begin
+                                ;
+                                ; No points to plot, so draw a box, 
+                                ; but no axes
+                                ;
+
+   !y.charsize = 0.001
+
+   dummy      = LABEL_DATE(DATE_FORMAT=['%H:%I'])
+
+   time_hours = TIMEGEN(START=time_samples[0], FINAL=(time_samples[0]+1.), UNITS='Hours')
+
+   ;print, ' NO IRRADIANCE DATA ... skipping plot, averages, ymax=', averages, ymax
+
+   plot, time_hours, l_dirnor_f5, $
+         ystyle = 1, $ 
+         xtitle = xtitl, $
+         ytitle = 'Direct Normal', $
+         background = white, $
+         color      = black, $
+         yrange     = [0.0, ymax], $
+         xstyle = 1, $
+         XTICKFORMAT='LABEL_DATE', XTICKS=8, XTICKUNITS=['Hours'], $
+;      xrange     = [xmin, xmax], $
+         /NODATA
+
+   xyouts, 0.35, 0.8, $
+           'No valid direct normal data for this day.', $
+           /normal, charsize=1.5, color=red
+
+endif else begin
+
+   ymax = ymax*2.5
+
+   xmin = time_samples[0] - .1
+   xmax = time_samples[n_elements(l_dirnor_f2)-1] + .1
+
+   ;print, '  ... IRRAD plot: xmin, xmax=', xmin, xmax
+
+   xtitl   = ' '
+   if (which_mfr EQ 'nimfr') then begin
+      !x.charsize     = 3.0
+      xtitl           = 'Time of Day, GMT'
+   endif else begin
+      xtitl           = ' '
+      !x.charsize     = .001
+   endelse
+   
+                                ; DQ office prefers time axis labels in hh:mm form
+
+   dummy      = LABEL_DATE(DATE_FORMAT=['%H:%I'])
+
+   time_hours = TIMEGEN(START=time_samples[0], FINAL=(time_samples[0]+1.), UNITS='Hours')
+
+   plot, time_hours, l_dirnor_f5, $
+      ystyle = 1, $ 
+      xtitle = xtitl, $
+      ytitle = 'Direct Normal', $
+      background = white, $
+      color      = black, $
+      yrange     = [0.0, ymax], $
+      xstyle = 1, $
+      XTICKFORMAT='LABEL_DATE', XTICKS=8, XTICKUNITS=['Hours'], $
+   ;      xrange     = [xmin, xmax], $
+      /NODATA
+
+   ; Try to add the day/night background
+
+   day_night_background, l_lat, l_lon, /NOON
+ 
+   ; redraw the axes ...
+ 
+   ;print, 'IRRADIANCE PLOT ... yrange=', ymax
+   ;print, 'IRRADIANCE PLOT ...   ymax=', ymax
+ 
+   axis,color=black,charsize=0.01, xticks=1, $
+        XTICKFORMAT='LABEL_DATE', XTICKUNITS=['Hours'], xticklen=0.0
+   axis,color=black,yaxis=0,yrange=[0,ymax],ystyle=1,charsize=0.01,xtickunits='Time'
+   axis,color=black,yaxis=1,yrange=[0,ymax],ystyle=1,charsize=0.01,xtickunits='Time'
+
+   ; print, ' npts_f1, npts_f2, npts_f3, npts_f4, npts_f5: ', npts_f1, npts_f2, npts_f3, npts_f4, npts_f5 
+ 
+   if (npts_f1 GT 0) then begin
+     ;print, 'plotting filter1'
+     oplot, time_f1[0:npts_f1-1], curve_f1[0:npts_f1-1], color = blue, $
+        psym=1, symsize=1.5
+   endif
+   if (npts_f2 GT 0) then begin
+     ;print, 'plotting filter2'
+     oplot, time_f2[0:npts_f2-1], curve_f2[0:npts_f2-1], color = green, $
+        psym=1, symsize=1.5
+   endif
+   if (npts_f3 GT 0) then begin
+     ;print, 'plotting filter3'
+     oplot, time_f3[0:npts_f3-1], curve_f3[0:npts_f3-1], color = yellow, $
+        psym=1, symsize=1.5
+   endif
+   if (npts_f4 GT 0) then begin
+     ;print, 'plotting filter4'
+     oplot, time_f4[0:npts_f4-1], curve_f4[0:npts_f4-1], color = orange, $
+        psym=1, symsize=1.5
+   endif
+   if (npts_f5 GT 0) then begin
+     ;print, 'plotting filter5'
+     oplot, time_f5[0:npts_f5-1], curve_f5[0:npts_f5-1], color = red, $
+        psym=1, symsize=1.5
+   endif
+ 
+   if (which_mfr EQ 'mfrsr') then  begin
+     ; next plot (diffuse exponent)
+ 
+     !p.region = [0.1, 0.15, 1.0, 0.45]
+     !y.ticks = 2
+     !y.charsize = 3
+     !x.charsize = 3
+ 
+     curve           = l_diffuse_f1
+     num_pts         = n_elements(curve)
+     start_line      = 0
+     max_time_diff   = 0.0055
+ 
+     curve_f1               = l_dirnor_f2
+     curve_f2               = l_dirnor_f2
+     curve_f3               = l_dirnor_f2
+     curve_f4               = l_dirnor_f2
+     curve_f5               = l_dirnor_f2
+     curve_f1[0:npts-1]     = 0
+     curve_f2[0:npts-1]     = 0
+     curve_f3[0:npts-1]     = 0
+     curve_f4[0:npts-1]     = 0
+     curve_f5[0:npts-1]     = 0
+     time_f1                = time_samples
+     time_f1                = time_samples
+     time_f1                = time_samples
+     time_f1                = time_samples
+     time_f1                = time_samples
+     time_f1[0:npts-1]      = 0
+     time_f2[0:npts-1]      = 0
+     time_f3[0:npts-1]      = 0
+     time_f4[0:npts-1]      = 0
+     time_f5[0:npts-1]      = 0
+     npts_f1                = -1
+     npts_f2                = -1
+     npts_f3                = -1
+     npts_f4                = -1
+     npts_f5                = -1
+ 
+     sum_f1                 = 0
+     sum_f2                 = 0
+     sum_f3                 = 0
+     sum_f4                 = 0
+     sum_f5                 = 0
+ 
+     for n  = 1, npts-1 do begin
+ ;
+       w_value = l_qc_diffuse_f1[n]
+ 
+       if ( w_value EQ 0 AND (l_diffuse_f1[n] > (-9000.0)) ) then begin
+         npts_f1           = npts_f1 + 1
+         time_f1[npts_f1]  = time_samples[n]
+         curve_f1[npts_f1] = l_diffuse_f1[n]
+         sum_f1            = sum_f1 + l_diffuse_f1[n]
+       endif
+     endfor
+     for n  = 1, npts-1 do begin
+ ;
+       w_value = l_qc_diffuse_f2[n]
+   
+       if ( w_value EQ 0 AND (l_diffuse_f2[n] > (-9000.0)) ) then begin
+         npts_f2           = npts_f2 + 1
+         time_f2[npts_f2]  = time_samples[n]
+         curve_f2[npts_f2] = l_diffuse_f2[n]
+         sum_f2            = sum_f2 + l_diffuse_f2[n]
+       endif
+     endfor
+     for n  = 1, npts-1 do begin
+ 
+       w_value = l_qc_diffuse_f3[n]
+   
+       if ( w_value EQ 0 AND (l_diffuse_f3[n] > (-9000.0)) ) then begin
+         npts_f3           = npts_f3 + 1
+         time_f3[npts_f3]  = time_samples[n]
+         curve_f3[npts_f3] = l_diffuse_f3[n]
+         sum_f3            = sum_f3 + l_diffuse_f3[n]
+       endif
+     endfor
+     for n  = 1, npts-1 do begin
+ ;
+       w_value = l_qc_diffuse_f4[n]
+ 
+       if ( w_value EQ 0 AND (l_diffuse_f4[n] > (-9000.0)) ) then begin
+         npts_f4           = npts_f4 + 1
+         time_f4[npts_f4]  = time_samples[n]
+         curve_f4[npts_f4] = l_diffuse_f4[n]
+         sum_f4            = sum_f4 + l_diffuse_f4[n]
+       endif
+     endfor
+     for n  = 1, npts-1 do begin
+ ;
+       w_value = l_qc_diffuse_f5[n]
+ 
+       if ( w_value EQ 0 AND (l_diffuse_f5[n] > (-9000.0)) ) then begin
+         npts_f5           = npts_f5 + 1
+         time_f5[npts_f5]  = time_samples[n]
+         curve_f5[npts_f5] = l_diffuse_f5[n]
+         sum_f5            = sum_f5 + l_diffuse_f5[n]
+       endif
+     endfor
+
+     averages         = fltarr(5)
+     averages[0:4]    = 0
+
+     if (npts_f1 GT 0) then begin
+       averages[0] = sum_f1 / npts_f1
+   ;    print, 'sum_f1=', sum_f1, ' npts_f1=', npts_f1
+     endif
+     if (npts_f2 GT 0) then begin
+       averages[1] = sum_f2 / npts_f2
+     endif
+     if (npts_f3 GT 0) then begin
+       averages[2] = sum_f3 / npts_f3
+     endif
+     if (npts_f4 GT 0) then begin
+       averages[3] = sum_f4 / npts_f4
+     endif
+     if (npts_f5 GT 0) then begin
+       averages[4] = sum_f5 / npts_f5
+     endif
+
+     ;  print, 'averages[0]=', averages[0]
+
+     ymax = max(averages)
+     ymax = ymax*1.5
+
+     ;print, ' IRRAD DIFFUSE: ymax=', ymax
+
+     if (ymax GT 0.0) then begin
+ 
+       plot, time_hours, l_diffuse_f2, ystyle=1, $
+           ytitle = 'Diffuse Hemispheric', $
+           xtitle='Time of Day, GMT', $
+           yrange=[0.0, ymax], $
+           background=white, color=black, $
+           xstyle=1, $
+           XTICKFORMAT='LABEL_DATE', XTICKS=8, XTICKUNITS=['Hours'], $
+     ;        xrange=[xmin,xmax], $
+           /NODATA
+
+     ; Try to add the day/night background
+
+       day_night_background, l_lat, l_lon, /NOON
+ 
+     ; redraw the axes ...
+ 
+       ;print, 'IRRADIANCE PLOT2 ... yrange=', ymax
+ 
+       axis,color=black,charsize=0.01, xticks=1, $
+          XTICKFORMAT='LABEL_DATE', XTICKUNITS=['Hours'], xticklen=0.0
+       axis,color=black,yaxis=0,yrange=[0,ymax],ystyle=1,charsize=0.01,xtickunits='Time'
+       axis,color=black,yaxis=1,yrange=[0,ymax],ystyle=1,charsize=0.01,xtickunits='Time'
+
+       ;print, ' IRRAD DIFFUSE ... npts_f1, npts_f2, npts_f3, npts_f4, npts_f5: ', npts_f1, npts_f2, npts_f3, npts_f4, npts_f5 
+ 
+       if (npts_f1 GT 0) then begin
+         oplot, time_f1[0:npts_f1-1], curve_f1[0:npts_f1-1], color=blue, $ 
+            psym=1, symsize=1.5
+       endif
+       if (npts_f2 GT 0) then begin
+         oplot, time_f2[0:npts_f2-1], curve_f2[0:npts_f2-1], color=green, $ 
+            psym=1, symsize=1.5
+       endif
+       if (npts_f3 GT 0) then begin
+         oplot, time_f3[0:npts_f3-1], curve_f3[0:npts_f3-1], color=yellow, $ 
+            psym=1, symsize=1.5
+       endif
+       if (npts_f4 GT 0) then begin
+         oplot, time_f4[0:npts_f4-1], curve_f4[0:npts_f4-1], color=orange, $ 
+            psym=1, symsize=1.5
+       endif
+       if (npts_f5 GT 0) then begin
+         oplot, time_f5[0:npts_f5-1], curve_f5[0:npts_f5-1], color=red, $ 
+              psym=1, symsize=1.5
+      endif
+    endif
+   endif
+ endelse  ; end of big if block (ymax LE 0)
+ 
+; annotate all plots in the plot at the bottom
+
+!p.region = [0., 0., 1., .1]
+
+;xyouts, 0.3, 0.075, $
+;      'Optical depths calculated from ' + data_used_str, $
+;       /normal, charsize=1.0, color=black
+xyouts, 0.3, 0.055, $
+      'Created on: ' + create_date_str, $
+      /normal, charsize=1.0, color=black
+xyouts, 0.3, 0.035, $
+      'Quicklook Version: ' + l_idlVersion, $ 
+      /normal, charsize=1.0, color=black
+xyouts, 0.3, 0.015, $
+      'mfrod1barnmich Version: ' + create_ver_str, $
+      /normal, charsize=1.0, color=black
+
+if (which_mfr EQ 'nimfr') then begin
+  xyouts, 450, 730, title_str, color=black, /device, charsize=1.0, charthick=1.25
+endif
+if (which_mfr EQ 'mfrsr') then begin
+  xyouts, 450, 730, title_str, color=black, /device, charsize=1.0, charthick=1.25
+endif
+
+
+xyouts, .08,  0.115,'Nominal Wavelengths', color=black, /normal, charsize=0.9
+
+xyouts, .08, 0.095, '+' , $
+         color=blue, /normal, charsize=0.9
+xyouts, .10, 0.095, '= 415 nm (filter 1)', $
+          color=black, /normal, charsize=0.9
+ 
+xyouts, .08, 0.075, '+' , $
+         color=green, /normal, charsize=0.9
+xyouts, .10, 0.075, '= 500 nm (filter 2)', $
+          color=black, /normal, charsize=0.9
+ 
+xyouts, .08, 0.055, '+' , $
+         color=yellow, /normal, charsize=0.9
+xyouts, .10, 0.055, '= 615 nm (filter 3)', $
+          color=black, /normal, charsize=0.9
+ 
+xyouts, .08, 0.035, '+' , $
+         color=orange, /normal, charsize=0.9
+xyouts, .10, 0.035,  '= 673 nm (filter 4)', $
+          color=black, /normal, charsize=0.9
+ 
+xyouts, .08, 0.015, '+' , $
+         color=red, /normal, charsize=0.9
+xyouts, .10, 0.015, '= 870 nm (filter 5)', $
+          color=black, /normal, charsize=0.9
+ 
+
+end

src/do_od_plots.pro

+pro do_od_plots, which_mfr, $
+                 date, time_samples, create_ver_str, create_date_str, l_idlVersion, $
+                 Langley_used_str, data_used_str, $
+                 l_aod_f1,       l_aod_f2,       l_aod_f3,       l_aod_f4,       l_aod_f5, $
+                 l_qc_aod_f1,    l_qc_aod_f2,    l_qc_aod_f3,    l_qc_aod_f4,    l_qc_aod_f5, $
+                 l_dirnor_f1,    l_dirnor_f2,    l_dirnor_f3,    l_dirnor_f4,    l_dirnor_f5, $
+                 l_qc_dirnor_f1, l_qc_dirnor_f2, l_qc_dirnor_f3, l_qc_dirnor_f4, l_qc_dirnor_f5, $
+                 l_diffuse_f1, l_diffuse_f2, l_diffuse_f3, l_diffuse_f4, l_diffuse_f5, $
+                 l_qc_diffuse_f1, l_qc_diffuse_f2, l_qc_diffuse_f3, l_qc_diffuse_f4, l_qc_diffuse_f5, $
+;                 sunrise, sunset, $
+	         l_angstrom, l_qc_angstrom
+
+@color_syms.include
+usersym, d_circle, /fill
+;dloadct, 42
+
+
+Compile_Opt DEFINT32
+
+str_tmp         = string(format='(I0)',date)
+str_year        = strmid(str_tmp, 0, 4)
+str_month       = strmid(str_tmp, 4, 2)
+str_day         = strmid(str_tmp, 6, 2)
+
+str_date2       = str_year+'/'+str_month+'/'+str_day
+
+!P.Multi        = [0, 1, 4, 0, 1]  ; erase the page, 1 column of 5 plots, top to bottom
+
+; This plot region is only for getting a legible title, since we set the
+; x.charsize in the first real plot to a very small value (so as to not get
+; label on our x-axis ticks.
+
+!p.region       = [0.0, .9, 1., 1.]
+
+title_str       = strupcase(which_mfr) + ' Aerosol Optical Depths for '+str_date2
+
+
+!p.region       = [0.1, .5, 1, .95]
+!x.charsize     = .001
+!y.charsize     = 2.0
+!x.ticklen      = -0.02
+!y.ticklen      = -0.005
+!y.tickformat   = '(f6.1)'
+
+; First real plot (aerosol_optical_depth filter 2)
+
+npts = n_elements(l_aod_f2)
+w_qc = intarr(npts)
+w_qc = l_aod_f2
+
+bits    = intarr(10)
+qc_bits = intarr(10)
+qc_bits = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512]
+
+allbad  = 0L
+w_value = 0L
+
+;
+; Connor wants to plot points as yellow if qc is indeterminate
+; and to plot points as green if qc is good
+; not plotting obviously bad points
+;
+; for AOD, QC bit 6 to 10 are an indeterminate bit
+;              
+
+curve_yellow               = l_aod_f2
+curve_green                = l_aod_f2
+curve_yellow[0:npts-1]     = 0
+curve_green[0:npts-1]      = 0
+time_yellow                = time_samples
+time_green                 = time_samples
+time_yellow[0:npts-1]      = 0
+time_green[0:npts-1]       = 0
+npts_yellow                = -1
+npts_green                 = -1
+
+for n  = 1, npts-1 do begin
+;
+  w_value = l_qc_aod_f2[n]
+
+  bits[0:9] = 0
+  for i = 9, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if any of bits[6:9] are set, the value is indeterminate
+;
+  if ( (bits[0] > 0 OR bits[1] > 0 OR bits[2] > 0 OR bits[3] > 0 OR bits[4] > 0 OR bits[5] > 0) ) then continue
+  
+  if ( (bits[6] > 0 OR bits[7] > 0 OR bits[8] > 0 OR bits[9] > 0) AND (l_aod_f2[n] > (-9000.0)) ) then begin
+    npts_yellow               = npts_yellow + 1
+    time_yellow[npts_yellow]  = time_samples[n]
+    curve_yellow[npts_yellow] = l_aod_f2[n]
+  endif else begin
+    npts_green                = npts_green + 1
+    time_green[npts_green]    = time_samples[n]
+    curve_green[npts_green]   = l_aod_f2[n]
+  endelse
+endfor
+
+;max_y           = 2.0
+
+xmin = time_samples[0] - .1
+xmax = time_samples[n_elements(l_aod_f2)-1] + .1
+plot, time_samples, l_aod_f2, $
+      ystyle = 1, xstyle = 1, $
+      ytitle = 'Aerosol Optical Depth Filter 2', $
+      background = d_white, $
+      color      = d_black, $
+      yrange     = [0, 0.5], $
+      xrange     = [xmin, xmax], $
+      /NODATA
+
+if (npts_yellow GE 0) then begin
+  oplot, time_yellow[0:npts_yellow-1], curve_yellow[0:npts_yellow-1], color = d_orange, $
+     psym=3, symsize=3.
+endif
+
+;print, npts_green 
+;print, time_green[0:npts_green-1]
+;print, curve_green[0:npts_green-1]
+
+if (npts_green GE 0) then begin
+  oplot, time_green[0:npts_green-1], curve_green[0:npts_green-1], color = d_green, $
+     psym=3, symsize=3.
+endif
+
+;print, npts_yellow  
+;print, time_yellow[0:npts_yellow-1]
+;print, curve_yellow[0:npts_yellow-1]
+
+; next plot (angstrom exponent)
+
+!p.region = [0.1, 0.15, 1.0, 0.45]
+!y.ticks = 2
+!y.charsize = 2
+!x.charsize = 2
+
+curve           = l_angstrom
+num_pts         = n_elements(curve)
+start_line      = 0
+max_time_diff   = 0.0055
+
+curve_yellow[0:npts-1]     = 0
+curve_green[0:npts-1]      = 0
+time_yellow                = time_samples
+time_green                 = time_samples
+time_yellow[0:npts-1]      = 0
+time_green[0:npts-1]       = 0
+npts_yellow                = -1
+npts_green                 = -1
+
+for n  = 1, npts-1 do begin
+;
+  w_value = l_qc_angstrom[n]
+
+  bits[0:9] = 0
+  for i = 5, 0, -1 do begin
+    if (w_value GE qc_bits[i]) then begin
+       bits[i] = 1
+       w_value = w_value - qc_bits[i]
+    endif
+  endfor
+
+; if bits[4] is set, the value of angstrom is indeterminate
+;
+  if ( (bits[0] > 0 OR bits[1] > 0 OR bits[2] > 0 OR bits[3] ) ) then continue
+  
+  if ( (bits[4] > 0 ) AND (l_angstrom[n] > (-9000.0)) ) then begin
+    npts_yellow               = npts_yellow + 1
+    time_yellow[npts_yellow]  = time_samples[n]
+    curve_yellow[npts_yellow] = l_angstrom[n]
+  endif else begin
+    npts_green                = npts_green + 1
+    time_green[npts_green]    = time_samples[n]
+    curve_green[npts_green]   = l_angstrom[n]
+  endelse
+endfor
+
+plot, time_samples, curve, ystyle=1, xstyle=1, $
+      ytitle = 'Angstrom Exponent', $
+      xtitle='Julian Day', $
+      yrange=[-0.5, 1.5], $
+      background=d_white, color=d_black, $
+      xrange=[xmin,xmax], /NODATA
+
+if (npts_yellow GE 0) then begin
+  oplot, time_yellow, curve_yellow, color=d_orange, $ 
+       psym=3, symsize=3.
+endif
+if (npts_green GE 0) then begin
+  oplot, time_green, curve_green, color=d_green, $ 
+       psym=3, symsize=3.
+endif
+
+
+; annotate all plots in the plot at the bottom
+
+!p.region = [0., 0., 1., .1]
+
+;xyouts, 0.7, 0.08, 'Langley information taken from '+Langley_used_str, $
+;      /normal, charsize=0.8, color=d_black
+;xyouts, 0.7, 0.06, $
+;      'Optical depths calculated from ' + data_used_str, $
+;       /normal, charsize=0.8, color=d_black
+xyouts, 0.7, 0.04, $
+      'Created on: ' + create_date_str, $
+      /normal, charsize=0.8, color=d_black
+xyouts, 0.7, 0.02, $
+      'Quicklook Version: ' + l_idlVersion, $ 
+      /normal, charsize=0.8, color=d_black
+xyouts, 0.7, 0.00, $
+      'mfrod1barnmich Version: ' + create_ver_str, $
+      /normal, charsize=0.8, color=d_black
+
+xyouts, 450, 730, title_str, color=d_black, /device, charsize=1., charthick=1.25
+
+
+end

src/forgan.h

+#ifndef FORGAN_H
+#define FORGAN_H
+
+typedef struct {
+	double	val;
+	double	lowess;
+} pointType;
+
+typedef struct
+{
+  int   	windowSize;	/* how big a window ? */
+  int	        slideSize;	/* how much to slide a window ? */
+  int	        winCount;	/* current window number */
+  int	        firstIo;	/* index of first Io for current window */
+  int	        lastIo;		/* index of last Io for current window */
+  int	        numSamples;     /* number of points */
+  int	        numChannels;    /* number of columns in Io */
+  double	Ios[2000][7];	/* array of Io's by channel  .... MFRs spit out 4320 samples per day, we're doubling that*/
+  double        IosTimes[2000];
+
+  int 	        ratio1;	/* column index used for numerator of forgan ratio */
+  int 	        ratio2;	/* column index used for denomenator of forgan ratio */
+  double 	sigmaFactor;
+  double        minratio;
+  double        maxratio;   /* just kept around for debugging output */
+} forganType;
+
+typedef enum { good, bad } point_type;
+
+typedef struct {
+	point_type status;
+	double	   val;
+} ratio_type;
+
+#ifndef TRUE
+#define TRUE 1
+#define FALSE 0
+#endif
+
+#endif

src/gauss_filter.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "forgan.h"
+
+void convert_time_to_year_frac(double  in_time,
+			       double *out_time)
+{
+  
+  double ask_time;
+  int    num_years;
+  int    yr_tmp;
+  int    kk;
+  double time_ask;
+
+  ask_time  = in_time;
+  ask_time  = ask_time - 25568.0 + 1; /* ask_time is now days since 1970 */
+
+  num_years = ask_time / 365.0;  /* an approximate number of years since 1970 */
+  yr_tmp    = 1970;
+
+  for (kk = 0; kk < num_years; kk++)
+    {
+      if (yr_tmp == 1972 || yr_tmp == 1976 || yr_tmp == 1980 ||
+	  yr_tmp == 1984 || yr_tmp == 1988 || yr_tmp == 1992 ||
+	  yr_tmp == 1996 || yr_tmp == 2000 || yr_tmp == 2004 ||
+	  yr_tmp == 2008 || yr_tmp == 2012 || yr_tmp == 2016 ||
+	  yr_tmp == 2020 || yr_tmp == 2024 || yr_tmp == 2028 ||
+	  yr_tmp == 2032)
+	{
+	  ask_time = ask_time - 366; /* Subtract the number of days in a leap year */
+	}
+      else
+	{
+	  ask_time = ask_time - 365; /* non leap year */
+	}
+      yr_tmp++;
+    }
+	      
+  time_ask = yr_tmp;
+  if (time_ask == 1972 || time_ask == 1976 || time_ask == 1980 ||
+      time_ask == 1984 || time_ask == 1988 || time_ask == 1992 ||
+      time_ask == 1996 || time_ask == 2000 || time_ask == 2004 ||
+      time_ask == 2008 || time_ask == 2012 || time_ask == 2016 ||
+      time_ask == 2020 || time_ask == 2024 || time_ask == 2028 ||
+      time_ask == 2032)
+    {
+      time_ask = time_ask + (ask_time / 366.0); /* Subtract the number of days in a leap year */
+    }
+  else
+    {
+      time_ask = time_ask + (ask_time / 365.0); /* non leap year */
+    }
+
+  /* printf ("time_in: %lf ... time_out: %lf\n", in_time, time_ask); */
+
+  *out_time = time_ask;
+} /* End of convert_time_to_year_frac() function */
+
+void convert_year_frac_to_mfr_time(double  yr_frac,
+				   double *mfr_time)
+{
+  
+  double ask_time;
+  double frac;
+  int    num_years;
+  int    yr_tmp;
+  int    tmp;
+  int    kk;
+ 
+  yr_tmp    = (int) yr_frac;    /* We enter this function with a value like 2011.08219178, so "yr_tmp" is 2011 */
+  frac      = yr_frac - yr_tmp; /* frac is 0.08219178 */
+  tmp       = yr_tmp;           /* tmp is 2011 */
+
+  num_years = yr_tmp - 1970;
+
+  ask_time  = 0.0;
+
+  /* We need to compute the number of days (total) between 1970 and our current year (2011 in this case) */
+
+  /* This kk loop is summing up the number of days for the years between 1970 and 2011 */
+  for (kk = 1970; kk < (int) tmp; kk++)
+    {
+      if (kk == 1972 || kk == 1976 || kk == 1980 ||
+	  kk == 1984 || kk == 1988 || kk == 1992 ||
+	  kk == 1996 || kk == 2000 || kk == 2004 ||
+	  kk == 2008 || kk == 2012 || kk == 2016 ||
+	  kk == 2020 || kk == 2024 || kk == 2028 ||
+	  kk == 2032)
+	{
+	  ask_time = ask_time + 366; /* Add the number of days in a leap year */
+	}
+      else
+	{
+	  ask_time = ask_time + 365; /* Add the number of days in a non leap year */
+	}
+      /* printf ("kk=%d ask_time=%lf\n", kk, ask_time); */
+    }
+
+  /* At this point, ask_time is the number of days since 1970.
+   * Now, add the portion accounted for by "frac"
+   */
+
+  if (yr_tmp == 1972 || yr_tmp == 1976 || yr_tmp == 1980 ||
+      yr_tmp == 1984 || yr_tmp == 1988 || yr_tmp == 1992 ||
+      yr_tmp == 1996 || yr_tmp == 2000 || yr_tmp == 2004 ||
+      yr_tmp == 2008 || yr_tmp == 2012 || yr_tmp == 2016 ||
+      yr_tmp == 2020 || yr_tmp == 2024 || yr_tmp == 2028 ||
+      yr_tmp == 2032)
+    {
+      ask_time = ask_time + (frac * 366.0); 
+    }
+  else
+    {
+      ask_time = ask_time + (frac * 365.0); 
+    }
+
+  /* Now, add the number of days between 1900 and 1970 */
+
+  ask_time = ask_time + 25568.0;
+
+  *mfr_time = ask_time;
+}
+	      
+void gauss_filter(int     vap_FirstChannel,
+		  int     vap_LastChannel,
+		  long    StartDate,
+		  long    EndDate,
+		  double  final_times[3000],
+		  double  final_ios[7][3000])
+{
+  /* Jim Barnard's "smoother_arm.f" rewritten in C */
+
+  double  new_time[10000];
+  double  var[7][10000];
+  double  sum;
+  double  arg;
+  double  time_ask[3000];
+  double  my_time;
+  
+  extern  forganType ft;               /* Defined in forgan.h */
+  extern  double     ForgIos[3000][7];  /* Also defined in forgan.c */
+  extern  double     ForgIosTimes[3000]; /* Also defined in forgan.c */
+
+  int     day_count;
+  int     day;
+  int     i, j;
+  int     ip, ip1;
+  int     total_days;
+  int     l;
+
+  double  start;
+  
+  double  syr;
+  double  afactor;
+  double  tmp_time;
+
+  double  x1, x2, x3, y1, y2, y3;
+  double  a, b, c, delta;
+
+  for (day=0; day < 3000; day++)
+    {
+      final_times[day] = 0.0;
+      for (j=0; j <= 6; j++)
+	{
+	  final_ios[j][day] = 0.0;
+	}
+    }
+  for (day=0; day < 10000; day++)
+    {
+      for (j=0; j <= 6; j++)
+	{
+	  var[j][day] = 0.0;
+	}
+    }
+
+  /* Start here */
+
+  syr     = 0.10;  /* Based on "experience", per Jim Barnard's comments in original Fortran */
+  afactor = syr / 2.0 / sqrt(log(2.0));
+
+  convert_time_to_year_frac(ForgIosTimes[0], 
+			    &start);
+
+  /* Convert the times in the "ft" structure to year.fraction times */
+
+  for (i=0; i < (int) (ft.winCount-1); i++)
+    {
+      /* Use convert_time_to_year_frac instead */
+
+      convert_time_to_year_frac(ForgIosTimes[i],
+				&time_ask[i]);
+    }
+  
+  for(j=vap_FirstChannel; j <= vap_LastChannel; j++) 
+    {
+      /* Start the smoothing loop */
+
+      day_count = 0;
+
+      for (day = StartDate; day <= EndDate; day++) /* These dates are in MFRSR time, like 35491 to 35641 */
+	{
+
+	  var[j][day_count]   = 0.0;
+	  sum                 = 0.0;
+
+	  convert_time_to_year_frac ( (double) day,
+				      &tmp_time);
+	  new_time[day_count] = tmp_time; /* new_time is in a fractional year AD (40618 in MFR time becomes 2011.2054 */
+
+	  for (i=0; i < (int) (ft.winCount-1); i++) 
+	    {
+	      /* Jim B. says to multiply by 365 to use times in MFR-time form */
+
+	      /* Compute in time in a manner that will allow this gaussian filter
+	       * logic to play nicely with Jim's Gaussian filter logic. Jim's
+	       * logic seems to depend on having the date expressed as a floating
+	       * point number, which is the year, and the fraction of the year.
+	       *
+	       * That is, 1997.125 means a date that is .125 (1/8th) of the way
+	       * through 1997. This is 45 days approximately past Jan 1, 1997.
+	       * 
+	       * The following logic attempts to convert the ft.IosTimes values
+	       * from MFRSR times, which are expressed as floating point days since
+	       * Jan. 1, 1900, in terms of year.fract.
+	       */
+
+	      arg               = fabs(new_time[day_count] - time_ask[i]) / (afactor); 
+
+	      var[j][day_count] = var[j][day_count] + 
+                        	  (exp(-1.0 * (arg * arg)) * ForgIos[i][j]);
+	      sum               = sum + (exp(-1.0 * (arg * arg)));
+	    } /* time loop */
+
+	  if (sum > 0.0)
+	    {
+	      var[j][day_count] = var[j][day_count] / sum;
+	    }
+	  else
+	    {
+	      var[j][day_count] = -9999.0;
+	    }
+
+	  day_count++;
+
+	} /* End of "day" loop */
+
+    } /* End of "j" (channel loop) */
+
+  /* From Jim: new end effects logic */
+
+  total_days = day_count-1;
+  ip         = total_days * 0.06;
+  l          = 1.0 * ip;
+
+  for (i=0; i < ip; i++)
+    {
+      for (j=0; j < 5; j++)
+	{
+	  x1=0.0;
+	  x2=new_time[ip+2*l]-new_time[ip];
+	  x3=new_time[ip+3*l]-new_time[ip];
+	  y1=var[j][ip];
+	  y2=var[j][ip+2*l];
+	  y3=var[j][ip+3*l];
+
+	  a = 0.0;
+	  b = 0.0;
+	  c = 0.0;
+	  /* c coefficients from Mathematica - saves me work: */
+	  a=(x3*(-y1 + y2) + x2*(y1 - y3) + x1*(-y2 + y3))/
+	    ((x1 - x2)*(x1 - x3)*(x2 - x3));
+	  b=(x3*x3*(y1 - y2) + x1*x1*(y2 - y3) + x2*x2*(-y1 + y3))/
+	    ((x1 - x2)*(x1 - x3)*(x2 - x3));
+	  c=(x3*(x2*(x2 - x3)*y1 + x1*(-x1 + x3)*y2) + x1*(x1 - x2)*x2*y3)/
+	    ((x1 - x2)*(x1 - x3)*(x2 - x3));
+	  delta=new_time[i]-new_time[ip];
+	  var[j][i]=delta*delta*a+delta*b+c;
+	}
+    }
+
+  /* c other end */
+
+  ip1 = day_count - ip;
+
+  for (i = ip1; i < day_count; i++) /*  other end */
+    {
+      for (j=0; j < 5; j++)
+	{
+	  x1=new_time[ip1-3*l] - new_time[ip1];
+	  x2=new_time[ip1-2*l] - new_time[ip1];
+	  x3=0.0;	
+	  y1=var[j][ip1-3*l];
+	  y2=var[j][ip1-2*l];
+	  y3=var[j][ip1];
+	  a = 0.0;
+	  b = 0.0;
+	  c = 0.0;
+	  /* c coefficients from Mathematica - saves me work: */
+	  a=(x3*(-y1 + y2) + x2*(y1 - y3) + x1*(-y2 + y3))/
+	    ((x1 - x2)*(x1 - x3)*(x2 - x3));
+	  b=(x3*x3*(y1 - y2) + x1*x1*(y2 - y3) + x2*x2*(-y1 + y3))/
+	    ((x1 - x2)*(x1 - x3)*(x2 - x3));
+	  c=(x3*(x2*(x2 - x3)*y1 + x1*(-x1 + x3)*y2) + x1*(x1 - x2)*x2*y3)/
+	    ((x1 - x2)*(x1 - x3)*(x2 - x3));
+	  delta=new_time[i] - new_time[ip1];
+	  var[j][i]=delta*delta*a+delta*b+c;
+	}
+    }
+  
+  /* Now, write out the smoothed results */
+
+  for (day=0; day < day_count; day++)
+    {
+      convert_year_frac_to_mfr_time (new_time[day], 
+				     &my_time);
+
+      final_times[day] = my_time;
+
+      for (j=vap_FirstChannel; j <= vap_LastChannel; j++)
+	{
+	  final_ios[j][day] = var[j][day];
+	}
+
+    }
+
+} /* End of gauss_filter() */
+ 
+
+

src/io_C1.c

+#include <stdio.h>
+#include <stdlib.h>
+/* define function */
+
+/* returns Io for the MFRSR C1 unit */
+/* currently valid for 1996/09/17 (35324.81) through 1998/08/03 (36009.79) */
+/* ichan is the channel number 1 = 415 nm,..., 5 = 873 nm */
+
+void get_io_C1(double mfrsr_time, float *Ios)
+{
+  double      jdate[10000],dummy;
+  float       io[10000][5];
+  FILE       *fptr;
+  int         IEOF,i,j,k;
+  static int  index=0;
+  static int  numrecords;
+  static int  counter;
+
+  char        filename[256];
+
+  /* end of declarations */
+
+  if(index==0)
+    {    
+
+      sprintf (filename, "%s/%s", getenv("VAP_CONF_HOME"), "mfrod1barnmich_special_Io.info");
+
+      fptr=fopen(filename,"r"); /* open file; assume it exists */
+      /* read data */
+
+      numrecords = 0;
+
+      for(i=0;i<=10000;i++)
+	{	
+	  IEOF=fscanf(fptr,"%lf",&jdate[i]);
+	  for(k=0;k<=4;k++)
+	    {
+	      fscanf(fptr,"%f",&io[i][k]);
+	    }
+	  IEOF=fscanf(fptr,"%lf",&dummy); 
+	  if(IEOF==-1)
+	    {
+	      break;
+	    }
+	  numrecords++;
+	}
+
+      numrecords--;
+      counter = 0;
+      index=1; /* don't read the dang file again */
+    }
+
+  /* OK grab the right io */
+  /* check that the MFRSR time lies within the span of time in which we have Ios */
+
+  if(mfrsr_time < jdate[0] || mfrsr_time > jdate[numrecords])
+    {
+
+      printf ("ERROR: io_C1 function will not handle time=%lf\n",
+	      mfrsr_time);
+
+      return;
+    }
+  /* find closest time greater than the mfrsr_time - inefficient code but computers are powerful */
+  for (j=counter; j<=numrecords; j++)
+    {
+      if(jdate[j] > mfrsr_time)
+	{
+
+	  Ios[0] = io[j][0];
+	  Ios[1] = io[j][1];
+	  Ios[2] = io[j][2];
+	  Ios[3] = io[j][3];
+	  Ios[4] = io[j][4];
+
+	  counter = j; /* save for next time in here */
+	  return;
+	}
+    }
+}
+
+
+

src/lfit.c

+/*
+ * $Id: lfit.c,v 1.1 2002-03-06 13:36:39 koontz Exp $
+ *
+ * lfit() does a weighted linear fit, returning all sorts of error
+ * numbers.  You can unweight by using NULL for the w array on input.
+ *
+ *
+ * $Log: not supported by cvs2svn $
+ * Revision 1.2  1997/07/22 16:27:45  d3a230
+ * Added some sanity checking when we have two or fewer non-zero-weighted
+ * points.  Pathological, of course, but if we don't do it we get math errors
+ * all the time.
+ *
+ * Revision 1.1  1997/07/18 15:58:43  d3a230
+ * Initial revision
+ *
+ */
+
+static char *rcsid="$Id: lfit.c,v 1.1 2002-03-06 13:36:39 koontz Exp $";
+static char *rcsstate="$State: Exp $";
+
+#include <math.h>
+/*
+ * Do a linear fit y = a+bx, weighted by w
+ *
+ * Right now the siga and sigb assume no measurement uncertainty - i.e.
+ * we estimate via sampling theory - sigy =~ sqrt(npr*syx2)
+ */
+
+void lfit(float *x, float *y, float *w, int n, float *a, float *b,
+	  float *siga, float *sigb, float *sigfit, int *count)
+{
+  int i;
+
+  double sumx=0, sumy=0, sumxy=0, sumx2=0, sumy2=0, sumw=0;
+  double denom, slope, intercept, npr, m, syx2;
+
+  m=0;
+
+  /* This lets us set w=NULL on input to do an unweighted fit */
+#define W (w ? w[i] : 1)
+
+  for (i=0; i<n; i++) {
+
+    sumx += W*x[i];
+    sumx2 += W*x[i]*x[i];
+    
+    sumy += W*y[i];
+    sumy2 += W*y[i]*y[i];
+
+    sumxy += W*x[i]*y[i];
+
+    sumw += W;
+
+    /* This is a little goofy - m holds the number of non-zero        */
+    /* weighted points, which we need to get our n/(n-2) factor right */ 
+    if (w[i] > 0.0) m += 1.0;
+
+  }
+#undef W
+
+  /* simple formulas */
+  denom = sumw*sumx2 - sumx*sumx;
+  slope = (sumw*sumxy - sumx*sumy)/denom;
+  intercept = (sumy*sumx2 - sumx*sumxy)/denom;
+
+  /* I'm not sure this is right, especially with the weights */
+  /* syx2 is standard error of estimate */
+  syx2 = (sumy2 - intercept*sumy - slope*sumxy)/sumw;
+
+  /* sampling error */
+  npr = m > 2 ? m/(m-2) : 1;
+
+  *a = intercept;
+  *b = slope;
+  if (siga) *siga = m>2 ? sqrt(syx2*npr*sumx2/denom) : 1e8;
+  if (sigb) *sigb = m>2 ? sqrt(syx2*npr*sumw/denom) : 1e8;
+  if (sigfit) *sigfit = m>2 ? sqrt(syx2) : 1e8;
+  if (count) *count = m;
+
+  return;
+}
+
+char *get_fit_id()
+{ 
+  return(rcsid);
+}
+
+char *get_fit_state()
+{ 
+  return(rcsstate);
+}

src/lowess.c

+#include <stdlib.h>
+#include <math.h>
+
+/*
+ * history
+ *
+ * 'round about 1/20/97
+ * Rob Seals hand translates Cleveland's RATFOR code to C.
+ *
+ * 'round about 1/30/97
+ * Jim Schlemmer takes off some rough edges
+ *
+ */
+
+#define	max(x,y)	((x) > (y) ? (x) : (y))
+#define	min(x,y)	((x) < (y) ? (x) : (y))
+
+
+/*****************************************************************************
+ *
+ * cube() function
+ *
+ *****************************************************************************
+ */
+double cube(double x)
+{
+    return x*x*x;
+} /* End of cube() function */
+
+
+/*****************************************************************************
+ *
+ * dcmp() function
+ *
+ *****************************************************************************
+ */
+
+int dcmp(const void *v1, const void *v2)
+{
+  double	*d1 = (double *) v1, 
+                *d2 = (double *) v2;
+
+  if (*d1 < *d2)
+    {
+      return -1;
+    }
+  else if (*d1 > *d2)
+    {
+      return 1;
+    }
+  else
+    {
+      return 0;
+    }
+} /* End of dcmp() function */
+
+
+/*****************************************************************************
+ *
+ * lowess() function
+ *
+ *    uses lowest() function
+ *    uses dcmp() via qsort()
+ *
+ *****************************************************************************
+ */
+
+int lowess(double *x, 
+	   double *y, 
+	   int     n, 
+	   double  f, 
+	   int     nsteps, 
+	   double  delta,
+	   double *ys)
+{
+  extern double sqr();
+  /*extern void  bw_abort(); */
+
+  double	*rw, *res; 
+  double	d1, d2, alpha, denom,
+                cut, c1, c9, r, cmad;
+
+  int		i, j, k, ns, iter, nleft,
+                nright, last, m1, m2;
+
+  int 	        dcmp(const void *, const void *);
+  int		lowest(double *x, double *y, int n, double xs, double *ys,
+		       int nleft, int nright, double *w, int userw,
+		       double *rw);
+
+  if (n < 2) 
+    {
+      ys[0] = y[0];
+      return 0;
+    }
+
+  /* "x" is the times, and "y" is the Ios values.
+   * 
+   * Modifying lowess() so that if a negative Ios is found
+   * it will not be used.
+   *
+   * This function uses the input times and Ios (x and y)
+   * and iterates (nsteps times) to compute a new
+   * value for the y's (which come back in the ys array)
+   *
+   */
+
+  /* NOTE: the forgan() allocates an "rw" and a "res"
+   *       arrays, too. 
+   */
+
+  rw  = calloc(n, sizeof(double));
+  if (rw == NULL)
+    {
+      /* bw_abort(-1, "Memory allocation error in lowess() "); */
+    }
+  res = calloc(n, sizeof(double));
+  if (res == NULL)
+    {
+      /*bw_abort(-1, "Memory allocation error in lowess() "); */
+    }
+
+  /* "at least 2, at most n points" */
+  ns = max(min((int) (f*n), n), 2);
+
+  for (iter=1; iter<=(nsteps+1); iter++) 
+    {
+      nleft  = 0;
+      nright = ns - 1;
+      last   = -1;
+      i      = 0;
+      do 
+	{
+	  
+	  while (nright < (n-1)) 
+	    {
+	      /* correction - used to be '(nright < n)' */
+	      d1 = x[i] - x[nleft];
+	      d2 = x[nright+1] - x[i];
+
+	      if (d1 <= d2)
+		{
+		  break;
+		}
+	      nleft++;
+	      nright++;
+	    }
+
+	  if (!lowest(x, 
+		      y, 
+		      n, 
+		      x[i], 
+		      ys + i, 
+		      nleft, 
+		      nright, 
+		      res, 
+		      iter>1, 
+		      rw))
+	    {
+	      ys[i] = y[i];
+	    }
+	    
+	  if (last < i-1) 
+	    {
+	      denom = x[i] - x[last];	/* which ought to be > -1 by now */
+	      for (j=last+i; j<i; j++) 
+		{
+		  alpha = (x[j] - x[last])/denom;
+		  ys[j] = alpha*ys[i] + (1.0 - alpha)*ys[last];
+		}
+	    }
+
+	  last = i;
+	  cut = x[last] + delta;
+	  for (k=last+1; k<n; k++) 
+	    {  /* correction - used to be 'i<=n'.  also changed 'i' to 'k' */
+	      if (x[k] > cut)
+		{
+		  break;
+		}
+	      if (x[k] == x[last]) 
+		{
+		  ys[k] = ys[last];
+		  last = k;
+		}
+	    }
+
+	  /* js - changed this comp from 1.0 to .01 because some calibrated Ios are
+	     less than 1.0 and this would zero them out.  Not sure what the function of
+	     this is anyway.
+	  */
+
+	  if(ys[i] < .01)
+	    {
+	      ys[i] = -9999.0;
+	    }
+
+	  i = max(last+1, i-1);
+	} while(last < (n-1));   /* correction - used to be (last >= n) */
+
+      for (k=0; k<n; k++)	
+	{
+	  res[k] = y[k] - ys[k];
+	}
+
+      if (iter > nsteps)
+	{
+	  break;
+	}
+
+      for (k=0; k<n; k++)
+	{
+	  rw[k] = abs(res[k]);
+	}
+
+      qsort(rw, n, sizeof(double), dcmp);
+
+      m1 = n/2;
+      m2 = n - m1 - 1;   /* correction - used to be 'm2 = n - m1' */
+      cmad = 3.0*(rw[m1] + rw[m2]);
+      c9 = 0.999*cmad;
+      c1 = 0.001*cmad;
+
+      for (k=0; k<n; k++) 
+	{
+	  r = abs(res[k]);
+	  if (r <= c1)
+	    {
+	      rw[k] = 1.0;
+	    }
+	  else if (r > c9)
+	    {
+	      rw[k] = 0.0;
+	    }
+	  else
+	    {
+	      rw[k] = sqr(1.0 - sqr(r/cmad));
+	    }
+	}
+    }
+
+  free(rw);
+  free(res);
+
+  return 1;
+} /* End of lowess() function */
+
+
+/*****************************************************************************
+ *
+ * lowest() function
+ *
+ *    uses cube() function
+ *    uses sqr() via qsort()
+ *
+ *****************************************************************************
+ */
+
+int lowest(double *x, 
+	   double *y, 
+	   int     n, 
+	   double  xs, 
+	   double *ys,
+	   int     nleft, 
+	   int     nright, 
+	   double *w, 
+	   int     userw, 
+	   double *rw)
+{
+
+  extern double sqr();
+
+  double	h, range, h1, h9, a, b, c, r;
+  int		j, nrt;
+
+  range = x[n-1] - x[0];
+  h     = max(xs - x[nleft], x[nright] - xs);
+  h9    = 0.999*h;
+  h1    = 0.001*h;
+
+  a     = 0.0;
+  for (j=nleft; j<n; j++) 
+    {
+      w[j] = 0.0;
+      r = abs(x[j] - xs);
+      if (r <= h9) 
+	{
+	  if (r > h1)
+	    {
+	      w[j] = cube(1.0 - cube(r/h));
+	    }
+	  else
+	    {
+	      w[j] = 1.0;
+	    }
+	  if (userw)
+	    {
+	      w[j] = rw[j]*w[j];
+	    }
+	  a += w[j];
+	} 
+      else if (x[j] > xs)
+	{
+	  break;
+	}
+    }
+
+    nrt = j-1;
+    if (a <= 0.0)
+      {
+	return 0;
+      }
+
+    for (j=nleft; j<=nrt; j++)
+      {
+	w[j] /= a;
+      }
+
+    if (h > 0.0) 
+      {
+	for (a=0.0, j=nleft; j<=nrt; j++)
+	  {
+	    a += w[j]*x[j];
+	  }
+	b = xs - a;
+	for (c=0.0, j=nleft; j<=nrt; j++)
+	  {
+	    c += w[j]*sqr(x[j] - a);
+	  }
+	if (sqrt(c) > 0.001*range) 
+	  {
+	    b /=c;
+	    for (j=nleft; j<=nrt; j++)
+	      {
+		w[j] *= (1.0 + b*(x[j] - a));
+	      }
+	  }
+      }
+
+    for (*ys=0.0, j=nleft; j<=nrt; j++)
+      {
+	if (y[j] > 0.0)
+	  {
+	    *ys += w[j]*y[j];
+	  }
+      }
+
+    return 1;
+} /* End of lowest() function */
+
+
+
+/*****************************************************************************
+ *
+ * sqr() function
+ *
+ *****************************************************************************
+ */
+
+double sqr(double x)
+{
+    return x*x;
+}
+
+

src/mfrod1barnmich_quicklooks.h

+/* Miscellaneous flags, etc. */
+
+int quicklook_flag = TRUE;
+int dev_flag       = FALSE;
+int save_flag      = FALSE;
+int use_x_device;

src/mfrod1barnmich_retriever.h

+#ifndef _MFROD1BARNMICH_RETRIEVER_H
+#define _MFROD1BARNMICH_RETRIEVER_H
+
+
+#include "bw_adi.h"
+
+#ifndef NO_RETRIEVER_H
+const char *Platform_List[]=
+  {
+    "mfrsr.a0",  /* 0 */
+    "mfrsr.b1",  /* 1 */
+    "nimfr.a0",  /* 2 */
+    "nimfr.b1",  /* 3 */
+    "mfrsrlangley.c1", /* 4 */
+    "omi.a1",    /* 5 */
+    "toms.a1",   /* 6 */
+    "met.b1",    /* 7 */
+    NULL
+  };
+
+const char *Field_List[][BW_MAX_FIELDS]=
+  {
+    {
+      //mfrsr.a0
+      "hemisp_broadband_raw",
+      "hemisp_narrowband_filter1_raw",
+      "hemisp_narrowband_filter2_raw",
+      "hemisp_narrowband_filter3_raw",
+      "hemisp_narrowband_filter4_raw",
+      "hemisp_narrowband_filter5_raw",
+      "hemisp_narrowband_filter6_raw",
+      "diffuse_hemisp_broadband_raw",
+      "diffuse_hemisp_narrowband_filter1_raw",
+      "diffuse_hemisp_narrowband_filter2_raw",
+      "diffuse_hemisp_narrowband_filter3_raw",
+      "diffuse_hemisp_narrowband_filter4_raw",
+      "diffuse_hemisp_narrowband_filter5_raw",
+      "diffuse_hemisp_narrowband_filter6_raw",
+      "lat_a0",
+      "lon_a0",
+      "alt_a0",
+      NULL
+    },
+    {
+      //mfrsr.b1
+      "hemisp_broadband",                      /* fn=0 */
+      "qc_hemisp_broadband",                   /* fn=1 */
+      "hemisp_narrowband_filter1",             /* fn=2 */
+      "qc_hemisp_narrowband_filter1",          /* fn=3 */
+      "hemisp_narrowband_filter2",             /* fn=4 */
+      "qc_hemisp_narrowband_filter2",          /* fn=5 */
+      "hemisp_narrowband_filter3",             /* fn=6 */
+      "qc_hemisp_narrowband_filter3",          /* fn=7 */
+      "hemisp_narrowband_filter4",             /* fn=8 */
+      "qc_hemisp_narrowband_filter4",          /* fn=9 */
+      "hemisp_narrowband_filter5",             /* fn=10 */
+      "qc_hemisp_narrowband_filter5",          /* fn=11 */
+      "hemisp_narrowband_filter6",             /* fn=12 */
+      "qc_hemisp_narrowband_filter6",          /* fn=13 */
+      "diffuse_hemisp_broadband",              /* fn=14 */
+      "qc_diffuse_hemisp_broadband",           /* fn=15 */
+      "diffuse_hemisp_narrowband_filter1",     /* fn=16 */
+      "qc_diffuse_hemisp_narrowband_filter1",  /* fn=17 */
+      "diffuse_hemisp_narrowband_filter2",     /* fn=18 */
+      "qc_diffuse_hemisp_narrowband_filter2",  /* fn=19 */
+      "diffuse_hemisp_narrowband_filter3",     /* fn=20 */
+      "qc_diffuse_hemisp_narrowband_filter3",  /* fn=21 */
+      "diffuse_hemisp_narrowband_filter4",     /* fn=22 */
+      "qc_diffuse_hemisp_narrowband_filter4",  /* fn=23 */
+      "diffuse_hemisp_narrowband_filter5",     /* fn=24 */
+      "qc_diffuse_hemisp_narrowband_filter5",  /* fn=25 */
+      "diffuse_hemisp_narrowband_filter6",     /* fn=26 */
+      "qc_diffuse_hemisp_narrowband_filter6",  /* fn=27 */
+      "direct_normal_broadband",	       /* fn=28 */
+      "qc_direct_normal_broadband",	       /* fn=29 */
+      "direct_normal_narrowband_filter1",      /* fn=30 */
+      "qc_direct_normal_narrowband_filter1",   /* fn=31 */
+      "direct_normal_narrowband_filter2",      /* fn=32 */
+      "qc_direct_normal_narrowband_filter2",   /* fn=33 */
+      "direct_normal_narrowband_filter3",      /* fn=34 */
+      "qc_direct_normal_narrowband_filter3",   /* fn=35 */
+      "direct_normal_narrowband_filter4",      /* fn=36 */
+      "qc_direct_normal_narrowband_filter4",   /* fn=37 */
+      "direct_normal_narrowband_filter5",      /* fn=38 */
+      "qc_direct_normal_narrowband_filter5",   /* fn=39 */
+      "direct_normal_narrowband_filter6",      /* fn=40 */
+      "qc_direct_normal_narrowband_filter6",   /* fn=41 */
+      "alltime_hemisp_broadband",	       /* fn=42 */
+      "qc_alltime_hemisp_broadband",           /* fn=43 */
+      "alltime_hemisp_narrowband_filter1",     /* fn=44 */
+      "qc_alltime_hemisp_narrowband_filter1",  /* fn=45 */
+      "alltime_hemisp_narrowband_filter2",     /* fn=46 */
+      "qc_alltime_hemisp_narrowband_filter2",  /* fn=47 */
+      "alltime_hemisp_narrowband_filter3",     /* fn=48 */
+      "qc_alltime_hemisp_narrowband_filter3",  /* fn=49 */
+      "alltime_hemisp_narrowband_filter4",     /* fn=50 */
+      "qc_alltime_hemisp_narrowband_filter4",  /* fn=51 */
+      "alltime_hemisp_narrowband_filter5",     /* fn=52 */
+      "qc_alltime_hemisp_narrowband_filter5",  /* fn=53 */
+      "alltime_hemisp_narrowband_filter6",     /* fn=54 */
+      "qc_alltime_hemisp_narrowband_filter6",  /* fn=55 */
+      "direct_horizontal_broadband",	       /* fn=56 */
+      "qc_direct_horizontal_broadband",	       /* fn=57 */
+      "direct_horizontal_narrowband_filter1",  /* fn=58 */
+      "qc_direct_horizontal_narrowband_filter1"/* fn=59 */,
+      "direct_horizontal_narrowband_filter2",   /* fn=60 */
+      "qc_direct_horizontal_narrowband_filter2",/* fn=61 */
+      "direct_horizontal_narrowband_filter3",	/* fn=62 */
+      "qc_direct_horizontal_narrowband_filter3",/* fn=63 */
+      "direct_horizontal_narrowband_filter4",	/* fn=64 */
+      "qc_direct_horizontal_narrowband_filter4",/* fn=65 */
+      "direct_horizontal_narrowband_filter5",	/* fn=66 */
+      "qc_direct_horizontal_narrowband_filter5",/* fn=67 */
+      "direct_horizontal_narrowband_filter6",	/* fn=68 */
+      "qc_direct_horizontal_narrowband_filter6",/* fn=69 */
+      "direct_diffuse_ratio_broadband",         /* fn=70 */
+      "qc_direct_diffuse_ratio_broadband",      /* fn=71 */
+      "direct_diffuse_ratio_filter1",		/* fn=72 */
+      "qc_direct_diffuse_ratio_filter1",	/* fn=73 */
+      "direct_diffuse_ratio_filter2",		/* fn=74 */
+      "qc_direct_diffuse_ratio_filter2",	/* fn=75 */
+      "direct_diffuse_ratio_filter3",		/* fn=76 */
+      "qc_direct_diffuse_ratio_filter3",	/* fn=77 */
+      "direct_diffuse_ratio_filter4",		/* fn=78 */
+      "qc_direct_diffuse_ratio_filter4",	/* fn=79 */
+      "direct_diffuse_ratio_filter5",           /* fn=80 */
+      "qc_direct_diffuse_ratio_filter5",        /* fn=81 */
+      "direct_diffuse_ratio_filter6",		/* fn=82 */
+      "qc_direct_diffuse_ratio_filter6",	/* fn=83 */
+      "head_temp",				/* fn=84 */
+      "qc_head_temp",				/* fn=85 */
+      "head_temp2",				/* fn=86 */
+      "qc_head_temp2",				/* fn=87 */
+      "logger_temperature",	        	/* fn=88 */
+      "qc_logger_temperature",   		/* fn=89 */
+      "logger_volt",                            /* fn=90 */
+      "qc_logger_volt",                         /* fn=91 */
+      "solar_zenith_angle",			/* fn=92 */
+      "cosine_solar_zenith_angle",		/* fn=93 */
+      "elevation_angle",			/* fn=94 */
+      "airmass",				/* fn=95 */
+      "azimuth_angle",				/* fn=96 */
+      "computed_cosine_correction_broadband",	/* fn=97 */
+      "computed_cosine_correction_filter1",	/* fn=98 */
+      "computed_cosine_correction_filter2",	/* fn=99 */
+      "computed_cosine_correction_filter3",     /* fn=100 */
+      "computed_cosine_correction_filter4",     /* fn=101 */
+      "computed_cosine_correction_filter5",	/* fn=102 */
+      "computed_cosine_correction_filter6",	/* fn=103 */
+      "bench_angle",				/* fn=104 */
+      "cosine_correction_sn_broadband",		/* fn=105 */
+      "cosine_correction_sn_filter1",		/* fn=106 */
+      "cosine_correction_sn_filter2",		/* fn=107 */
+      "cosine_correction_sn_filter3",		/* fn=108 */
+      "cosine_correction_sn_filter4",		/* fn=109 */
+      "cosine_correction_sn_filter5",           /* fn=110 */
+      "cosine_correction_sn_filter6",           /* fn=111 */
+      "cosine_correction_we_broadband",		/* fn=112 */
+      "cosine_correction_we_filter1",		/* fn=113 */
+      "cosine_correction_we_filter2",		/* fn=114 */
+      "cosine_correction_we_filter3",		/* fn=115 */
+      "cosine_correction_we_filter4",		/* fn=116 */
+      "cosine_correction_we_filter5",		/* fn=117 */
+      "cosine_correction_we_filter6",		/* fn=118 */
+      "wavelength_filter1",			/* fn=119 */
+      "normalized_transmittance_filter1",       /* fn=120 */
+      "wavelength_filter2",                     /* fn=121 */
+      "normalized_transmittance_filter2",	/* fn=122 */
+      "wavelength_filter3",			/* fn=123 */
+      "normalized_transmittance_filter3",	/* fn=124 */
+      "wavelength_filter4",			/* fn=125 */
+      "normalized_transmittance_filter4",	/* fn=126 */
+      "wavelength_filter5",			/* fn=127 */
+      "normalized_transmittance_filter5",	/* fn=128 */
+      "wavelength_filter6",			/* fn=129 */
+      "normalized_transmittance_filter6",       /* fn=130 */
+      "offset_broadband",                       /* fn=131 -- now time varying */  
+      "offset_filter1",				/* fn=132 -- now time varying */
+      "offset_filter2",				/* fn=133 -- now time varying */
+      "offset_filter3",				/* fn=134 -- now time varying */
+      "offset_filter4",				/* fn=135 -- now time varying */
+      "offset_filter5",				/* fn=136 -- now time varying */
+      "offset_filter6",				/* fn=137 -- now time varying */
+      "diffuse_correction_broadband",		/* fn=138 */
+      "diffuse_correction_filter1",		/* fn=139 */
+      "diffuse_correction_filter2",             /* fn=140 */
+      "diffuse_correction_filter3",             /* fn=141 */
+      "diffuse_correction_filter4",		/* fn=142 */
+      "diffuse_correction_filter5",		/* fn=143 */
+      "diffuse_correction_filter6",		/* fn=144 */
+      "nominal_calibration_factor_broadband",	/* fn=145 */
+      "nominal_calibration_factor_filter1",	/* fn=146 */
+      "nominal_calibration_factor_filter2",	/* fn=147 */
+      "nominal_calibration_factor_filter3",	/* fn=148 */
+      "nominal_calibration_factor_filter4",	/* fn=149 */
+      "nominal_calibration_factor_filter5",     /* fn=150 */
+      "nominal_calibration_factor_filter6",     /* fn=151 */
+      "lat_b1",                                 /* fn=152 */
+      "lon_b1",                                 /* fn=153 */
+      "alt_b1",                                 /* fn=154 */
+      NULL
+    },
+    {
+      //nimfr.a0
+      "direct_normal_broadband_raw",
+      "direct_normal_narrowband_filter1_raw",
+      "direct_normal_narrowband_filter2_raw",
+      "direct_normal_narrowband_filter3_raw",
+      "direct_normal_narrowband_filter4_raw",
+      "direct_normal_narrowband_filter5_raw",
+      "direct_normal_narrowband_filter6_raw",
+      "direct_normal_broadband",
+      "direct_normal_narrowband_filter1",
+      "direct_normal_narrowband_filter2",
+      "direct_normal_narrowband_filter3",
+      "direct_normal_narrowband_filter4",
+      "direct_normal_narrowband_filter5",
+      "direct_normal_narrowband_filter6",
+      "lat_a0",
+      "lon_a0",
+      "alt_a0",
+      NULL
+    },
+    {
+      //nimfr.b1
+      "direct_normal_broadband_b1",                      /* fn=0 */   
+      "qc_direct_normal_broadband_b1",			 /* fn=1 */   
+      "direct_normal_narrowband_filter1_b1",		 /* fn=2 */   
+      "qc_direct_normal_narrowband_filter1_b1",		 /* fn=3 */   
+      "direct_normal_narrowband_filter2_b1",		 /* fn=4 */   
+      "qc_direct_normal_narrowband_filter2_b1",		 /* fn=5 */   
+      "direct_normal_narrowband_filter3_b1",		 /* fn=6 */   
+      "qc_direct_normal_narrowband_filter3_b1",		 /* fn=7 */   
+      "direct_normal_narrowband_filter4_b1",		 /* fn=8 */   
+      "qc_direct_normal_narrowband_filter4_b1",		 /* fn=9 */   
+      "direct_normal_narrowband_filter5_b1",		 /* fn=10 */  
+      "qc_direct_normal_narrowband_filter5_b1",		 /* fn=11 */  
+      "direct_normal_narrowband_filter6_b1",		 /* fn=12 */  
+      "qc_direct_normal_narrowband_filter6_b1",		 /* fn=13 */  
+      "direct_horizontal_broadband",			 /* fn=14 */  
+      "qc_direct_horizontal_broadband",			 /* fn=15 */  
+      "direct_horizontal_narrowband_filter1",		 /* fn=16 */  
+      "qc_direct_horizontal_narrowband_filter1",	 /* fn=17 */  
+      "direct_horizontal_narrowband_filter2",		 /* fn=18 */  
+      "qc_direct_horizontal_narrowband_filter2",	 /* fn=19 */  
+      "direct_horizontal_narrowband_filter3",		 /* fn=20 */  
+      "qc_direct_horizontal_narrowband_filter3",	 /* fn=21 */  
+      "direct_horizontal_narrowband_filter4",		 /* fn=22 */  
+      "qc_direct_horizontal_narrowband_filter4",	 /* fn=23 */  
+      "direct_horizontal_narrowband_filter5",		 /* fn=24 */  
+      "qc_direct_horizontal_narrowband_filter5",	 /* fn=25 */  
+      "direct_horizontal_narrowband_filter6",		 /* fn=26 */  
+      "qc_direct_horizontal_narrowband_filter6",	 /* fn=27 */  
+      "head_temp",					 /* fn=28 */  
+      "qc_head_temp",					 /* fn=29 */  
+      "logger_volt",					 /* fn=30 */  
+      "qc_logger_volt",					 /* fn=31 */  
+      "tube_temp",					 /* fn=32 */  
+      "qc_tube_temp",					 /* fn=33 */  
+      "head_temp2",					 /* fn=34 */  
+      "qc_head_temp2",					 /* fn=35 */  
+      "solar_zenith_angle",				 /* fn=36 */  
+      "cosine_solar_zenith_angle",			 /* fn=37 */  
+      "elevation_angle",				 /* fn=38 */  
+      "airmass",					 /* fn=39 */  
+      "azimuth_angle",					 /* fn=40 */  
+      "wavelength_filter1",				 /* fn=41 */  
+      "normalized_transmittance_filter1",		 /* fn=42 */  
+      "wavelength_filter2",				 /* fn=43 */  
+      "normalized_transmittance_filter2",		 /* fn=44 */  
+      "wavelength_filter3",				 /* fn=45 */  
+      "normalized_transmittance_filter3",		 /* fn=46 */  
+      "wavelength_filter4",				 /* fn=47 */  
+      "normalized_transmittance_filter4",		 /* fn=48 */  
+      "wavelength_filter5",				 /* fn=49 */  
+      "normalized_transmittance_filter5",		 /* fn=50 */  
+      "wavelength_filter6",				 /* fn=51 */  
+      "normalized_transmittance_filter6",		 /* fn=52 */  
+      "offset_broadband",				 /* fn=53 */  
+      "offset_filter1",					 /* fn=54 */  
+      "offset_filter2",					 /* fn=55 */  
+      "offset_filter3",					 /* fn=56 */  
+      "offset_filter4",					 /* fn=57 */  
+      "offset_filter5",					 /* fn=58 */  
+      "offset_filter6",					 /* fn=59 */ 
+      "nominal_calibration_factor_broadband",		 /* fn=60 */ 
+      "nominal_calibration_factor_filter1",		 /* fn=61 */ 
+      "nominal_calibration_factor_filter2",		 /* fn=62 */ 
+      "nominal_calibration_factor_filter3",		 /* fn=63 */ 
+      "nominal_calibration_factor_filter4",		 /* fn=64 */ 
+      "nominal_calibration_factor_filter5",		 /* fn=65 */ 
+      "nominal_calibration_factor_filter6",		 /* fn=66 */ 
+      "lat_b1",						 /* fn=67 */ 
+      "lon_b1",						 /* fn=68 */ 
+      "alt_b1",						 /* fn=69 */ 
+      NULL			
+    },				
+    {				
+      //mfrsrlangley.c1					 
+      "barnard_optical_depth_broadband",		 
+      "barnard_optical_depth_filter1",			 
+      "barnard_optical_depth_filter2",			 
+      "barnard_optical_depth_filter3",			 
+      "barnard_optical_depth_filter4",			 
+      "barnard_optical_depth_filter5",			 
+      "barnard_optical_depth_filter6",			 
+      "michalsky_optical_depth_broadband",		 
+      "michalsky_optical_depth_filter1",		 
+      "michalsky_optical_depth_filter2",		 
+      "michalsky_optical_depth_filter3",		 
+      "michalsky_optical_depth_filter4",		 
+      "michalsky_optical_depth_filter5",		 
+      "michalsky_optical_depth_filter6",		 
+      "barnard_solar_constant_sdist_broadband",		 
+      "barnard_solar_constant_sdist_filter1",		 
+      "barnard_solar_constant_sdist_filter2",		 
+      "barnard_solar_constant_sdist_filter3",		 
+      "barnard_solar_constant_sdist_filter4",		 
+      "barnard_solar_constant_sdist_filter5",		 
+      "barnard_solar_constant_sdist_filter6",		 
+      "michalsky_solar_constant_sdist_broadband",	 
+      "michalsky_solar_constant_sdist_filter1",		 
+      "michalsky_solar_constant_sdist_filter2",		 
+      "michalsky_solar_constant_sdist_filter3",		 
+      "michalsky_solar_constant_sdist_filter4",		 
+      "michalsky_solar_constant_sdist_filter5",		 
+      "michalsky_solar_constant_sdist_filter6",		 
+      "barnard_badflag_broadband",		         
+      "barnard_badflag_filter1",    		         
+      "barnard_badflag_filter2",        		 
+      "barnard_badflag_filter3",	        	 
+      "barnard_badflag_filter4",		         
+      "barnard_badflag_filter5",		         
+      "barnard_badflag_filter6",		         
+      "michalsky_badflag_broadband",		         
+      "michalsky_badflag_filter1",		         
+      "michalsky_badflag_filter2",        		 
+      "michalsky_badflag_filter3",	        	 
+      "michalsky_badflag_filter4",		         
+      "michalsky_badflag_filter5",	        	 
+      "michalsky_badflag_filter6",		         
+      "lat_langley",
+      "lon_langley",
+      "alt_langley",
+      NULL
+    },
+    {
+      //omi.a1
+      "ozone_omi",
+      "lat_omi",
+      "lon_omi",
+      NULL
+    },
+    {
+      //toms.a1
+      "ozone_toms",
+      "lat_toms",
+      "lon_toms",
+      NULL
+    },
+    {
+      //met.b1
+      "atmos_pressure",
+      NULL
+    },
+  };
+
+#endif // NO_RETRIEVER_H
+
+#endif // _MFROD1BARNMICH_RETRIEVER_H

src/sw.h

+#ifndef H_SW
+#define H_SW
+
+#ifndef TRUE
+#define TRUE 1
+#define FALSE 0
+#endif
+
+typedef enum 
+{
+  northern, 
+  southern
+} hemis;
+
+typedef enum 
+{
+  tauonly, 
+  tauio, 
+  taur, 
+  tauior
+} output_type;
+
+typedef struct 
+{
+  FILE	*iofp, *infp, *outfp, *outintfp;	/* file pointers for Ios and output */
+  char	*iofname, *infname, *outfname, *outintfname; 
+  char 	**infilelist;
+  double        am[8640];		/* airmass vector [numsamples * 2]*/
+
+  double        ios[7][366];		/* Ios vector for "epoch" ... these are output of Forgan analysis */
+  double        iosTimes[366];
+
+  double	direct[7][8640];	/* pointer into **data for direct columns*/
+  double	dataTimes[8640];	/* vector of times for current file */
+
+  double        taus[8640];			/* taus for 'key' wavelength */
+  double	currIos[7];		/* Ios for current day */
+  double        logCurrIos[7];		/* log(currIos) */
+  double	soldst;		        /* current day's solar distance */
+  double	tol;			/* tolerance for lsfit */
+  double	longitude;
+  double        latitude;
+  double        azimuthSolarNoon;	/* what's the azimuth of solar noon ? */
+  double	minAM;
+  double        maxAM;		        /* airmass limits */
+  long	        resultsWinBegin[8640];	/* lists of sliding window results */
+  long	        resultsWinEnd[8640];		
+  int	        am_begin;
+  int           am_end;           	/* indicies of data array corresp. to AM range */
+  int		amVectSize;		/* current memory allocation for airmass vector */
+  int		tauVectSize;		/* current memory allocation for tau vector */
+  int		resVectSize;		/* current memory allocation for results vector */
+  int		numTaus;		/* current number of taus = am_end-am_begin+1 */
+  int		numRes;			/* # of results (good windows) coming back */
+  int        	numinfiles;
+  int 	   	totalRowsIos;		/* # of rows in Ios file */
+  int		totalColsIos;		/* # of cols in Ios file */
+  int		totalRowsData;		/* # of rows in beam input file */
+  int		totalColsData;		/* # of cols in beam input file */
+  int       	sampling;		/* the sampling rate of the data, in seconds */
+  int	   	winsize;		/* sliding window size (in samples) */
+  int	   	winminutes;		/* sliding window size (in minutes) */
+  int		wincount;		/* total number of windows encountered */
+  int	   	dirColOffsets[7];	/* list of beam columns */
+  int 		ioColOffsets[7];	/* list of Io columns */
+  int		numDirColOffsets;	/* number beam columns */
+  int		numIoColOffsets;	/* number of Io columns */
+  int		key;			/* key wavelength col */
+  char 	        integrate;
+  char		debug; 
+  char		verbose;
+  hemis	        hemisphere;
+} sw_type;
+
+sw_type sw; 
+
+#endif

src/time_funcs.h

+#define ISLEAP(y) (y<1900 \
+                   ? (((y+1900) % 4) == 0 && ((y+1900) % 100) != 0 \
+                     || ((y+1900) % 400) == 0) \
+                   : (((y) % 4) == 0 && ((y) % 100) !=0 || ((y) % 400) == 0))
+
+#define DYSIZE(y) 365 + ISLEAP(y)
+
+#define DAYSTO1970 25568
+#define SECSINDAY  86400
+
+int vap_Rundate;