/*
 *  (C) 2000 Alan Donovan
 *
 *  Author: Alan Donovan <alan.donovan@arm.com>
 *
 *  newt-term.c -- A Quartic Newton-Raphson-method fractal for terminals!
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 *  See the GNU General Public License for more details.
 *  <http://www.gnu.org/copyleft/gpl.html>.
 *
 *  $Id: newt-term.c,v 1.3 2000/10/27 10:46:13 adonovan Exp adonovan $
 *
 */


#include <stdio.h>
#include <stdlib.h>
#include <math.h>

const double macheps = 1E-10;
const double sin60 = 0.86602496151913422;
const double cos60 = 0.5;

void usage()
{
    fprintf(stderr,
	    "usage: newton [-size <width>x<height>]\n"
	    "              [-scale <num>]\n"
	    "              [-cubic] [-shading]\n");
    exit(1);
}

int main(int argc, char **argv)
{
    double	scale = 2.0;
    int		W=80, H=24; // standard terminal size
    int		ii, sy, sx;
    int 	cubic = 0;
    int 	shading = 0;

    for(ii=1; ii<argc; ++ii)
    {
	if(strcmp(argv[ii], "-size") == 0)
	{
	    if(++ii == argc) 
		usage();
	    if(sscanf(argv[ii], "%ix%i", &W, &H) != 2)
		usage();
	}
	else if(strcmp(argv[ii], "-scale") == 0)
	{
	    if(++ii == argc) 
		usage();
	    if((scale = atof(argv[ii])) == 0.0)
		usage();
	}
	else if(strcmp(argv[ii], "-cubic") == 0)
	{
	    cubic = 1;
	}
	else if(strcmp(argv[ii], "-shading") == 0)
	{
	    shading = 1;
	}
	else if(argv[ii][0] == '-')
	{
	    fprintf(stderr, "newton: unknown option `%s'.\n", argv[ii]);
	    exit(1);
	}
	else
	    usage();
    }

    for(sy = 0; sy < H; ++sy)
    {
	double iy = (scale * (double)(sy - H/2) / (double)H);

	for(sx = 0; sx < W; ++sx)
	{
	    double ix = (scale * (double)(sx - W/2) / (double)W);
	    
	    unsigned root = 0;
	    	
	    double x = ix;
	    double y = iy;

	    for(ii=0; ii<50; ++ii)
	    {
		double x2 = x*x;
		double y2 = y*y;
		double x2_y2 = x2 - y2;

		double a,b,c,d,e,f,g;

		if(cubic)
		{
		    a = x * (x2 - 3*y2) - 1;
		    b = y * (3*x2 - y2);
		    c = 3*x2_y2;
		    d = 6*x*y;
		}
		else // quartic
		{
		    a = x2*x2 - 6*x2*y2 + y2*y2 - 1;
		    b = 4*x*y * x2_y2;
		    c = 4*x * (x2_y2 - 2*y2);
		    d = 4*y * (x2_y2 + 2*x2);
		}

		e = c*c + d*d;
		e = (e == 0) ? 1E10 : (1.0 / e);
		
		f = a*c + b*d;
		g = b*c - a*d;

		f *= e;
		g *= e;

		x -= f;
		y -= g;

		// check roots
		
		if(cubic)
		{
		    if((fabs(x - 1) < macheps) && // 1
		       (fabs(y + 0) < macheps))
		    {
			root = 1;
		    }
		    else if((fabs(x + cos60) < macheps) && // -cos60, sin60
			    (fabs(y - sin60) < macheps))
		    {
			root = 2;
		    }
		    else if((fabs(x + cos60) < macheps) && // -cos60, -sin60
			    (fabs(y + sin60) < macheps))
		    {
			root = 3;
		    }	
		} 
		else // quartic
		{
		    if((fabs(x - 1) < macheps) && // 1
		       (fabs(y + 0) < macheps))
		    {
			root = 1;
		    }
		    else if((fabs(x + 1) < macheps) && // -1
			    (fabs(y + 0) < macheps))
		    {
			root = 2;
		    }
		    else if((fabs(x + 0) < macheps) && // i
			    (fabs(y - 1) < macheps))
		    {
			root = 3;
		    }
		    else if((fabs(x + 0) < macheps) && // -i
			    (fabs(y + 1) < macheps))
		    {
			root = 4;
		    }
		}

		if(root) break;
	    } 
	    fputc(root[" #$.,"], stdout);

	    if(shading)
		fputc((ii%16)["ABCDEFGHIJKLMNOP"], stdout);
	}
    }

    return 0;
}