#include <exec/types.h>
#include <graphics/gfxbase.h>
#include <intuition/intuition.h>
#include <math.h>
#include <libraries/iff.h>
#include <ctype.h>
#include <stdio.h>

#define BITS2SHIFT 27

double atof();
long atol();
int height=256;
void *OpenLibrary(),*GetBMHD(),*OpenIFF(),*AllocRaster();
struct Window *OpenWindow();
struct Screen *OpenScreen();
struct IntuiMessage *GetMsg();
struct IntuitionBase *IntuitionBase;
struct GfxBase *GfxBase;
struct Library *IFFBase=NULL;

struct Window *w;
struct Screen *s;
struct IntuiMessage *msg;
struct BitMap *SBitMap=NULL;

struct NewScreen ns = {0,0,640,512,1,0,0,HIRES|LACE,CUSTOMSCREEN,NULL,NULL,NULL,NULL};

struct NewWindow nw = {
	0,1,640,511,-1,-1,VANILLAKEY|CLOSEWINDOW,WINDOWCLOSE|NOCAREREFRESH|SIMPLE_REFRESH|
	ACTIVATE,NULL,NULL,NULL,NULL,NULL,0,0,0,0,CUSTOMSCREEN};

short ColorMap[32] = {     /* format 0x0RGB */
 /*  0-7  */  0x0000,0x0FFF,0x0EEE,0x0DDD,0x0CCC,0x0BBB,0x0AAA,0x0999 ,
 /*  8-15 */  0x0888,0x0777,0x0666,0x0555,0x0444,0x0333,0x0222,0x0111 ,
 /* 16-23 */  0x066F,0x0FFF,0x0EEE,0x0DDD,0x0CCC,0x0BBB,0x0AAA,0x0999 ,
 /* 24-31 */  0x0888,0x0777,0x0666,0x0555,0x0444,0x0333,0x0222,0x0111 };

OpenAll()
{
	if (!(IntuitionBase = (struct IntuitionBase *) OpenLibrary("intuition.library", 0L)))
		abort("No Intuition");
	if (!(GfxBase = (struct GfxBase *) OpenLibrary("graphics.library",0L)))
		abort("No GFXBase");
	ns.Height=height=GfxBase->NormalDisplayRows;
	nw.Height=height-1;
	
	if (!(IFFBase = (struct Library *) OpenLibrary(IFFNAME,IFFVERSION)))
		abort("Somebody's using IFF-library - or it is absent");
	if ((s=(struct Screen *)OpenScreen(&ns))==NULL)
		abort("No screen here");

	LoadRGB4(&s->ViewPort,&ColorMap[0],(long)1<<s->BitMap.Depth);
	nw.Screen=s;
	if ((w = (struct Window *) OpenWindow(&nw)) == NULL)
		abort("Can't (and won't) open window");
}

abort(exitcode)
char *exitcode;
{
	if (exitcode) puts(exitcode);
	if (w) CloseWindow(w);
	if (s) CloseScreen(s);
	if (IFFBase) CloseLibrary(IFFBase);
	if (GfxBase) CloseLibrary(GfxBase);
		OpenWorkBench();
	if (IntuitionBase) CloseLibrary(IntuitionBase);
	exit(0);
}

max(x,y)
double x,y;
{
	if (x>y) return(x); else return(y);
}

long LoadPicture(filename)
char *filename;
{
	APTR ifffile = NULL;
	struct BitMapHeader *bmhd;
	long count,bd,bw,bh;
	register long x,y,inity;
	UWORD colortable[128];
	char *iffinfile=":Pics/WBBackPic.ILBM";
	
	if(!(ifffile=OpenIFF(filename))) {
		return(0L);
	}
	if(!(bmhd=GetBMHD(ifffile)))
	{
		abort("Bad BitMap Header on IFF file");
	}
	count = GetColorTab(ifffile,colortable);
/* 	LoadRGB4(&(s->ViewPort),colortable,count); */

	if(!(DecodePic(ifffile,&s->BitMap)))
	{
		abort("Decoding IFF - bad things happened!");
	}
	CloseIFF(ifffile);
	printf("Analyzing picture for resume...");
	fflush(stdout);
	for (y=0,x=0;(y<w->Height)&&(x<w->Height);y++)
		for (x=0;(x<w->Width)&&(ReadPixel(w->RPort,x,y)==0);x++);
	if (y>1) y-=2;
	else y=0;
	printf(" resuming at line %ld.\n",y);
	DisplayBeep(0L);
	return(y);
}

