/* genplot - use "-l" output from mkfilter to generate .png mag/phase plots
   A.J. Fisher, University of York   <fisher@minster.york.ac.uk>
   April 1995 

   16 July 2000 by Nikolai Golovchenko
		1) Replaced GIF format by PNG generation (http://www.boutell.com/gd)
		2) Compiled for Win32
		3) Corrected set_new_handler evocation
*/

/* Uses the "gd" gif manipulation library, originally from
   Quest Protein Database Center, and now available from:
   http://www.boutell.com/gd
*/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <new.h>

extern "C"
{
   #include <gd.h>
   #include <gdfonts.h>
};

#include "mkfilter.h"
#include "complex.h"
#include "genplot.h"

char fmt[10];

global void main(int argc, char *argv[])
{
	char junkv[MAXSTRING+1];
	double junkd;

	_set_new_handler(newhandler);
	readcmdline(argv);
	readdata(junkv, junkd, nzeros, xcoeffs, npoles, ycoeffs);
	drawgraph();
	exit(0);
}

void readcmdline(char *argv[])
{
	int ap = 0;
	char *s;
	outfn = NULL; gtype = unknown; logopt = dflag = false;
	alpha1 = 0.0; alpha2 = 0.5;
	progname = (argv[ap] != NULL) ? argv[ap++] : "???";
	until (argv[ap] == NULL)
	{
		s = argv[ap++];
		if (seq(s,"-a"))
		{
			alpha1 = getfarg(argv[ap++]);
			alpha2 = getfarg(argv[ap++]);
			unless (alpha2 > alpha1) giveup("Invalid values for -a: a1 must be .lt. a2");
		}
		else if (gtype == unknown && seq(s,"-l"))
		{
			gtype = locus;
		}
		else if (gtype == unknown && (seq(s,"-i") || seq(s,"-s")))
		{
			gtype = (s[1] == 's') ? step : impulse;
			nirsteps = getiarg(argv[ap++]);
			if (nirsteps < 0 || nirsteps > MAXSTEPS) giveup("Invalid value for num steps: limits are 0 .. %d", MAXSTEPS);
		}
		else if (seq(s,"-log"))
		{
			logmin = getfarg(argv[ap++]);
			if (logmin >= 0.0) giveup("Invalid value for log min (%s): must be .lt. 0", argv[ap-1]);
			logopt = true;
		}
		else if (seq(s,"-d")) 
		{
			dflag = true;
		}
		else
		{
			if (outfn != NULL) 
			{
				usage();
			}
			outfn = s;
		}
	}
	if (outfn == NULL)
	{
		usage();
	}
	if (gtype == unknown) 
	{
		gtype = phmag;    /* default */
	}
}

double getfarg(char *s)
{
	if (s == NULL)
	{
		usage();
	}
   return atof(s);
}

int getiarg(char *s)
  { if (s == NULL) usage();
    return atoi(s);
  }

void usage()
  { fprintf(stderr, "Genplot V.%s from <fisher@minster.york.ac.uk>\n", VERSION);
    fprintf(stderr, "Usage: genplot [-i <n> | -s <n> -a <a1> <a2> | -l] [-log <min>] [-d] out.png\n");
    exit(1);
  }

void drawgraph()
  { switch (gtype)
      { default:
	    giveup("bug (gtype)");

	case locus:
	  { int nsteps = MINSTEPS;
	    bool ok = false; double stepl;
	    until (ok || nsteps > MAXSTEPS)
	      { complex *fr = new complex[nsteps+1];
		computefr(fr, nsteps);
		locusgraph(outfn, fr, nsteps, stepl);
		if (stepl < 0.1) ok = true;
		nsteps *= 10;
		delete fr;
	      }
	    unless (ok) fprintf(stderr, "genplot: warning: step length %g in locus plot\n", stepl);
	    break;
	  }

	case phmag:
	  { complex *fr = new complex[FRSTEPS+1];
	    computefr(fr, FRSTEPS);
	    frgraph(outfn, fr, FRSTEPS);
	    delete fr;
	    break;
	  }

	case impulse:
	  { double *ir = new double[nirsteps+1];
	    computeir(ir, nirsteps, false);
	    irgraph(outfn, ir, nirsteps);
	    delete ir;
	    break;
	  }

	case step:
	  { double *ir = new double[nirsteps+1];
	    computeir(ir, nirsteps, true);
	    irgraph(outfn, ir, nirsteps);
	    delete ir;
	    break;
	  }
      }
  }

