``````{
>   Thanks for writing a working Pascal source.  Hopefully it will
>   work with 640x480 resolution (320x200 is a bit grainy, specieally
>   with the default palette.)

I changed Norbert's source a little. Now it looks nicer, and I believe it's
even a fraction faster (not sure, though, didn't time it):
}

{\$G+,N+,E-} { if you have no CoPro, set E+ }

{ Reals   Complex
-1        0
-0.1      0.8
0.3     -0.5
-1.139    0.238
}

program Julia;
const Gseg : word = \$a000;
Type real = double;
var Cx,Cy,Xo,Yo,X1,Y1 : real; Mx,My,A,B,I,Orb : word;

procedure Pset(X,Y : word; C : byte); assembler;
asm
mov es,Gseg
mov ax,[Y]
shl ax,6
mov di,ax
shl ax,2
mov al,[C]
mov [es:di],al
end;

function keypressed : boolean; assembler; asm
mov ah,0bh; int 21h; and al,0feh; end;

procedure Setpalette;
var I : byte;
begin
for I := 1 to 64 do begin
port[\$3c8] := I;
port[\$3c9] := 10+I div 3;
port[\$3c9] := 10+I div 3;
port[\$3c9] := 15+round(I/1.306122449);
end;
end;

begin
write('Real part: '); readln(Cx);
write('Imaginary part: '); readln(Cy);
asm mov ax,13h; int 10h; end;
Setpalette;
Mx := 319; My := 199;
for A := 1 to Mx  do
for B := 1 to My do begin
Xo := -2+A/(Mx/4); { X complex plane coordinate }
Yo :=  2-B/(My/4); { Y complex plane coordinate }
Orb := 0; I := 0;
repeat
X1 := Xo*Xo-Yo*Yo+Cx;
Y1 := 2*Xo*Yo+Cy;
Xo := X1;
Yo := Y1;
inc(I);
until (I = 64) or (X1*X1+Y1*Y1 > 4);
if I <> 64 then Orb := I;
Pset(A,B,Orb); { Plot orbit }
end;
while not keypressed do;
asm mov ax,3; int 10h; end;
end.

``````