/**---------------------------------------------------------------------------*/
// tchart: A program to generate temperature charts
// See the help text for documentation.
// Copyright (c) 2009, Richard Heurtley. All rights reserved.
/**---------------------------------------------------------------------------*/
#include <ctype.h>
#include <limits.h>
#include <malloc.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "image.h"
#include "matrix.h"
#include "mem.h"

// Tell the kiss_fft package to use doubles instead of floats.

#define kiss_fft_scalar double

#include "kiss_fft.h"
/**---------------------------------------------------------------------------*/
#define HICOLOR		(LRED)
#define LOCOLOR		(LBLUE)

#define COLOR_RED	(1)
#define COLOR_BLUE	(2)

#define OPTION_TYPE_NONE	(0)
#define OPTION_TYPE_INT		(1)
#define OPTION_TYPE_STRING	(2)

typedef struct
{
 int row;
 int color;
} PIX;

typedef struct
{
 PIX *pix0p;
 PIX *pix1p;
} GAP;

// daily data

typedef struct
{
 double ahi;		// ((double *)dayp)[0]
 double alo;		// ((double *)dayp)[1]
 double sdhi;		// ((double *)dayp)[2]
 double sdlo;		// ((double *)dayp)[3]
 double xahi;		// ((double *)dayp)[0+4]
 double xalo;		// ((double *)dayp)[1+4]
 double xsdhi;		// ((double *)dayp)[2+4]
 double xsdlo;		// ((double *)dayp)[3+4]
 int month;
 int day;
 int left;
 int right;
 int mid;
} DAY;

typedef struct
{
 char *name;		// option name
 int type;		// option type
 int *specp;		// option specified: 0 = no, 1 = yes
 void *valuep;		// option value
} OPTION;

// options

typedef struct
{
 int debug;
 int quiet;
 int hispec;
 int hi;
 int lospec;
 int lo;
 int help;
 int dpispec;
 int dpi;
 int lsq;
 int nterms;
 int fft;
 int fftradius;
 int ave;
 int averadius;
 int med;
 int medradius;
 int o;
 char filename[64];	// output filename
} LOCAL;

