Welcome to Fractal Forums

function CalcDE(It3Dex: TPIteration3Dext; mctp: PMCTparameter): Double;
var dCT: TVec3D;
    Rtmp, RD1, RD2, RD3, Rst: Double;
    stmp, stmp2: Single;
    ftmp, itt, mi: Integer;
    CS, bIsC: LongBool;
begin
    if mctp.IsCustomDE then   //analytical DE method, based on one additional parameter
    with mctp^ do
    begin
      if DEoption = 20 then It3Dex.RStopD := msDEstop * StepWidth * 1.03;
      Result := mMandFunctionDE(@It3Dex.C1) * dDEscale;
    end
    else with It3Dex^ do
    begin                     //3d gradient
      CS := CalcSIT;
      mctp.mMandFunction(@C1);  
      CalcSIT := False;
      mi    := MaxIt;
      MaxIt := ItResultI;
      RStop := RStop * 8;
      Rtmp  := Rout;
      Rst   := mctp.mctDEoffset;
      mCopyVec(@dCT, @C1);
      mAddVecWeight(@C1, @mctp.Vgrads[2], Rst);
      mctp.mMandFunction(@C1);
      RD1 := Sqr(Rtmp - Rout);

      mCopyAddVecWeight(@C1, @dCT, @mctp.Vgrads[1], Rst);
      mctp.mMandFunction(@C1);
      RD2 := Sqr(Rtmp - Rout);

      mCopyAddVecWeight(@C1, @dCT, @mctp.Vgrads[0], Rst);
      mctp.mMandFunction(@C1);
      RD3 := Sqr(Rtmp - Rout);

      Result := Rtmp * Ln(Rtmp) * mctp.dDEscale / (Sqrt(RD1 + RD2 + RD3) + (Rst * 0.06));

      Rout := Rtmp;
      MaxIt := mi;
      RStop := mctp.dRStop;
      mCopyVec(@C1, @dCT);
      CalcSIT := CS;
    end;
    if Result < mctp.msDEstop * 0.25 then Result := mctp.msDEstop * 0.25;
    if mctp.FormulaType > 0 then   //DEcombinate, second formula
    begin
      itt := It3Dex.ItResultI;
      Rst := It3Dex.Otrap;
      stmp2 := It3Dex.SmoothItD;
      ftmp := mctp.FormulaType;
      mctp.FormulaType := 0;
      stmp := mctp.dDEscale;
      mctp.dDEscale := mctp.dDEscale2;
      bIsC := mctp.IsCustomDE;
      mctp.IsCustomDE := mctp.IsCustomDE2;
      It3Dex.nHybrid[0] := 999999;
      Rtmp := calcDE(It3Dex, mctp);
      case ftmp of
        1:  if Rtmp < Result then Result := Rtmp else begin It3Dex.ItResultI := itt; It3Dex.SmoothItD := stmp2; It3Dex.Otrap := Rst; end;
        2:  if Rtmp > Result then Result := Rtmp else begin It3Dex.ItResultI := itt; It3Dex.SmoothItD := stmp2; It3Dex.Otrap := Rst; end;
        3:  begin
              Result := Exp((Ln(Result) + Ln(Rtmp)) * 0.5);   
              It3Dex.ItResultI := Min(itt, It3Dex.ItResultI);
              It3Dex.SmoothItD := (It3Dex.SmoothItD + stmp2) * 0.5;
              It3Dex.Otrap := (It3Dex.Otrap + Rst) * 0.5;
            end;
      end;
      It3Dex.nHybrid[0] := 0;
      mctp.FormulaType := ftmp;
      mctp.dDEscale := stmp;
      mctp.IsCustomDE := bIsC;
    end;
end;


Raymarching part:


            RSFmul := 1;
            RLastStepWidth := dTmp * sZstepDiv;  //dTmp is the distance estimation, here on the first calculation point

            repeat

              if Iteration3Dext.ItResultI >= iMaxIt2 then // if inside while stepping, step back
              begin
                dT1 := -0.5 * RLastStepWidth;
                mZZ := mZZ + dT1;
                mAddVecWeight(@Iteration3Dext.C1, @VgradsFOV, dT1);
                msDEstop := DEstop * (1 + mZZ * mctDEstopFactor);
                dTmp := CalcDE(@Iteration3Dext, @MCTparas);
                RLastStepWidth := - dT1;
              end;

              if (Iteration3Dext.ItResultI < iMinIt) or
                ((Iteration3Dext.ItResultI < iMaxIt2) and (dTmp >= msDEstop)) then    //##### next step ######
              begin

                RLastDE := dTmp;

                dTmp := dTmp * sZstepDiv * RSFmul;
                dT1 := MaxCS(msDEstop, 0.4) * mctMH04ZSD;
                if dT1 < dTmp then
                begin
                  if not bFirstStep then StepCount := StepCount + dT1 / dTmp else StepCount := StepCount + Random;
                  dTmp := dT1;
                end
                else StepCount := StepCount + 1;
                
                if bFirstStep then
                begin
                  seed := 214013 * seed + 2531011;
                  bFirstStep := False;
                  dTmp := ((seed shr 16) and $7FFF) * 0.000030517578125 * dTmp;
                end;

                RLastStepWidth := dTmp;
                mZZ := mZZ + dTmp;
                mAddVecWeight(@Iteration3Dext.C1, @VgradsFOV, dTmp);
                msDEstop := DEstop * (1 + mZZ * mctDEstopFactor);

                dTmp := CalcDE(@Iteration3Dext, @MCTparas);

                if dTmp > RLastDE + RLastStepWidth then dTmp := RLastDE + RLastStepWidth;  //helps a little on julias etc

                if RLastDE > dTmp + s1em30 then  //a little down regulation
                begin
                  dT1 := RLastStepWidth / (RLastDE - dTmp);
                  if dT1 < 1 then
                    RSFmul := maxCS(0.5, dT1)
                  else
                    RSFmul := 1;
                end
                else RSFmul := 1;

              end
...
            until (mZZ > Zend) or PCalcThreadStats.pLBcalcStop^;