Fraktalai.

Paskalio programa Mandelbroto aibei gauti

Programa, kurios dėka galite gauti Mandelbroto aibės vaizdą, turėtų atrodyti apytiksliai taip:

      FOR EACH pixel:
         REPEAT
           CALCULATE Z(n)=Z(n-1)^2+C
           iterations=iterations+1
         UNTIL value=infinite OR iterations>max_iterations
         SET COLOR=ITERATIONS

(T.y., su kiekvienu ekrano tašku atliekame procesą Z(n)=(Z(n-1))²+C, kurį aptarėme puslapyje Mandelbroto aibė, iki tol, kol nustatysime, kad Z(n) tolsta į begalybę arba kol procesą pakartosime tam tikrą skaičių, kurį pažymėjome max_iterations. Tada taško spalva bus tokia, kurios numeris atitinka proceso pakartojimų (iteracijų) skaičių).
Kaip minėta, Z(n)=(Z(n-1))²+C ir .Tai reiškia:

Z(1)=C
       Z(2)=C²+C
                            Z(3)=(C²+C)²+C    ir t.t.     Kaip žinia C=a+ib. Tuomet skaičiuosime pagal formules:
            Re Z(n)=(Re Z(n-1)) ²-(Im Z(n-1)) ²+Re C
            Im Z(n)=2* Re Z(n-1)* Im Z(n-1)+Im C
    (Tai iš formulių (x+iy)²+a+ib=(x²-y²+a)+(2xy+b)i, kur Z(n-1)=x+iy).
    Dabar turime nustatyti, kada Z(n) tolsta į begalybę.
    Z(n) reiškia vektorių, išvestą iš koordinačių pradžios ir turintį koordinates {Re Z(n), Im Z(n)}. Tuomet, mes galime rasti jo ilgį (ILGIS=kvadratinė šaknis iš (Re Z(n)) ²+(Im Z(n)) ²).
    Taigi, kada ilgis begalinis? Bloga žinia: Jūs negalite patikrinti, ar jis begalinis. Gera žinia: jeigu tik kada ilgis viršys 2, tai jis tikrai taps (kurią tai dieną) begaliniu.
    Patikslinsime programą:

        FOR EACH pixel:
             REPEAT
                CALCULATE Z(n)=Z(n-1)^2+C
                Length=(Re Z(n)) ^2 +(Im Z(n)) ^2
                iterations=iterations+1
             UNTIL Length>4 OR iterations>max_iterations
        SET COLOR=ITERATIONS

    Pastebėkite, kad skaičiavome, ar ilgis viršija 4, o ne 2, ir kad skaičiuodami ilgį netraukėme šaknies.
    Turėtumėte žinoti, kad šaknies traukimas užima daug laiko, tad kodėl gi negalime skaičiuoti 'Length'>4, kai SQRT(Length) >2 ?!?
    Kad būtų patogiau, parašykime funkciją, kur C yra jos kintamasis, o iteracijų skaičius yra jos reikšmė.

    FUNCTION CALC_PIXEL(CA,CBi: REAL): INTEGER;
    CONST MAX_ITERATION=128;
    VAR
        OLD_A:REAL                     {pagalbinis kintamasis 'a' reikšmei saugoti}
        ITERATION:INTEGER      {skaičiuoja iteracijas}
        A,B:REAL                           {z realioji ir menamoji dalys}
        LENGTH_Z:REAL              {z ilgis kvadratu}
    BEGIN
        A:=0; B:=0;                         {Z(0):=0
        ITERATION:=0;
        REPEAT
            OLD_A:=A;
            A:=A*A-B*B+CA;
            B:=2*OLD_A*B+CBi;
            length_z:=a*a+b*b;          {netraukėme šaknies !}
        UNTIL (length_z>4) OR (iteration>max_iteration);
        Calc_Pixel:=iteration;
    END;

Mandelbroto aibės radimui jums reikės tokios programos:

    FOR y:=-1.25 to 1.25 DO
          FOR x:=2 to 1.25 DO
         color:=CALC_PIXEL(x,y);

Bet jūs dar negalite aibės pavaizduoti ekrane. Prisiminę apie ekrano koordinates (mūsų atveju naudojamės 640*480VGA):

  FOR y:=0 to 480-1 DO
   FOR x:=0 to 640-1 DO
    BEGIN
       color:=CALC_PIXEL(Re(x),Im(y));
       PUTPIXEL(X,Y, color);
    END;

Ši programa jau gali pavaizduoti Mandelbroto aibę ekrane, bet dar reikia apibrėžti, kaip skaičiuoti realią ir menamą kiekvieno taško dalį.
Kad gauti visą Mandelbroto aibę, x turi kisti nuo -2 iki 1.25, o y nuo -1.25 iki1.25. Šiuos rėžius aprašysime konstantų dalyje:

PROGRAM Mandelbrot;
CONST MinX=-2;
    MaxX=1.25;
    MinX=-1.25;
    MaxX=1.25;
   [likusi programos dalis...]

    Jei viršutinis kairysis kampas (MinX, MinY), o dešinysis apatinis (MaxX, MaxY), gausime tokias formules:
          Re=MinX+x*(MaxX-MinX)/screenwidth;
          Im=MinY+y*(MaxY-MinY)/screenheight;
arba, kad išvengtume dažno dalijimo (kas užimtu daugiau laiko):
          dx= MaxX-MinX)/screenwidth;
          dy=(MaxY-MinY)/screenheight;
          Re=MinX+x*dx;
          Im=MinY+y*dy;