static LOCAL l;
/**---------------------------------------------------------------------------*/
int main(int npars, char **pars)
{
 int rv = 0;
 int n;
 int month;
 int day;
 int fh;
 int fw;
 int left;
 int bottom;
 int right;
 int top;
 int ileft;
 int iright;
 int nicols;
 int itop;
 int ibottom;
 int nirows;
 int row;
 int row0;
 int ngaps;
 int npix;
 int pix;
 int itemp;
 int doy;
 int index;
 int i;
 int j;
 int best;
 int rc;
 int term;
 int col;
 int nrows;
 int ncols;
 int par;
 int opt;

 double hi;
 double lo;
 double bhi;
 double blo;
 double bdelta;
 double rdelta;
 double temp;
 double ahi;
 double alo;
 double precip;
 double sdhi;
 double sdlo;
 double nsamp1;
 double nsamp2;
 double nsamp3;
 double ddoy;
 double v;
 double x;

 char *p;
 char *q;
 char *parp;
 char *dotp;
 OPTION *optp;
 double *arrayp;
 XRGB xrgb;
 FILE *filep;
 IMAGE *ip;
 DAY *dayp;
 PIX *pix0p;
 PIX *pix1p;
 kiss_fft_cfg kfc;
 MATRIX *mp;
 double *cp;
 TERM *terms;
 TERM *termp;

 static GAP gaps[256];
 static PIX pixs[256];
 static DAY days[366];
 static XRGB xrgb_red  = {0, 255, 127, 127};
 static XRGB xrgb_blue = {0, 127, 127, 255};
 static XRGB xrgb_mag  = {0, 255, 127, 255};

 static char location[256];
 static char description[256];
 static char buffer[256];
 static char string[80];
 static char *months[12] =
 {
  "Jan",
  "Feb",
  "Mar",
  "Apr",
  "May",
  "Jun",
  "Jul",
  "Aug",
  "Sep",
  "Oct",
  "Nov",
  "Dec",
 };

 static OPTION options[] =
 {
  {"debug",	OPTION_TYPE_NONE,	&l.debug,	0},
  {"q",		OPTION_TYPE_NONE,	&l.quiet,	0},
  {"help",	OPTION_TYPE_NONE,	&l.help,	0},
  {"o",		OPTION_TYPE_STRING,	&l.o,		&l.filename},
  {"dpi",	OPTION_TYPE_INT,	&l.dpispec,	&l.dpi},
  {"average",	OPTION_TYPE_INT,	&l.ave,		&l.averadius},
  {"median",	OPTION_TYPE_INT,	&l.med,		&l.medradius},
  {"lsq",	OPTION_TYPE_INT,	&l.lsq,		&l.nterms},
  {"fft",	OPTION_TYPE_INT,	&l.fft,		&l.fftradius},
  {"hi",	OPTION_TYPE_INT,	&l.hispec,	&l.hi},
  {"lo",	OPTION_TYPE_INT,	&l.lospec,	&l.lo},
 };

 static int noptions = sizeof(options) / sizeof(OPTION);

 static kiss_fft_cpx cx_in[366];
 static kiss_fft_cpx cx_out[366];
/**---------------------------------------------------------------------------*/
// initialize handles and pointers

 ip = 0;		// image
 mp = 0;		// matrix
 filep = 0;		// file
 terms = 0;		// terms array
 arrayp = 0;		// median array

 if (!l.quiet)
  printf("\n");
/**---------------------------------------------------------------------------*/
// parse parameters

// if no parameters then render help

 if (npars == 1)
 {
  l.help = 1;
  goto HELP;
 }

// process parameters

 for (par = 1; par < npars; par++)
 {
  parp = pars[par];
 
  if (parp[0] != '-')
  {
   if (filep)
   {
    printf("Error: More than one file specified.\n");
    rv = 1;
    goto EGRESS;
   }

   strcpy(string, parp);
   dotp = strchr(string, '.');
   if (!dotp)
    strcat(string, ".html");

   filep = fopen(string, "r");
   if (!filep)
   {
    printf("Error: Unable to open file \"%s\".\n", string);
    rv = 1;
    goto EGRESS;
   }
   if (!l.quiet)
    printf("Reading file \"%s\".\n", string);
   continue;
  }

  for (opt = 0; opt < noptions; opt++)
  {
   optp = &options[opt];
   if (!strcmp(optp->name, &parp[1]))
    break;
  }
  if (opt >= noptions)
  {
   printf("Error: \"%s\" is not a parameter.", parp);
   rv = 1;
   goto EGRESS;
  }

  if (optp->specp)
   *(optp->specp) = 1;		// note parameter specified

  switch (optp->type)
  {
   case OPTION_TYPE_NONE:
   {
    if (!l.quiet)
     printf("Option -%s specified.\n", optp->name);
    break;			// don't do anything
   }
   case OPTION_TYPE_INT:
   {
    par++;
    if (par == npars)		// dangling parameter
    {
     printf("Error: No value for \"%s\".", parp);
     rv = 1;
     goto EGRESS;
    }

    parp = pars[par];
    n = sscanf(parp, "%d", (int *)optp->valuep);
    if (n != 1)
    {
     printf("Error: Unable to parse \"%s\".", parp);
     rv = 1;
     goto EGRESS;
    }
    if (!l.quiet)
     printf("Option -%s %d specified.\n", optp->name, *(int *)optp->valuep);
    break;
   }
   case OPTION_TYPE_STRING:
   {
    par++;
    if (par == npars)		// dangling parameter
    {
     printf("Error: No value for \"%s\".", parp);
     rv = 1;
     goto EGRESS;
    }

    parp = pars[par];
    strcpy(optp->valuep, parp);

    if (!l.quiet)
     printf("Option -%s %s specified.\n", optp->name, (char *)optp->valuep);

    break;
   }
  }

 }		// for (par = 1; par < npars; par++)
/**---------------------------------------------------------------------------*/
HELP:

 if (l.help)
 {
  printf
  (
   "tchart: A program to generate temperature charts.\n"
   "\n"
   "Usage: tchart filename{.html} [options]\n"
   "\n"
   "If the -o option is not specified then tchart generates a\n"
   "file with the same name as the HTML input file but with a\n"
   ".BMP filetype.\n"
   "\n"
   "The chart is landscape 7.5 inches by 10.0 inches. The\n"
   "default DPI is 150.\n"
   "\n"
   "Multiple processing steps may be specified. Processing is\n"
   "done in this order: averaging, median filtering, polynomial\n"
   "fitting, FFT filtering.\n"
   "\n"
   "Options:\n"
   "\n"
   "-average radius\n"
   "\n"
   "Smooths the input data by replacing each value with the\n"
   "average of the values plus and minus the specified radius\n"
   "from the center value. A radius of 1 means that three points\n"
   "are averaged. A radius of 2 means that five points are\n"
   "averaged.\n"
   "\n"
   "-debug\n"
   "\n"
   "If the -fft option is specified then the -debug option logs\n"
   "the FFT coefficients to file tchart.txt.\n"
   "\n"
   "-dpi n\n"
   "\n"
   "Sets the chart's resolution. The default DPI is 150.\n"
   "\n"
   "-fft nterms\n"
   "\n"
   "Smooths the input data by calculating a Fast Fourier\n"
   "Transform (FFT), zeroing high frequency terms, and\n"
   "calculating the inverse FFT. The \"nterms\" parameter\n"
   "specifies the number of low frequency terms to keep. A value\n"
   "of 1 keeps just the first term which is the average value. A\n"
   "value of 2 keeps the first two terms which are the average\n"
   "value plus the full-width wave function.\n"
   "\n"
   "-help\n"
   "\n"
   "Displays this help message.\n"
   "\n"
   "-hi temp\n"
   "-lo temp\n"
   "\n"
   "By default tchart calculates the Y axis temperature range\n"
   "automatically. The Y axis high and low temperatures can be\n"
   "be specified with the -hi and -lo options. The high and low\n"
   "temperatures are integer values.\n"
   "\n"
   "-lsq nterms\n"
   "\n"
   "Smooths the input data by calculating an optimal\n"
   "least-squares polynomial equation. Parameter \"nterms\" is the\n"
   "number of polynomial terms. A value of 3 fits a three term\n"
   "parabola with equation v = ax^2 + bx + c. A value of 4 fits\n"
   "four term curve with equation v = ax^3 + bx^2 + cx + d.\n"
   "\n"
   "-median radius\n"
   "\n"
   "Smooths the input data by replacing each value with the\n"
   "median of the values plus and minus the specified radius\n"
   "from the center value. A radius of 1 means that three points\n"
   "are sorted to determine the median. A radius of 2 means that\n"
   "five points are sorted.\n"
   "\n"
   "-o filename{.bmp}\n"
   "\n"
   "Specifies the output filename. The default output filename\n"
   "is the input filename with a .BMP filetype.\n"
   "\n"
   "-q\n"
   "\n"
   "Quiet mode. Don't explain what's going on.\n"
   "\n"
  );
  rv = 0;
  goto EGRESS;
 }
/**---------------------------------------------------------------------------*/
// input file sample

#if 0

<html><head>
<meta http-equiv="content-type" content="text/html; charset=ISO-8859-1"><title> ENOSBURG FALLS, VERMONT 30 Year Daily Summary ,1971-2000 (WRCC)</title>
</head><body bgcolor="#ffffff"><h1> ENOSBURG FALLS, VERMONT </h1>
<h3>30 Year Daily Temperature and Precipitation Summary </h3>
<pre> STATION 432769 AVERAGES FROM AVAILABLE YEARS IN PERIOD 1971 TO 2000 .

 DOY  MON DY  TMAX  #YRS  TMIN  #YRS PRECIP #YRS SD MAX SD MIN
    1  1  1   28.4  29.    7.3  29.  0.090  28. 11.921 16.210

#endif
/**---------------------------------------------------------------------------*/
// get location

 if (!filep)
 {
  printf("Error: No file specified.\n");
  rv = 1;
  goto EGRESS;
 }

 if (!l.quiet)
  printf("Parsing station location.\n");

 for (;;)
 {
  p = fgets(buffer, sizeof(buffer), filep);
  if (!p)
   break;
  p = strstr(buffer, "<h1>");
  if (p)
   break;
 }
 if (!p)
 {
  printf("Error: No location.\n");
  rv = 1;
  goto EGRESS;
 }

 q = strstr(buffer, "</h1>");
 if (!q)
 {
  printf("Error: No end of header.\n");
  rv = 1;
  goto EGRESS;
 }

 *q = 0;
 strcpy(location, p+4);
/**---------------------------------------------------------------------------*/
// get description

 if (!l.quiet)
  printf("Parsing data description location.\n");

 for (;;)
 {
  p = fgets(buffer, sizeof(buffer), filep);
  if (!p)
   break;
  p = strstr(buffer, "<pre>");
  if (p)
   break;
 }
 if (!p)
 {
  printf("Error: No description.\n");
  rv = 1;
  goto EGRESS;
 }

 strcpy(description, p+5);
/**---------------------------------------------------------------------------*/
// get data

 if (!l.quiet)
  printf("Parsing data.\n");

 for (;;)
 {
  p = fgets(buffer, sizeof(buffer), filep);
  if (!p)
   break;

  n = sscanf
  (
   buffer,
   "%d %d %d %lf %lf %lf %lf %lf %lf %lf %lf\n",
   &doy, &month, &day,
   &ahi, &nsamp1,
   &alo, &nsamp2,
   &precip, &nsamp3,
   &sdhi, &sdlo
  );
  if (n != 11)
   continue;

  dayp = &days[doy-1];
  dayp->month = month - 1;		// zero based
  dayp->day = day;
  dayp->ahi = ahi;
  dayp->alo = alo;
  dayp->sdhi = sdhi;
  dayp->sdlo = sdlo;
 }

 fclose(filep);
 filep = 0;
/**---------------------------------------------------------------------------*/
// average data if requested

 if (!l.ave)
  goto NOAVERAGE;

 if (!l.quiet)
  printf("Averaging data.\n");

 for (doy = 0; doy < 366; doy++)
 {
  int i;
  double dn;
  double sumalo;
  double sumahi;
  double sumsdlo;
  double sumsdhi;

  sumalo = 0.0;
  sumahi = 0.0;
  sumsdlo = 0.0;
  sumsdhi = 0.0;
  dn = 0.0;
  for (i = doy - l.averadius; i <= doy + l.averadius; i++)
  {
   if (i < 0 || 366 <= i)
    continue;
   dayp = &days[i];
   sumalo += dayp->alo;
   sumahi += dayp->ahi;
   sumsdlo += dayp->sdlo;
   sumsdhi += dayp->sdhi;
   dn++;
  }
  dayp = &days[doy];
  dayp->xalo = sumalo / dn;
  dayp->xahi = sumahi / dn;
  dayp->xsdlo = sumsdlo / dn;
  dayp->xsdhi = sumsdhi / dn;
 }

 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];
  dayp->alo = dayp->xalo;
  dayp->ahi = dayp->xahi;
  dayp->sdlo = dayp->xsdlo;
  dayp->sdhi = dayp->xsdhi;
 }

