#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();
void *OpenLibrary();
int height=256;
struct Window *OpenWindow();
struct Screen *OpenScreen();
struct IntuiMessage *GetMsg();
struct IntuitionBase *IntuitionBase=NULL;
struct GfxBase *GfxBase=NULL;
struct Library *IFFBase=NULL;

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

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

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

short ColorMap[32] = {     /* format 0x0RGB */
 /*  0-7  */  0x0000,0x0F00,0x0FB0,0x0FE4,0x0ED9,0x0DDC,0x0CCD,0x0BBC ,
 /*  8-15 */  0x0AAB,0x099A,0x0889,0x0778,0x0667,0x0556,0x0445,0x0334 ,
 /* 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);
}

double MSetPot(cx,cy,maxiter)
double cx,cy;
long maxiter;
{
	long iter=0,;
	register double x,y,x2,y2,temp,pot=0.0;

	for (x=cx,y=cy,x2=x*x,y2=y*y;(iter<maxiter)&&(x2+y2<10000.0);iter++)
	{
		temp=x2-y2+cx;
		y=2*x*y+cy;
		x=temp;
		x2=x*x;
		y2=y*y;
	}

	if (iter<maxiter) pot=0.5*log(x2+y2)/pow(2.0,(double)iter);
	return(pot);
}

long oldh[640];

MSetCPM3D(xmin,xmax,ymin,ymax,maxiter,angle,scale,lift)
double xmin,xmax,ymin,ymax,scale,lift;
long maxiter,angle;
{
	register long ix,iy,shx,shh,shdep,col,h,maxh,cut,allsteps=0;
	double delta,cx,cy,pot,cutf;

	maxh=w->Height;
	if (angle<0)
	{
		allsteps=1;
		angle=-angle;
	}
	cut=w->Height/2;
	if (ymax-ymin!=xmax-xmin) ymax=ymin+xmax-xmin;
	SetAPen(w->RPort,1L);
	for (iy=0;(iy<angle*w->Height/2)&&(msg==NULL);)
	{
		cy=ymin+iy*(ymax-ymin)/s->Width;
		for (ix=0;(ix<s->Width)&&((msg=GetMsg(w->UserPort))==NULL);ix++)
		{
			cx=xmin+ix*(xmax-xmin)/s->Width;
			pot=(double)MSetPot(cx,cy,maxiter);
			h=maxh*(1-pow(pow(2.0,scale)*pot,lift));
			if (iy>0)
			{
				col=(oldh[ix]+h)/2-oldh[ix-1]+8;
				if (col<4) col=8;
				if (col<6) col=7;
				if ((oldh[ix]>=cut)&&(h>=cut))
				{
					col=8;
					if (oldh[ix]>cut+(maxh-cut)/8) col=7;
					if (oldh[ix]>cut+(maxh-cut)/6) col=6;
					if (oldh[ix]>cut+(maxh-cut)/4) col=5;
					if (oldh[ix]>cut+(maxh-cut)/2) col=4;
					if (oldh[ix]>cut+3*(maxh-cut)/4) col=3;
					if (oldh[ix]>cut+5*(maxh-cut)/6) col=2;
					if (oldh[ix]>cut+7*(maxh-cut)/8) col=1;
					oldh[ix]=cut;
				}
				if (oldh[ix]>=cut) oldh[ix]=cut;
				if (oldh[ix]<cut)
				{
					shh=oldh[ix]-ix;
					for (shdep=0,shx=ix;(shx<w->Width)&&(shx+shh<cut)&&(shdep<2);shx++)
						if (shx+shh<oldh[shx]) shdep++;
					col+=shdep;
				}
				if (col>15) col=15;
				if (pot==0.0) col=0;
				if (ix>0)
				{
					Move(w->RPort,ix,(long)(w->Height/2+(iy-1)/angle));
					SetAPen(w->RPort,col);
					Draw(w->RPort,ix,(long)(w->Height/2+(iy-1)/angle-oldh[ix]));
				}
			}
			oldh[ix]=h;
		}
		if (allsteps==0) iy+=angle;
		else iy++;
	}
	if (msg==NULL)
	{
		for (ix=1;ix<s->Width;ix++)
		{
			SetAPen(w->RPort,9L);
			if (oldh[ix]>cut) oldh[ix]=cut;
			Move(w->RPort,ix,(long)(w->Height/2+(iy-1)/angle));
			Draw(w->RPort,ix,(long)(w->Height/2+(iy-1)/angle-oldh[ix]+1));
		}
	}
}

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("CPM - 3D Mandelbrot program using Continuous Potiential Method");
	puts("      1989 by Lars R. Clausen (2:230/22.34)");
	puts("      Algorithm from 'The Science of Fractal Images'.");
	puts("Usage: CPM [datafile] or");
	puts("       CPM xmin xmax ymin ymin maxit angle scale lift H|L outputfile");
	puts("This program is PD, copy it as you see fit. (No fun for hackers!)");
	if (fromcli==0) gets();
	exit(0);
}

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

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

	if (argc==0) DisplayUsage(0);
	if ((argc==2)&&(strcmp(argv[1],"?")==0)) DisplayUsage(1);
	switch (argc)
	{
		case 2: dataname=argv[1]; break;
		case 11:
			xmin=atof(argv[1]);
			xmax=atof(argv[2]);
			ymin=atof(argv[3]);
			ymax=atof(argv[4]);
			maxit=atoi(argv[5]);
			angle=atoi(argv[6]);
			scale=atof(argv[7]);
			lift=atof(argv[8]);
			size=*argv[9];
			strcpy(savefile,argv[10]);
			fromfile=0;
			size=toupper(size);
			if (size!='H') DoSize(size);
		case 1: break;
		default: DisplayUsage(argc);
	}
	OpenAll();

/* format xmin;xmax;ymin;ymax;maxit;angle;scale;buttom;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);
			getlong(&angle,datafile);
			getdoub(&scale,datafile);
			getdoub(&lift,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],(long)1<<s->BitMap.Depth);
				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");}
		}
		Move(w->RPort,0L,0L);
		ClearScreen(w->RPort);
		printf("Doing: %s Maxit %ld Angle %ld Scale %2.0f Lift %1.5f Size %c\n",savefile,maxit,angle,scale,lift,size);
		MSetCPM3D(xmin,xmax,ymin,ymax,maxit,angle,scale,lift);
		if (msg==NULL) SaveBitMap(savefile,&s->BitMap,&s->ViewPort.ColorMap->ColorTable[0],1L);
		Move(w->RPort,0L,0L);
		ClearScreen(w->RPort);
		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);
}