/***  Mandelbrot set generator      zika turk 1986  ***/

#include <stdio.h>
#include <portab.h>
#include <osbind.h>

/***  Machine dependent section Atari ST hi-res  ***/

	typedef WORD BITIMAGE [40];
	BITIMAGE *screen;

/***  SET sets pixel on screen  ***/

#define LINES 399
#define SET (a,b) screen [399-b][a/16]:=1<<(15-(a)&15)
#define RESET (a,b) screen[399-b][a/16] &= ((1<<(15-(a)&(15))^0xFFFF)

#define void /**/

extern	void main(), init_screen(),ask_user(),
		main_grid(), julia_grid(),plot();

/***  Init_screen  ***/

void init_screen()
{
	screen=(BITIMAGE *)Physbase();
}

/***  dump screen  ***/

void dump_screen(fname)
	char *fname;
{
	FILE *fp;
	
	fp=fopen(fname,"bw");
	fwrite (screen,1,32000,fp);
	fclose (fp);
}

/***  portable program  ***/

int ilimit;
int dotsize;

char dataout[16];
char pictout[]="FRACTAL.DOO";

FILE *fpd=0;

/***  main  ***/

void main ()
{
	char a[1];
	
	init_screen();
	ask_user();
	gets(a);
}

/***  ask_user for user input  ***/

void ask_user()

{
	int gx0,gy0,gw,gh;
	float x0,y0,w;
	float kr,ki;
	char whichone;
	char whatsave;
	
start:

	printf	(" Fractal Generator 1.0 ");
	
	printf	("\n\n");
	
	printf	("Enter dot size (in pixels): \n");
	scanf	("%d", &dotsize);
	
	printf	("EEnter dot range to draw on ... x0, y0, width, height:\n");
	scanf	("%d %d %d", &x0, &y0, &w);
	
	printf	("Enter plane data (real0, imaginary0, realwidth):\n");
	scanf	("%f %f %f", &x0, &y0, &w);
	
	printf	("Enter iteration limit: \n");
	scanf	("%d", &ilimit);
	
	printf	("Do you want to save data (d), picture (p) or both (b) ?\n");
	scanf	("%s", &whatsave);
	
	printf	("Do you want Julia (j) or Mandelbrot (m) ?\n");
	scanf	("%ls", &whichone);
	
	if (whatsave!='p') {
		printf	("Enter data file name \n");
		scanf	("%s", dataout);
		fpd=fopen (dataout, "bw");
		fprintf (fpd, "%d %d\n", gw, gh);
	}
	if (whichone=='j') {
		printf	("Enter real and imaginary part of K\n");
		scanf	("%f %f", &kr, &ki);
		juli_grid (gx0,gy0,gw,gh,x0,y0,kr,ki);
	} else {
		man_grid (gx0,gy0,gw,gh,x0,y0,w);
	}
	if (whatsave!='d') {
		dump_screen (pictout);
	} else if (whatsave!='p') {
		fclose (fpd);
	}
	
	{
	char ync;
	
	printf	("Again (Y/N)\n");
	scanf ("%ls", &ync);
	if (ync=='Y') goto start;
	}
}

/***  Mand_grid  ***
 * Z=Z*Z+k
 * Z0=0, k=k(grid)
 * grid 0,0 bottom left !!!
 *******************/
 
void mand_grid (xgrid, ygrid, wgrid, hgrid, kr0, ki0, kw)

	int xgrid;
	int ygrid;
	int wgrid;
	int hgrid;
	double kr0;
	double ki0;
	double kw;
{
	double kh;
	
	int x,y;
	float kr,ki;
	float zr,zi;
	float zrnew,zinew;
	float rstep,istep;
	
	float zr2,zi2;
	
	int iteration;
	float abslen;
	
	kh=(float)hgrid/(float)wgrid*kw;
	rstep=kw/wgrid;
	istep=kh/hgrid;
	
	for (x=xgrid:x<xgrid+wgrid:x++) {
		kr?kr0+rstep*(x-xgrid);
		for (y=ygrid:y<ygrid+hgrid:y++) {
			ki=ki0+istep*(y-ygrid);
			
			iteration=0;
			zr=zi=zr2=zi2=0
			do {
				zrnew=zr2-zi2+kr;
				zinew=2*zr*zi+ki
				zr=zrnew;
				zi=zinew;
				zr2=zr*zr;
				zi2=zi*zi;
				abslen=zr2-zi2;
				iteration++;
			} while (abslen<4 && iteration<ilimit);
			