NOAVERAGE:
/**---------------------------------------------------------------------------*/
// median filter if requested

 if (!l.med)
  goto NOMEDIAN;

 if (!l.quiet)
  printf("Median filtering data.\n");

 arrayp = MemGet((l.medradius+1+l.medradius) * sizeof(double));
 if (!arrayp)
 {
  printf("Error: Unable to allocate median array.\n");
  rv = 1;
  goto EGRESS;
 }

 for (index = 0; index < 4; index++)
 {
  for (doy = 0; doy < 366; doy++)
  {

// get range of samples

   n = 0;
   for (i = doy - l.medradius; i <= doy + l.medradius; i++)
   {
    if (i < 0 || 366 <= i)
     continue;
    dayp = &days[i];
    arrayp[n++] = ((double *)dayp)[index];
   }

// sort range (could do just first half)

   for (i = 0; i < n-1; i++)
   {
    lo = arrayp[i];
    best = i;
    for (j = i + 1; j < n; j++)
    {
     if (arrayp[j] < lo)
     {
      lo = arrayp[j];
      best = j;
     }
    }
    arrayp[best] = arrayp[i];
    arrayp[i] = lo;
   }

// get median value

   dayp = &days[doy];
   ((double *)dayp)[index+4] = arrayp[n/2];

  }
 }

 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];
  dayp->alo = dayp->xalo;
  dayp->ahi = dayp->xahi;
  dayp->sdlo = dayp->xsdlo;
  dayp->sdhi = dayp->xsdhi;
 }

 MemFree(arrayp);
 arrayp = 0;