Gavome beveik baigtą programą:

PROGRAM Mandelbrot;
  CONST MinX=-2;
              MaxX=1.25;
              MinX=-1.25;
              MaxX=1.25;
  VAR dx,dy : REAL;
           x,y : INTEGER;
BEGIN
     dx= MaxX-MinX)/640;
     dy=(MaxY-MinY)/480;
     FOR y:=0 to 480-1 DO
     FOR x:=0 to 640-1 DO
         BEGIN
          color:=CALC_PIXEL(MinX+x*dx, MinY+y*dy);
          PUTPIXEL(X,Y, color);
        END;
END.

Galų gale, mūsų programa atrodo taip:

program Mandelbrot;
USES Crt, Graph;
Const MinX=-2;
         MaxX=1.25;
         MinY=-1.25;
         MaxY=1.25;
var dx,dy:real;
    x,y:integer;
    color:integer;
    screen_x,screen_y:integer;
    grDriver:integer;
    grMode:integer;
    ErrCode:integer;

Function calc_pixel(CA,CB:real):integer;
const max_iteration=64;
var
   old_a:real;
   iteration:integer;
   a,b:real;
   length_z:real;
begin
    a:=0;b:=0;
    iteration:=0;
    repeat
         old_a:=a;
         a:=a*a-b*b+ca;
         b:=2*old_a*b+cb;
         iteration:=iteration+1;
         length_z:=a*a+b*b;
     until (length_z>4) or (iteration>max_iteration);
     Calc_Pixel:=iteration;
End;

Begin
    grdriver:=Detect;
    InitGraph(grDriver,grMode,'c:\mokymas\tp\bgi\');
    ErrCode:=GraphResult;
    if ErrCode<>grOk then
      begin
        writeln('could not');
        Writeln('do you have correct..??');
        halt;
      end;
    screen_x:=getmaxx;
    screen_y:=getmaxy;
    dx:=(MaxX-MinX)/screen_x;
    dy:=(MaxY-MinY)/screen_y;
    for y:=0 to screen_y-1 do
      for x:=0 to screen_x-1 do
         begin
          color:=calc_pixel(MinX+x*dx,MinY+y*dy);
          putpixel(x,y,color);
        end;
     repeat until keypressed;
End.