double xorbit[2001],yorbit[2001];

double MSetDist(cx,cy,maxiter)
double cx,cy;
long maxiter;
{
	long iter=0,i,flag;
	double huge=100000.0;
	double x=0.0,y=0.0,x2=0.0,y2=0.0,temp,xder,yder;
	double dist=0.0;

	xorbit[0]=yorbit[0]=0.0;
	do
	{
		temp=x2-y2+cx;
		y=2*x*y+cy;
		x=temp;
		x2=x*x;
		y2=y*y;
		iter++;
		xorbit[iter]=x;
		yorbit[iter]=y;
	}
	while ((iter<maxiter)&&(x2+y2<huge));

	if (x2+y2>huge)
	{
		xder=yder=0.0;
		i=flag=0;
		do
		{
			temp=2*(xorbit[i]*xder-yorbit[i]*yder)+1;
			yder=2*(yorbit[i]*xder+xorbit[i]*yder);
			xder=temp;
			if ((max(abs(xder),abs(yder)))>10000000.0) flag=1;
			i++;
		}
		while ((i<iter)&&(!flag));
		if ((!flag)&&((xder!=0)||(yder!=0)))
			dist=log(x2+y2)*sqrt(x2+y2)/sqrt(xder*xder+yder*yder);
	}
	return(dist);
}

MSetDEM(inity,xmin,xmax,ymin,ymax,maxiter,treshold,savefile)
long inity;
double xmin,xmax,ymin,ymax;
long maxiter;
double treshold;
char *savefile;
{
	long ix,iy;
	double delta,cx,cy,dist;

	delta=treshold*(xmax-xmin)/s->Width;
	ymax=ymin+(w->Width/w->Height)*(xmax-xmin);
	SetAPen(w->RPort,1L);
	msg=GetMsg(w->UserPort);
	for (iy=inity;(iy<w->Height)&&(msg->Code!=27)&&(msg->Class!=CLOSEWINDOW);iy++)
	{
		cy=ymin+iy*(ymax-ymin)/s->Width;
		for (ix=0;(ix<s->Width)&&(msg->Code!=27)&&(msg->Class!=CLOSEWINDOW);ix++)
		{
			cx=xmin+ix*(xmax-xmin)/s->Width;
			dist=(double)MSetDist(cx,cy,maxiter);
			if (dist>=delta) WritePixel(w->RPort,ix,iy);
			msg=GetMsg(w->UserPort);
			if (msg->Code==115) SaveBitMap(savefile,&s->BitMap,&s->ViewPort.ColorMap->ColorTable[0],1L);
		}
	}
}


getdoub(doub,datafile)
double *doub;
FILE *datafile;
{
	int i;
	char doublebuf[30];

	i=-1;
	doublebuf[0]=0;
	do fread(&doublebuf[++i],1,1,datafile);
	while ((doublebuf[i]!=';')&&(i<30));
	doublebuf[i]=0;
	*doub=atof(&doublebuf[0]);
}

getlong(lon,datafile)
long *lon;
FILE *datafile;
{
	int i;
	char doublebuf[30];

	i=-1;
	doublebuf[0]=0;
	do fread(&doublebuf[++i],1,1,datafile);
	while ((doublebuf[i]!=';')&&(i<30));
	*lon=atol(&doublebuf[0]);
}

DisplayUsage(FromCLI)
int FromCLI;
{
	puts("DEM - B & W MandelBrot generator using Distance Estimator Method");
	puts("      Algorithm from 'The Science of Fractal Images.");
	puts("(c) 1989 by Lars R. Clausen (2:230/22.34)");
	puts("Usage: DEM [datafilename] or");
	puts("       DEM xmin xmax ymin ymax maxit [L|W|H|G|F] savefilename");
	puts("This program is PD (no fun for hackers).");
	puts("Improvements, bugreports and $$$ are welcome.");
	if (!FromCLI) gets();
	exit(0);
}