NOMEDIAN:
/**---------------------------------------------------------------------------*/
// polynomial curve fitting if requested
// include 30 days preceeding and subsequent to the year

 if (!l.lsq)
  goto NOLSQ;

 if (!l.quiet)
  printf("Fitting polynomial equation.\n");

 mp = MatrixGet(l.nterms);
 if (!mp)
 {
  printf("Error: Unable to allocate matrix.\n");
  rv = 1;
  goto EGRESS;
 }

 terms = MemGet((l.nterms+1) * sizeof(TERM));
 if (!terms)
 {
  printf("Error: Unable to allocate terms array.\n");
  rv = 1;
  goto EGRESS;
 }

 for (index = 0; index < 4; index++)
 {
  MatrixZero(mp);

  for (doy = -30; doy < 366 + 30; doy++)
  {
   i = doy;
   if (i < 0)		// -1 maps to 365
    i += 366;
   if (i >= 366)	// 366 maps to 0
    i -= 366;

   dayp = &days[i];

   ddoy = (double)doy - 183.0;		// the middle of the year is 0.0

// add doy:value sample to least squares matrix

   terms[0].i = 0;
   terms[0].v = 1.0;

   for (term = 1; term < l.nterms; term++)
   {
    termp = &terms[term];
    terms[term].i = term;
    terms[term].v = terms[term-1].v * ddoy;
   }

   terms[l.nterms].i = l.nterms;
   terms[l.nterms].v = ((double *)dayp)[index];

   MatrixHang(mp, l.nterms+1, terms, 0, 0);
  }

// solve matrix

  rc = MatrixReduce(mp, EXACT);
  if (rc != -1)
  {
   printf("Error: Unable to solve matrix.\n");
   rv = 1;
   goto EGRESS;
  }

  cp = MatrixCoeffs(mp);		// get calculated coefficients

// calculate curves

  for (doy = 0; doy < 366; doy++)
  {
   dayp = &days[doy];

   x = (double)doy - 183.0;

   v = 0.0;
   for (term = l.nterms-1; term > 0; term--)
   {
    v += cp[term];
    v *= x;
   }
   v += cp[0];

   ((double *)dayp)[index] = v;
  }

 }

 MemFree(terms);
 terms = 0;

 MatrixFree(mp);
 mp = 0;