void computefr(complex fr[], int nsteps)
  { complex topcos[MAXPZ+1], botcos[MAXPZ+1]; int i;
    for (i=0; i <= nzeros; i++) topcos[i] = complex(xcoeffs[i]);
    for (i=0; i <= npoles; i++) botcos[i] = complex(-ycoeffs[i]);
    double maxmag = 0.0;
    for (i = 0; i <= nsteps; i++)
      { double a = (double) i / (double) nsteps;
	double theta = TWOPI * (alpha1 + a * (alpha2-alpha1));
	complex z = expj(theta);
	fr[i] = evaluate(topcos, nzeros, botcos, npoles, z);
	double mag = hypot(fr[i]);
	if (mag > maxmag) maxmag = mag;
	if (dflag) fr[i] = fr[i] * expj(0.5 * nzeros * theta);
      }
    if (maxmag > 0.0)
      { for (i=0; i <= nsteps; i++) fr[i] /= maxmag;  /* scale so max mag is 1.0 */
      }
  }

inline void shiftdown(double v[], int n)
  { memcpy(&v[0], &v[1], n * sizeof(double));
  }

void computeir(double ir[], int nsteps, bool stp)
  { double xv[MAXPZ+1], yv[MAXPZ+1];
    double maxmag = 0.0; int i, j;
    for (i=0; i <= nzeros; i++) xv[i] = 0.0;
    for (i=0; i <= npoles; i++) yv[i] = 0.0;
    for (i=0; i <= nsteps; i++)
      { shiftdown(xv, nzeros);
	xv[nzeros] = (i == 0 || stp) ? 1.0 : 0.0;	/* impulse or step */
	shiftdown(yv, npoles);
	yv[npoles] = 0.0;
	for (j=0; j <= nzeros; j++) yv[npoles] += (xcoeffs[j] * xv[j]);
	for (j=0; j < npoles; j++) yv[npoles] += (ycoeffs[j] * yv[j]);
	ir[i] = yv[npoles];
	double mag = fabs(ir[i]);
	if (mag > maxmag) maxmag = mag;
      }
    if (maxmag > 0.0)
      { for (i=0; i <= nsteps; i++) ir[i] /= maxmag;  /* scale to range +-1 */
      }
  }

void locusgraph(char *fn, complex fr[], int nsteps, double &stepl)
{
	stepl = 0.0;
	gdImagePtr im = gdImageCreate(XYPIXELS + 2*BORDER, XYPIXELS + 2*BORDER);
	gdImageColorAllocate(im, 255, 255, 255); /* white background */
	int red = gdImageColorAllocate(im, 255, 0, 0);
	int black = gdImageColorAllocate(im, 0, 0, 0);
	gdImageLine(im, xmap_lc(-1.0), ymap_lc(0.0), xmap_lc(+1.0), ymap_lc(0.0), black); /* X axis	 */
	gdImageLine(im, xmap_lc(0.0), ymap_lc(-1.0), xmap_lc(0.0), ymap_lc(+1.0), black); /* Y axis	 */
	for (int i = 1; i <= nsteps; i++)
	{ complex z0 = fr[i-1], z1 = fr[i];
	gdImageLine(im, xmap_lc(z0.re), ymap_lc(z0.im), xmap_lc(z1.re), ymap_lc(z1.im), red);
	double sl = hypot(z1 - z0);
	if (sl > stepl) stepl = sl;
	}
	FILE *fi = fopen(fn, "wb");
	if (fi == NULL) giveup("can't create %s", fn);
	//NG//gdImageGif(im, fi);
	gdImagePng(im, fi);
	fclose(fi);
	gdImageDestroy(im);
}

int xmap_lc(double x)
  { return (int) (XYPIXELS * 0.5 * (x + 1.0)) + BORDER;
  }

int ymap_lc(double y)
  { return (int) (XYPIXELS * 0.5 * (1.0 - y)) + BORDER;
  }

