unit PotentialU;
//this procedure should draw Hubbard-Duaddy potential of Mandelbrot set
// but did not
//  I dont know what it is , but it is nice (:-))
// see http://rgba.scenesp.org/iq/trastero/potential/
// for correct images by Inigo Quilez
//
// Inigo Quilez have found errors:
// 1-st in Procedure CmplxDiv
// second: Z0 should be 0 not c


interface
uses  graphics, // TColor
sysutils, //FDloatToStr
dialogs, // showmessage
        windows; // GetBValue

Procedure DrawPotential_M;

implementation
uses Math,  {SameValue,InRange}
     MainFormU, C_PlaneU, AlgorithmDlgU, BitmapU, MappingU, MathU;

//--------------------------------
Procedure DrawPotential_M;
  
  var iteration:integer;
     iY,iX:integer;
     eCy ,eCx:extended; // C:=eCx + eCy*i
     temp:extended;
     Zx,Zy:extended;
     //potential
     dx,dy,
     fac:extended;
     G:extended;  // Green function
     bG:byte;
     //
  begin
    For iY :=iYmin to iYMax do
      begin
        eCy:=Convert_iY_to_eY(iY);
        for iX:=iXmin to iXmax do
          begin
            eCx:=Convert_iX_to_eX(iX);
            //F(0)=c   first iteration
            Zx:=eCx; // error 
            Zy:=eCy;
            // for Green function
            G :=Ln(Hypot(eCx,eCy));
            fac :=0.5;

            //if PointIsInCardioid(eCx,eCy) or PointIsInComponent(eCx,eCy)
              //then Iteration:=IterationMax
              //else
               For iteration:=1 to iterationMax do
                      begin
                        // Z(n+1):= Sqr(Zn) + C
                        temp:=Sqr(Zx)-Sqr(Zy) +eCx;
                        Zy:=2*Zx*Zy+eCy;
                        Zx:=temp;
                        // d:=1+c/(z*z)
                        GreenSum (eCx,eCy,Zx,Zy,dx,dy);
                        G:= G + fac*ln(hypot(dx,dy));
                        //fac:=fac*0.5;
                        // instead check Sqrt(Sqr(Zx))>bailOut check Sqr(Zx)>sqr(bailOut))
                        // ************* bailout test ***************************888
                        if (Sqr(Zx)+Sqr(Zy))> bailout2    then break;
                      end;
              // potential
              //G:= + sum;

            // drawing procedure
            if iteration>=IterationMax
              then with Bitmap.FirstLine[iY*Bitmap.LineLength+iX] do
                  begin
                    B := 0;
                    G := 0;
                    R := 0;
                    //A := 0;
                  end
              //--------------

              else
                begin
                  bG:= round(g*25)  ;
                  with Bitmap.FirstLine[iY*Bitmap.LineLength+iX] do
                              begin
                                B := bG;
                                G := bG;
                                R := Bg;
                                //A := 0;
                              end ;
                end; // else
              // ----------
      end; // for iX
    end; //For iY

  end;
//------------------------
//********************************************************************************************
  Procedure CplxSqrt(Xin,Yin:extended; var Xout,Yout:extended);
     var r:extended;
     begin
      r:=Sqrt(Xin*Xin + Yin*Yin);
      Xout:=Sqrt(0.5 * (r + Xin));
      if Yin< 0.0 then Yout:= - Sqrt(0.5 * (r - Xin))
                  else Yout:=   Sqrt(0.5 * (r - Xin));
     end;



//----------------------------------------------
Procedure CmplxSum(X1in,Y1in,X2in,Y2in:extended; var Xout,Yout:extended);
begin
  Xout:=X1in+X2in;
  Yout:=Y1in+Y2in;

end;
//-------------------
Procedure CmplxDiv(X1in,Y1in,X2in,Y2in:extended; var Xout,Yout:extended);
 var denominator:extended;
 begin
  denominator:= sqr(X2in) + sqr(Y2in);
  Xout:=(X1in*X2in+Y1in*Y2in) / denominator;
  Yout:=(X2in*Y1in - X1in * Y2in)/ denominator; //error !!!!!!!!!!!!
end;
//--------------------------------------------
Procedure CmplxSqr(Xin,Yin:extended; var Xout,Yout:extended);
  begin
    Xout:=(Xin-Yin)*(Xin+Yin);  // (x*x) - (y*y) = (x-y)*(x+y)
    Yout:=2* Xin * Yin;
  end;
//---------------------------------
procedure GreenSum (Cx,Cy,Zx,Zy:extended; var Xout,Yout:extended);
  // = 1 + (c/z2)
  var xTemp,yTemp,Z2x,Z2y :extended;

  begin
  // z2 = z*z
   CmplxSqr(Zx,Zy,Z2x,Z2y);
   // temp = c /z2
   CmplxDiv(Cx,Cy,Z2x,Z2y,xTemp,yTemp);
   // out = 1+ temp
   Xout:=1 +xTemp;
   Yout:=yTemp;


  end;
//------------------



end.