NOLSQ:
/**---------------------------------------------------------------------------*/
// FFT curve fitting if requested

 if (!l.fft)
  goto NOFFT;

 if (!l.quiet)
  printf("FFT low frequency filtering.\n");

 if (l.debug)
  filep = fopen("tchart.txt", "w");

 l.fftradius = 183 - l.fftradius;	// nretained to nzeroed

 for (index = 0; index < 4; index++)
 {
  for (doy = 0; doy < 366; doy++)
  {
   dayp = &days[doy];
   cx_in[doy].r = ((double *)dayp)[index];
   cx_in[doy].i = 0.0;
  }

// do forward transform

  kfc = kiss_fft_alloc(366, 0, 0, 0);
  kiss_fft(kfc, cx_in, cx_out);
  free(kfc);

// log FFT values

  if (l.debug)
  {
   fprintf(filep, "----------------------------------------\n");

   doy = 183;
   fprintf
   (
    filep,
    "%3i %6.2f %6.2f\n",
    doy,
    log10(fabs(cx_out[doy].r)), 
    log10(fabs(cx_out[doy].i))
   );

   for (i = 0; i <= 182; i++)
   {
    doy = 182 - i;
    fprintf
    (
     filep,
     "%3i %6.2f %6.2f",
     doy,
     log10(fabs(cx_out[doy].r)), 
     log10(fabs(cx_out[doy].i))
    );

    doy = 184 + i;
    if (doy < 366)
    {
     fprintf(filep, "  -  ");
     fprintf
     (
      filep,
      "%3i %6.2f %6.2f",
      doy,
      log10(fabs(cx_out[doy].r)), 
      log10(fabs(cx_out[doy].i))
     );
    }

    fprintf(filep, "\n");
   }
  }

// zero high frequency terms to smooth curve

  for (doy = 183 - l.fftradius; doy <= 183 + l.fftradius; doy++)
  {
   cx_out[doy].r = 0.0;
   cx_out[doy].i = 0.0;
  }

// get complex conjugate and scale

  for (doy = 0; doy < 366; doy++)
  {
   cx_in[doy].r =  cx_out[doy].r / 366.0;
   cx_in[doy].i = -cx_out[doy].i / 366.0;
  }

// do reverse transform

  kfc = kiss_fft_alloc(366, 0, 0, 0);
  kiss_fft(kfc, cx_in, cx_out);
  free(kfc);

// get complex conjugate and scale

  for (doy = 0; doy < 366; doy++)
  {
   cx_in[doy].r =  cx_out[doy].r / 366.0;
   cx_in[doy].i = -cx_out[doy].i / 366.0;
  }

// get real result

  for (doy = 0; doy < 366; doy++)
  {
   dayp = &days[doy];
   ((double *)dayp)[index] = cx_out[doy].r;
  }

 }

 if (l.debug)
 {
  fclose(filep);
  filep = 0;
 }