void frgraph(char *fn, complex fr[], int nsteps)
  { double *ymvals = new double[nsteps+1], *ypvals = new double[nsteps+1]; int i; char *fmt;
    gdImagePtr im = gdImageCreate(XYPIXELS + 2*BORDER, XYPIXELS + 2*BORDER);
    gdImageColorAllocate(im, 255, 255, 255); /* white background */
    int red = gdImageColorAllocate(im, 255, 0, 0);
    int blue = gdImageColorAllocate(im, 0, 0, 255);
    int black = gdImageColorAllocate(im, 0, 0, 0);
    if (dflag)
      { char str[MAXSTRING+1];
	sprintf(str, "FILTER DELAY: %g SAMPLES.", 0.5 * nzeros);
	gdImageString(im, gdFontSmall, 20, 10, (unsigned char*)str, black);
      }
    gdImageLine(im, xmap_pm(0.0), ymap_pm(0.0), xmap_pm(1.0), ymap_pm(0.0), black); /* X axis	    */
    gdImageLine(im, xmap_pm(0.0), ymap_pm(0.0), xmap_pm(0.0), ymap_pm(1.0), black); /* left Y axis  */
    gdImageLine(im, xmap_pm(1.0), ymap_pm(0.0), xmap_pm(1.0), ymap_pm(1.0), black); /* right Y axis */
    for (i=0; i <= nsteps; i++)
      { ymvals[i] = hypot(fr[i]);
	ypvals[i] = (ymvals[i] > PHEPS) ? atan2(fr[i]) / TWOPI + 0.5 :
		    -1.0; /* out-of-range value will not be plotted */
      }
    /* ymvals, ypvals in range 0 .. 1 */
    if (logopt)
      { for (i=0; i <= nsteps; i++)
	  ymvals[i] = (ymvals[i] > 0.0) ? (20.0 * log10(ymvals[i]) - logmin) / (-logmin) :
					  -1.0; /* out-of-range value will not be plotted */
      }
    pgraph(im, ymvals, nsteps, red, false);    /* mag graph   */
    pgraph(im, ypvals, nsteps, blue, true);    /* phase graph */
    if (logopt)
      { fmt = choosefmt(logmin, 0.0);
	for (i = 0; i <= 10; i++) draw_ytick(im, 0.0, 0.1*i, red, -TICKLEN, fmt, (1.0 - 0.1*i) * logmin);
      }
    else
      { for (i = 0; i <= 10; i++) draw_ytick(im, 0.0, 0.1*i, red, -TICKLEN, "%3.1f", 0.1*i);
      }
    draw_ytick(im, 1.0, 0.0, blue, TICKLEN, "-pi");
    draw_ytick(im, 1.0, 0.5, blue, TICKLEN, "0");
    draw_ytick(im, 1.0, 1.0, blue, TICKLEN, "+pi");
    for (i = 1; i <= 9; i++) draw_ytick(im, 1.0, 0.1*i, blue, TICKLEN, "", 0.2*(i-5));
    fmt = choosefmt(alpha1, alpha2);
    for (i = 0; i <= 10; i++)
      { double a = (double) i / 10.0;
	double x = alpha1 + a * (alpha2-alpha1);
	draw_xtick(im, 0.1*i, 0.0, black, TICKLEN, fmt, x);
      }
    FILE *fi = fopen(fn, "wb");
    if (fi == NULL) giveup("can't create %s", fn);
    //NG//gdImageGif(im, fi);
	 gdImagePng(im, fi);
    fclose(fi);
    gdImageDestroy(im);
    delete ymvals; delete ypvals;
  }

void irgraph(char *fn, double ir[], int nsteps)
{
	double *ymvals = new double[nsteps+1]; int i;
	gdImagePtr im = gdImageCreate(XYPIXELS + 2*BORDER, XYPIXELS + 2*BORDER);
	gdImageColorAllocate(im, 255, 255, 255); /* white background */
	int red = gdImageColorAllocate(im, 255, 0, 0);
	int black = gdImageColorAllocate(im, 0, 0, 0);
	gdImageLine(im, xmap_pm(0.0), ymap_pm(0.0), xmap_pm(1.0), ymap_pm(0.0), black); /* X axis	    */
	gdImageLine(im, xmap_pm(0.0), ymap_pm(0.0), xmap_pm(0.0), ymap_pm(1.0), black); /* left Y axis  */
	gdImageLine(im, xmap_pm(1.0), ymap_pm(0.0), xmap_pm(1.0), ymap_pm(1.0), black); /* right Y axis */
	for (i=0; i <= nsteps; i++) ymvals[i] = 0.5 * (ir[i] + 1.0);
	/* ymvals in range 0 .. 1 */
	pgraph(im, ymvals, nsteps, red, false);
	for (i=0; i <= 10; i++) draw_ytick(im, 0.0, 0.1*i, red, -TICKLEN, "%3.1f", 0.2*(i-5));
	for (i=0; i <= 10; i++)
	{ double x = (double) nsteps * (double) i / 10.0;
	draw_xtick(im, 0.1*i, 0.0, black, TICKLEN, "%g", x);
	}
	FILE *fi = fopen(fn, "wb");
	if (fi == NULL) giveup("can't create %s", fn);
	//NG//gdImageGif(im, fi);
	gdImagePng(im, fi);
	fclose(fi);
	gdImageDestroy(im);
	delete ymvals;
}