DoSize(size)
char size;
{
	switch (size)
	{
	case 'L':	ns.Height=height; ns.Width=320;
				nw.Height=height-1; nw.Width=320;
				ns.ViewModes=0;
				break;
	case 'W':	ns.Height=height; ns.Width=640;
				nw.Height=height-1; nw.Width=640;
				ns.ViewModes=HIRES;
				break;
	case 'H':	ns.Height=2*height; ns.Width=640;
				nw.Height=2*height-1; nw.Width=640;
				ns.ViewModes=HIRES|LACE;
				break;
	case 'G':   ns.Height=1024; ns.Width=960;
				nw.Height=1023; nw.Width=960;
				ns.ViewModes=HIRES|LACE;
				break;
	case 'F':	ns.Height=ns.Width=nw.Width=1024;
				nw.Height=1023;
				ns.ViewModes=HIRES|LACE;
				break;
	default: 	puts("Illegal screen type");
	}
}

main(argc,argv)
int argc;
char *argv[];
{
	long i,end=0,fromfile=1,maxit,inity;
	FILE *datafile=NULL;
	double xmin,xmax,ymin,ymax;
	char size,oldsize='H',savefile[200],*dataname="DEMDataFile";

	if ((argc==2)&&(strcmp(argv[1],"?")==0)) DisplayUsage(1);
	switch (argc)
	{
		case 1: break;
		case 2: dataname=argv[1]; break;
		case 8:
			xmin=atof(argv[1]);
			xmax=atof(argv[2]);
			ymin=atof(argv[3]);
			ymax=atof(argv[4]);
			maxit=atoi(argv[5]);
			size=*argv[6];
			strcpy(savefile,argv[7]);
			fromfile=0;
			size=toupper(size);
			if (size!='H') DoSize(size);
			oldsize=size;
			break;
		default: DisplayUsage(argc);
	}
	OpenAll();
	if (GfxBase->DisplayFlags&NTSC) height=200;

/* format xmin;xmax;ymin;ymax;maxit;size;savefile[EOL] */
	if (fromfile)
	{
		datafile=fopen(dataname,"r");
		if (datafile==NULL) abort("Can't open file");
	}
	do
	{
		if (fromfile)
		{
			getdoub(&xmin,datafile);
			getdoub(&xmax,datafile);
			getdoub(&ymin,datafile);
			getdoub(&ymax,datafile);
			getlong(&maxit,datafile);
			do
				fread(&size,1,1,datafile);
			while (!(isalpha(size)));
			size=toupper(size);
			if (oldsize!=size)
			{
				CloseWindow(w);
				CloseScreen(s);
				DoSize(size);
				if ((s=(struct Screen *)OpenScreen(&ns))==NULL)
					abort("No screen here");

				LoadRGB4(&s->ViewPort,&ColorMap[0],1L);
				nw.Screen=s;
				if ((w = (struct Window *) OpenWindow(&nw)) == NULL)
					abort("Can't (and won't) open window");
			}
			i=-1;
			fread(&savefile[0],1,1,datafile);
			do
				fread(&savefile[++i],1,1,datafile);
			while ((savefile[i]!=EOF)&&(savefile[i]!='\n')&&(i<200));
			if (savefile[i]==EOF) end=1;
			savefile[i]=0;
			if (i==200) {fclose(datafile);abort("filename too long");}
		}
		oldsize=size;
		Move(w->RPort,0L,0L);
		ClearScreen(w->RPort);
		inity=LoadPicture(savefile);
		printf("Doing: %s Maxit %ld Size %c\n",savefile,maxit,size);
		MSetDEM(inity,xmin,xmax,ymin,ymax,maxit,0.5,savefile);
		if (msg==NULL) SaveBitMap(savefile,&s->BitMap,&s->ViewPort.ColorMap->ColorTable[0],1L);
		if (msg->Class==CLOSEWINDOW) end=1;
		for (msg=GetMsg(w->UserPort);msg!=NULL;ReplyMsg(msg),msg=GetMsg(w->UserPort));
		msg=NULL;
	}
	while ((end==0)&&(fromfile==1));
	if (fromfile) fclose(datafile);
	abort(NULL);
}