NOFFT:
/**---------------------------------------------------------------------------*/
// allocate image memory

 if (!l.quiet)
  printf("Rendering chart.\n");

 if (!l.dpispec)			// if no dpi specified
  l.dpi = 150;				// get default

 nrows = ( 75 * l.dpi + 5) / 10;	// 7.5" vertically
 ncols = (100 * l.dpi + 5) / 10;	// 10" horizontally

 rv = ImageInit
 (
  &ip,		// structure
  nrows,	// height
  ncols,	// width
  ncols * 3,	// pitch
  3,		// nbpp
  0xff0000,	// red (doesn't really matter)
  0x00ff00,	// green "
  0x0000ff,	// blue  "
  0		// datap
 );
 if (rv)
 {
  printf("Error: Unable to allocate image.\n");
  rv = 1;
  goto EGRESS;
 }
/**---------------------------------------------------------------------------*/
// setup and determine dimensions

 ImageRadiusSet(ip, 0);			// line width
 ImageMagSet(ip, 1);			// font size
 ImageColorSet(ip, BWHITE);		// background color
 ImageBlank(ip);			// blank image

 fh = ImageFontHeightGet(ip);
 fw = ImageFontWidthGet(ip);

// determine box

 left = (0) + (fw * 3 + 1 + 1);		// "-10" (blank) (tic)
 bottom = (nrows-1) - (3 * fh + 1 + 1);	// station, title, "Jan" (blank) (tic)
 right = (ncols-1) - (1);		// (tic)
 top = (0) + ((fh+1)/2 + 1);		// (1/2 font height) (tic)

// determine box interior

 ileft = left + 1;
 iright = right - 1;
 nicols = iright - ileft + 1;
 itop = top + 1;
 ibottom = bottom - 1;
 nirows = ibottom - itop + 1;
/**---------------------------------------------------------------------------*/
// determine day columns

 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];
  dayp->left  = ileft + ((double)(doy+0) * (double)nicols / 366.0 + 0.5);
  dayp->right = ileft + ((double)(doy+1) * (double)nicols / 366.0 + 0.5);
  dayp->right--;
  dayp->mid = (dayp->left + dayp->right + 1) / 2;
 }
/**---------------------------------------------------------------------------*/
// determine temperature rows

 dayp = &days[0];
 hi = dayp->ahi + dayp->sdhi;
 lo = dayp->alo - dayp->sdlo;
 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];
  if (dayp->ahi + dayp->sdhi > hi)
   hi = dayp->ahi + dayp->sdhi;
  if (dayp->alo - dayp->sdlo < lo)
   lo = dayp->alo - dayp->sdlo;
 }

 hi += 1000.0;
 bhi = (double)( (int)(hi / 5.0) );
 bhi *= 5.0;
 if (bhi != hi)
  bhi += 5.0;
 hi -= 1000.0;
 bhi -= 1000.0;

 lo += 1000.0;
 blo = (double)( (int)(lo / 5.0) );
 blo *= 5.0;
 lo -= 1000.0;
 blo -= 1000.0;

 if (l.hispec)
  bhi = (double)l.hi;

 if (l.lospec)
  blo = (double)l.lo;

 bdelta = bhi - blo;
 rdelta = (double)(ibottom - itop);
/**---------------------------------------------------------------------------*/
// draw box

 ImageColorSet(ip, BLACK);

 ImageMoveTo(ip, top, left);
 ImageLineTo(ip, top, right);
 ImageLineTo(ip, bottom, right);
 ImageLineTo(ip, bottom, left);
 ImageLineTo(ip, top, left);
/**---------------------------------------------------------------------------*/
// draw average high + sd

 ImageColorSet(ip, HICOLOR);

 dayp = &days[365];
 temp = dayp->ahi + dayp->sdhi;
 dayp = &days[0];
 temp += dayp->ahi + dayp->sdhi;
 temp /= 2.0;
 row0 = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
 ImageMoveTo(ip, row0, dayp->left);

 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];
  temp = dayp->ahi + dayp->sdhi;
  row = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
  ImageLineTo(ip, row, dayp->mid);
 }

 ImageLineTo(ip, row0, dayp->right);
/**---------------------------------------------------------------------------*/
// draw average high - sd

 ImageColorSet(ip, HICOLOR);

 dayp = &days[365];
 temp = dayp->ahi - dayp->sdhi;
 dayp = &days[0];
 temp += dayp->ahi - dayp->sdhi;
 temp /= 2.0;
 row0 = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
 ImageMoveTo(ip, row0, dayp->left);

 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];
  temp = dayp->ahi - dayp->sdhi;
  row = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
  ImageLineTo(ip, row, dayp->mid);
 }

 ImageLineTo(ip, row0, dayp->right);