void pgraph(gdImagePtr im, double vec[], int nsteps, int col, bool wrapy)
  { for (int i = 1; i <= nsteps; i++)
      { double x0 = (double) (i-1) / (double) nsteps,
	       x1 = (double) i / (double) nsteps;
	double y0 = vec[i-1], y1 = vec[i];
	if (y0 >= 0.0 && y1 >= 0.0)	/* suppress discontinuities */
	  { if (wrapy && y0-y1 > 0.5)
	      { gdImageLine(im, xmap_pm(x0), ymap_pm(y0), xmap_pm(x1), ymap_pm(1.0), col);
		gdImageLine(im, xmap_pm(x0), ymap_pm(0.0), xmap_pm(x1), ymap_pm(y1), col);
	      }
	    else if (wrapy && y1-y0 > 0.5)
	      { gdImageLine(im, xmap_pm(x0), ymap_pm(y0), xmap_pm(x1), ymap_pm(0.0), col);
		gdImageLine(im, xmap_pm(x0), ymap_pm(1.0), xmap_pm(x1), ymap_pm(y1), col);
	      }
	    else gdImageLine(im, xmap_pm(x0), ymap_pm(y0), xmap_pm(x1), ymap_pm(y1), col);
	  }
      }
  }

void draw_xtick(gdImagePtr im, double x, double y, int col, int tlen, char *fmt, double val)
  { int ix = xmap_pm(x), iy = ymap_pm(y);
    gdImageLine(im, ix, iy, ix, iy+tlen, col);
    char str[32]; sprintf(str, fmt, val);
    int swid = strlen(str) * (gdFontSmall -> w);
    gdImageString(im, gdFontSmall, ix - swid/2, iy + 2*tlen, (unsigned char*)str, col);
  }

void draw_ytick(gdImagePtr im, double x, double y, int col, int tlen, char *fmt, double val)
  { int ix = xmap_pm(x), iy = ymap_pm(y);
    gdImageLine(im, ix, iy, ix+tlen, iy, col);
    char str[32]; sprintf(str, fmt, val);
    int swid = strlen(str) * (gdFontSmall -> w);
    int shgt = gdFontSmall -> h;
    if (tlen < 0) gdImageString(im, gdFontSmall, ix + 2*tlen - swid, iy - shgt/2, (unsigned char*)str, col);
    if (tlen > 0) gdImageString(im, gdFontSmall, ix + 2*tlen, iy - shgt/2, (unsigned char*)str, col);
  }

int xmap_pm(double x)
  { return (int) (XYPIXELS * x) + BORDER;
  }

int ymap_pm(double y)
  { return (int) (XYPIXELS * (1.0 - y)) + BORDER;
  }


char *choosefmt(double a1, double a2)
  { /* choose smallest format such that all ticks are labelled distinctly */
     bool ok;
    int p = 0;
    do
      { sprintf(fmt, "%%.%df", p++);
	ok = true;
	char buf1[20], buf2[20]; buf2[0] = '\0';
	for (int i = 0; i <= 10 && ok; i++)
	  { double a = (double) i / 10.0;
	    double x = a1 + a * (a2-a1);
	    strcpy(buf1, buf2); sprintf(buf2, fmt, x);
	    if (seq(buf1, buf2)) ok = false;
	  }
      }
    until (ok || p > 10);
    unless (ok) strcpy(fmt, "?.?");
    return fmt;
  }

int
newhandler(size_t size)
{
   printf( "Allocation failed, " );
   printf( "%u bytes not available.\n", size );
	giveup("No room");
   // Tell new to stop allocation attempts.
   return 0;
}

void giveup(char *msg, word p1)
  { fprintf(stderr, "genplot: ");
    fprintf(stderr, msg, p1.i); putc('\n', stderr);
    exit(1);
  }