/**---------------------------------------------------------------------------*/
// draw average low + sd

 ImageColorSet(ip, LOCOLOR);

 dayp = &days[365];
 temp = dayp->alo + dayp->sdlo;
 dayp = &days[0];
 temp += dayp->alo + dayp->sdlo;
 temp /= 2.0;
 row0 = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
 ImageMoveTo(ip, row0, dayp->left);

 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];
  temp = dayp->alo + dayp->sdlo;
  row = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
  ImageLineTo(ip, row, dayp->mid);
 }

 ImageLineTo(ip, row0, dayp->right);
/**---------------------------------------------------------------------------*/
// draw average low - sd

 ImageColorSet(ip, LOCOLOR);

 dayp = &days[365];
 temp = dayp->alo - dayp->sdlo;
 dayp = &days[0];
 temp += dayp->alo - dayp->sdlo;
 temp /= 2.0;
 row0 = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
 ImageMoveTo(ip, row0, dayp->left);

 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];
  temp = dayp->alo - dayp->sdlo;
  row = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
  ImageLineTo(ip, row, dayp->mid);
 }

 ImageLineTo(ip, row0, dayp->right);
/**---------------------------------------------------------------------------*/
// draw colors

 for (col = ileft; col <= iright; col++)
 {

// compile pixel groups

  npix = 0;
  for (row = itop; row <= ibottom; row++)
  {
   xrgb = ImagePixelGet(ip, row, col);
   if (!memcmp(&xrgb, &HICOLOR, sizeof(xrgb)))
   {
    pixs[npix].row = row;
    pixs[npix].color = COLOR_RED;
    npix++;
   }
   else
   if (!memcmp(&xrgb, &LOCOLOR, sizeof(xrgb)))
   {
    pixs[npix].row = row;
    pixs[npix].color = COLOR_BLUE;
    npix++;
   }
  }		// for (row = itop; row <= ibottom; row++)

// compile gaps

  ngaps = 0;
  for (pix = 0; pix < npix-1; pix++)
  {
   pix0p = &pixs[pix+0];
   pix1p = &pixs[pix+1];
   if (pix1p->row - pix0p->row > 1)
   {
    gaps[ngaps].pix0p = pix0p;
    gaps[ngaps].pix1p = pix1p;
    ngaps++;
   }
  }

// color gaps

  while (ngaps)
  {
   if (ngaps == 1)					// single gap
   {
    pix0p = gaps[0].pix0p;
    pix1p = gaps[0].pix1p;
    if (pix0p->color == COLOR_RED && pix1p->color == COLOR_BLUE)
     xrgb = BWHITE;
    else
    if (pix0p->color == COLOR_BLUE && pix1p->color == COLOR_RED)
     xrgb = xrgb_mag;
    else
    if (pix0p->color == COLOR_RED && pix1p->color == COLOR_RED)
     xrgb = xrgb_red;
    else
    if (pix0p->color == COLOR_BLUE && pix1p->color == COLOR_BLUE)
     xrgb = xrgb_blue;

    if (memcmp(&xrgb, &BWHITE, sizeof(xrgb)))		// if not white
    {
     ImageColorSet(ip, xrgb);
     ImageMoveTo(ip, pix0p->row+1, col);
     ImageLineTo(ip, pix1p->row-1, col);
    }

    ngaps--;
   }
   else							// more than one gap
   {
    pix0p = gaps[0].pix0p;				// upper gap
    pix1p = gaps[0].pix1p;
    if (pix0p->color == COLOR_RED)			// upper color
     xrgb = xrgb_red;
    else
     xrgb = xrgb_blue;
    ImageColorSet(ip, xrgb);
    ImageMoveTo(ip, pix0p->row+1, col);
    ImageLineTo(ip, pix1p->row-1, col);

    pix0p = gaps[ngaps-1].pix0p;			// lower gap
    pix1p = gaps[ngaps-1].pix1p;
    if (pix1p->color == COLOR_RED)			// lower color
     xrgb = xrgb_red;
    else
     xrgb = xrgb_blue;
    ImageColorSet(ip, xrgb);
    ImageMoveTo(ip, pix0p->row+1, col);
    ImageLineTo(ip, pix1p->row-1, col);

    ngaps -= 2;
    memcpy(&gaps[0], &gaps[1], ngaps * sizeof(GAP));
   }
  }

 }		// for (col = ileft; col <= iright; col++)
/**---------------------------------------------------------------------------*/
// draw grid

 ImageColorSet(ip, LGRAY);

 for (temp = blo; temp <= bhi; temp++)
 {
  row = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);

  for (doy = 0; doy < 366; doy++)
  {
   dayp = &days[doy];
   ImagePixelDraw(ip, row, dayp->left);
  }
 }
/**---------------------------------------------------------------------------*/
// draw month labels

 ImageColorSet(ip, BLACK);

 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];

  ImageMoveTo(ip, bottom,     dayp->left);
  ImageLineTo(ip, bottom + 1, dayp->left);	// tic

  if (dayp->day == 1)
  {
   ImageMoveTo(ip, bottom + 3, dayp->left);
   ImageStringDrawJust(ip, months[dayp->month], JUST_UC);
  }

  if (dayp->day == 10 || dayp->day == 20)
  {
   ImageMoveTo(ip, bottom + 3, dayp->left);
   sprintf(string, "%i", dayp->day);
   ImageStringDrawJust(ip, string, JUST_UC);
  }

// draw dark grid lines

  if (dayp->day == 1 || dayp->day == 10 || dayp->day == 20)
  {
   for (temp = blo; temp <= bhi; temp++)
   {
    row = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
    ImagePixelDraw(ip, row, dayp->left);
   }
  }
 }
/**---------------------------------------------------------------------------*/
// draw location and description

 ImageColorSet(ip, BLACK);

 col = left + ((right - left + 1) - (strlen(location) * fw)) / 2;
 ImageMoveTo(ip, bottom + 3 + fh, col);
 ImageStringDraw(ip, location);

 col = left + ((right - left + 1) - (strlen(description) * fw)) / 2;
 ImageMoveTo(ip, bottom + 3 + fh * 2, col);
 ImageStringDraw(ip, description);
/**---------------------------------------------------------------------------*/
// draw temperature labels

 ImageColorSet(ip, BLACK);

 for (temp = blo; temp <= bhi; temp++)
 {
  row = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);

  ImageMoveTo(ip, row, left);
  ImageLineTo(ip, row, left - 1);		// tic

  itemp = (int)temp;
  if (!(itemp % 5))
  {
   ImageMoveTo(ip, row, left - 3);
   sprintf(string, "%i", itemp);
   ImageStringDrawJust(ip, string, JUST_CR);

// draw dark grid lines

   for (doy = 0; doy < 366; doy++)
   {
    dayp = &days[doy];
    ImagePixelDraw(ip, row, dayp->left);
   }
  }
 }
/**---------------------------------------------------------------------------*/
// draw average high

 ImageColorSet(ip, HICOLOR);

 dayp = &days[365];
 temp = dayp->ahi;
 dayp = &days[0];
 temp += dayp->ahi;
 temp /= 2.0;
 row0 = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
 ImageMoveTo(ip, row0, dayp->left);

 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];
  temp = dayp->ahi;
  row = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
  ImageLineTo(ip, row, dayp->mid);
 }

 ImageLineTo(ip, row0, dayp->right);
/**---------------------------------------------------------------------------*/
// draw average low

 ImageColorSet(ip, LOCOLOR);

 dayp = &days[365];
 temp = dayp->alo;
 dayp = &days[0];
 temp += dayp->alo;
 temp /= 2.0;
 row0 = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
 ImageMoveTo(ip, row0, dayp->left);

 for (doy = 0; doy < 366; doy++)
 {
  dayp = &days[doy];
  temp = dayp->alo;
  row = ibottom - (int)(((temp - blo) / bdelta) * rdelta + 0.5);
  ImageLineTo(ip, row, dayp->mid);
 }

 ImageLineTo(ip, row0, dayp->right);
/**---------------------------------------------------------------------------*/
 if (!l.o)				// if no output file specified
 {
  strcpy(l.filename, pars[1]);		// create one from input file
  dotp = strchr(l.filename, '.');
  if (dotp)
   strcpy(dotp, ".bmp");		// replace filetype
 }

 dotp = strchr(l.filename, '.');
 if (!dotp)
  strcat(l.filename, ".bmp");		// append default filetype

 if (!l.quiet)
  printf("Writing file \"%s\".\n", l.filename);

 ImageBMPSave(ip, l.filename);

 ImageExit(ip);
 ip = 0;
/**---------------------------------------------------------------------------*/
EGRESS:

 if (filep)
  fclose(filep);

 if (arrayp)
  MemFree(arrayp);

 if (mp)
  MatrixFree(mp);

 if (terms)
  MemFree(terms);

 if (ip)
  ImageExit(ip);

 return(rv);
}