

%% 0) Inicializace pralelních simulací
clear; clc; close all;
if isempty(gcp('nocreate'))
    parpool;                              % Spustí paralelní pool
end
dq = parallel.pool.DataQueue;
% po každém send(dq,...) se vytiskne zpráva na konzoli
afterEach(dq, @(msg) fprintf('%s\n', msg));


%% 1) Definice parametrů simulace a konstant
numSims    = 24;                          % počet paralelních běhů
Ip         = 20;                          % špičkový proud blesku [kA]

% Výpočet potenciálu horní elektrody z Bazelyanovy rovnice
V0_MV      = 3.6858 + 0.7282*Ip - 0.000818*Ip^2;  
V_up       = -V0_MV * 1e6;                % potenciál horní elektrody [V]
V_lo       = 0;                           % potenciál země [V]

% Geometrie mřížky
N          = 100;                         % počet bodů na osách
margin     = 15;                          % okraj kolem oblastí zájmu
xseed      = 50;                           % počáteční x-pozice leaderu
x1_up = -150; x2_up = 150; y_up = 300;     % poloha horní elektrody
y_lo  = 0;                                % poloha země

% souřadnice bodových emitérů
xbod  = [5.7, -5.7, -10.7, 10.7, -14.5, 14.5];  
ybod  = [41.5,41.5,30.519,30.519,19.019,19.019];  
M      = numel(xbod);                    % jejich počet

% vektor x, y a mřížka XX,YY
x_min = x1_up - margin;  x_max = x2_up + margin;
y_min = y_lo  - 5;       y_max = max([y_up,ybod]) + 5;
x = linspace(x_min, x_max, N);
y = linspace(y_min, y_max, N);
[XX,YY] = meshgrid(x,y);
dx = x(2) - x(1);                        % krok mřížky
dy = y(2) - y(1);

% Příprava Dirichletových okrajů (horní elektroda a zem)
th      = max(dx,dy);
mask_up = (abs(YY - y_up) < th/2) & (XX >= x1_up) & (XX <= x2_up);
mask_lo = abs(YY - y_lo) < th/2;

% Masky pro bodové emitéry (bodové kanály)
mask_pt = false(N,N,M);
for m = 1:M
    [~, ix] = min(abs(x - xbod(m)));
    [~, iy] = min(abs(y - ybod(m)));
    mask_pt(iy,ix,m) = true;
end

% Parametry DBM modelu
maxIter    = 1100;       % max. iterací SOR
tol        = 1e-4;       % tolerance konvergence
omega      = 1.8;        % relaxační faktor SOR
Ecr_neg    = 1e6;        % kritické pole pro negativní leader
Ecr_pos    = 500e3;      % kritické pole pro pozitivní strimér
eta_neg    = 8;          % exponent pro váhy negativního leaderu
eta_gnd    = 8;          % exponent pro ground-plane strimér
eta_pt     = 8;          % exponent pro vzestupné striméry
sigma_frac = 0.7;        % šířka rozložení pravděpodobnosti
nmaxSteps  = 2000;       % maximální počet kroků růstu

% Směry růstu: indexy posunu pro negativní leader a pozitivní strimér
offsets_neg = [ 0 -1; 0 1; -1 0; -1 -1; -1 1 ];
offsets_pos = [ 0 -1; 0 1; 1 0; 1 -1; 1 1 ];


%% 2) Inicializace výsledků
Larray     = zeros(N, N, numSims, 'uint8');  % třírozměrná matice fáze
connResult = strings(numSims, 1);            % typ spojení pro každý běh


%% 3) Paralelní smyčky
parfor sim = 1:numSims
    send(dq, sprintf('Sim %d startuje', sim));

    % reset masek a stavů pro tento běh
    mask_neg       = false(N);          % negativní leader maska
    mask_pos_gnd   = false(N);          % ground-plane strimér
    mask_pos_pt    = false(N,N,M);      % bodové striméry
    started_pt     = false(M,1);        % které body již aktivovány
    ground_started = false;             % zda byl aktivován zemní strimér
    connection     = "";                % typ dosaženého spojení

    % Místo vzniku sestupného výboje
    [~, ix0] = min(abs(x - xseed));
    [~, iy0] = min(abs(y - y_up));
    mask_neg(iy0, ix0) = true;

    % Indexace vodičů
    ix_pt = zeros(M,1); iy_pt = zeros(M,1);
    for m = 1:M
        [~, ix_pt(m)] = min(abs(x - xbod(m)));
        [~, iy_pt(m)] = min(abs(y - ybod(m)));
    end

    % hlavní smyčka růstu
    for step = 1:nmaxSteps
        % výpis každých 100 kroků
        if mod(step,100) == 0
            send(dq, sprintf('Sim %d: krok %d/%d', sim, step, nmaxSteps));
        end

        % vytvoříme masku okolí elektrod
        fixed = mask_up | mask_lo | any(mask_pt,3) | mask_neg ...
                | mask_pos_gnd | any(mask_pos_pt,3);

        %% 3.1) Sestavení potenciálové matice V
        V = zeros(N);
        V(mask_up)   = V_up;      % horní elektroda
        V(mask_lo)   = V_lo;      % zem
        for m = 1:M
            V(mask_pt(:,:,m)) = V_lo;   % bodoví emitory na zemi
        end
        V(mask_neg)  = V_up;     % negativní leader má stejný potenciál jako elektroda
        if ground_started
            V(mask_pos_gnd) = V_lo;    % zemní strimér na zemi
        end
        for m = 1:M
            if started_pt(m)
                V(mask_pos_pt(:,:,m)) = V_lo;  % aktivní bodové strimér
            end
        end

        %% 3.2) SOR‐solver pro Laplaceovu rovnici
        for it = 1:maxIter
            V_old = V;
            % aktualizace každého nenastaveného bodu
            for j = 2:N-1
              for i = 2:N-1
                if fixed(j,i), continue; end
                Vnew    = 0.25*(V(j+1,i) + V(j-1,i) + V(j,i+1) + V(j,i-1));
                V(j,i)  = V(j,i) + omega*(Vnew - V(j,i));
              end
            end
            % Neumannovy okraje (zero-normální derivace)
            V(1,:)   = V(2,:);
            V(end,:) = V(end-1,:);
            V(:,1)   = V(:,2);
            V(:,end) = V(:,end-1);
            if max(abs(V(:) - V_old(:))) < tol
                break;
            end
        end

        %% 3.3) Výpočet elektrického pole a spuštění strimérů
        [Vx, Vy] = gradient(V, dx, dy);
        Emag     = sqrt(Vx.^2 + Vy.^2);

        % spustíme zemní (ground-plane) strimér, až pole překročí Ecr_pos
        if ~ground_started && max(Emag(mask_lo)) >= Ecr_pos
            ground_started = true;
            tmp = Emag; tmp(~mask_lo) = 0;
            [~, idx] = max(tmp(:));
            [jr, ir] = ind2sub([N,N], idx);
            mask_pos_gnd(jr, ir) = true;
        end

        % spustíme bodové striméry, pokud na bodě přesáhlo Ecr_pos
        for m = 1:M
            if ~started_pt(m) && Emag(iy_pt(m), ix_pt(m)) >= Ecr_pos
                started_pt(m) = true;
                mask_pos_pt(iy_pt(m), ix_pt(m), m) = true;
            end
        end

        %% 3.4) Sestavení front: kde se může růst
        % fronta negativního leaderu
        if step == 1
            front_neg = false(N);
            jdown = iy0 - 1;
            if jdown >= 1 && ~fixed(jdown, ix0)
                front_neg(jdown, ix0) = true;
            end
        else
            front_neg = false(N);
            [rN, cN] = find(mask_neg);
            for k = 1:numel(rN)
              for o = 1:size(offsets_neg,1)
                rr = rN(k) + offsets_neg(o,1);
                cc = cN(k) + offsets_neg(o,2);
                if rr>=1 && rr<=N && cc>=1 && cc<=N && ~fixed(rr,cc)
                    front_neg(rr,cc) = true;
                end
              end
            end
        end

        % fronta zemního striméru
        front_gnd = false(N);
        if ground_started
            [rG, cG] = find(mask_pos_gnd);
            for k = 1:numel(rG)
              for o = 1:size(offsets_pos,1)
                rr = rG(k) + offsets_pos(o,1);
                cc = cG(k) + offsets_pos(o,2);
                if rr>=1 && rr<=N && cc>=1 && cc<=N && ~fixed(rr,cc)
                    front_gnd(rr,cc) = true;
                end
              end
            end
        end

        % fronty pro každý bodový strimér
        front_pt = false(N,N,M);
        for m = 1:M
          if started_pt(m)
            [rP, cP] = find(mask_pos_pt(:,:,m));
            for k = 1:numel(rP)
              for o = 1:size(offsets_pos,1)
                rr = rP(k) + offsets_pos(o,1);
                cc = cP(k) + offsets_pos(o,2);
                if rr>=1 && rr<=N && cc>=1 && cc<=N && ~fixed(rr,cc)
                  front_pt(rr,cc,m) = true;
                end
              end
            end
          end
        end

        %% 3.5) Probabilistický růst na frontě
        % negativní leader
        En = Emag; En(~front_neg) = 0;
        Emax = max(En(:));
        sigma = sigma_frac * Emax;
        Wn   = exp(-(Emax - En).^2 / (2*sigma^2));
        Pn   = Wn.^eta_neg; Pn(~front_neg) = 0;
        Pn   = Pn / sum(Pn(:));
        idx  = find(rand <= cumsum(Pn(:)),1);
        [jn, in_] = ind2sub([N,N], idx);
        mask_neg(jn, in_) = true;

        % zemní strimér
        if ground_started
            Eg = Emag; Eg(~front_gnd) = 0;
            Emax_g = max(Eg(:)); sigma_g = sigma_frac * Emax_g;
            Wg = exp(-(Emax_g - Eg).^2 / (2*sigma_g^2));
            Pg = Wg.^eta_gnd; Pg(~front_gnd) = 0;
            Pg = Pg / sum(Pg(:));
            idxG = find(rand <= cumsum(Pg(:)),1);
            [jg, ig] = ind2sub([N,N], idxG);
            mask_pos_gnd(jg, ig) = true;
        end

        % bodové strimér
        for m = 1:M
            if started_pt(m)
                Ep = Emag; Ep(~front_pt(:,:,m)) = 0;
                Emax_p = max(Ep(:)); sigma_p = sigma_frac * Emax_p;
                Wp = exp(-(Emax_p - Ep).^2 / (2*sigma_p^2));
                Pp = Wp.^eta_pt; Pp(~front_pt(:,:,m)) = 0;
                Pp = Pp / sum(Pp(:));
                idxP = find(rand <= cumsum(Pp(:)), 1);
                [jp, ip] = ind2sub([N,N], idxP);
                mask_pos_pt(jp, ip, m) = true;
            end
        end

        %% 3.6) Kontrola přímého spojení mezi kanály
        adj_gnd = false;
        [rN, cN] = find(mask_neg);
        for k = 1:numel(rN)
          for o = 1:size(offsets_neg,1)
            rr = rN(k) + offsets_neg(o,1);
            cc = cN(k) + offsets_neg(o,2);
            if rr>=1 && rr<=N && cc>=1 && cc<=N && mask_pos_gnd(rr,cc)
              adj_gnd = true; break;
            end
          end
          if adj_gnd, break; end
        end

        adj_pt = false(M,1);
        for m = 1:M
          for k = 1:numel(rN)
            for o = 1:size(offsets_neg,1)
              rr = rN(k) + offsets_neg(o,1);
              cc = cN(k) + offsets_neg(o,2);
              if rr>=1 && rr<=N && cc>=1 && cc<=N && mask_pos_pt(rr,cc,m)
                adj_pt(m) = true; break;
              end
            end
            if adj_pt(m), break; end
          end
        end

        % určíme typ spojení a přerušíme růstovou smyčku
        if adj_gnd || any(adj_pt)
          if adj_gnd && ~any(adj_pt)
            connection = 'ground-plane strimér';
          elseif ~adj_gnd && any(adj_pt)
            m0 = find(adj_pt,1);
            connection = sprintf('point %d strimér', m0);
          else
            E_vals = []; names = {};
            if adj_gnd
              vals = []; 
              for k = 1:numel(rN)
                for o = 1:size(offsets_pos,1)
                  rr = rN(k)+offsets_pos(o,1);
                  cc = cN(k)+offsets_pos(o,2);
                  if mask_pos_gnd(rr,cc)
                    vals(end+1) = Emag(rr,cc);
                  end
                end
              end
              E_vals(end+1)   = max(vals); 
              names{end+1}    = 'ground-plane strimér';
            end
            for m = 1:M
              if adj_pt(m)
                vals = [];
                for k = 1:numel(rN)
                  for o = 1:size(offsets_pos,1)
                    rr = rN(k)+offsets_pos(o,1);
                    cc = cN(k)+offsets_pos(o,2);
                    if mask_pos_pt(rr,cc,m)
                      vals(end+1) = Emag(rr,cc);
                    end
                  end
                end
                E_vals(end+1)   = max(vals);
                names{end+1}    = sprintf('point %d strimér',m);
              end
            end
            [~, iMax]   = max(E_vals);
            connection  = names{iMax};
          end
          break;
        end
    end  % konec růstové smyčky

    %% 4) Uložení výsledků pro danou simulaci
    L = zeros(N,'uint8');
    L(mask_neg)       = 1;        % negativní leader
    L(mask_pos_gnd)   = 2;        % ground-plane strimér
    for m = 1:M
      L(mask_pos_pt(:,:,m)) = 2 + m;  % body + jejich striméry
    end
    Larray(:,:,sim)  = L;
    connResult(sim)  = connection;
    send(dq, sprintf('Sim %d dokončena → %s', sim, connection));
end  % konec parfor

%% 5) Vyhodnocení výsledků
count_gnd = sum(connResult == 'ground-plane strimér');
count_pt  = zeros(M,1);
for m = 1:M
  count_pt(m) = sum(connResult == sprintf('point %d strimér', m));
end
fprintf('Po %d simulacích: ground-plane=%d, points=', numSims, count_gnd);
fprintf('%d ', count_pt); 
fprintf('\n');

%% 6) Vykreslení výsledků (skryté obrázky)
% Nastavení colormap: bílá, modrá, zelená pro body 1&2 a červená pro ostatní
cmap = [
  1 1 1;                         % pozadí
  0 0 1;                         % negativní leader – modrá
  repmat([0 1 0], 3, 1);         % ground+body1,2 – zelená
  repmat([1 0 0], M-2, 1)        % ostatní body – červená
];

figHandles = gobjects(numSims,1);
for sim = 1:numSims
  figHandles(sim) = figure('Visible','off');
  imagesc(x, y, Larray(:,:,sim));
  axis xy equal tight;
  colormap(cmap);
  caxis([-0.5, 2 + M + 0.5]);
  hold on;

  % horní elektroda (černá čára)
  plot([x1_up x2_up], [y_up y_up], 'k-', 'LineWidth', 2);
  % zem (černá čára)
  plot([x_min x_max], [y_lo y_lo], 'k-', 'LineWidth', 2);

  % vykreslení bodových emitérů
  for m = 1:M
    mk = 'o'; if mod(m,2)==0, mk = 's'; end
    if m <= 2
      col = 'g';  % první dva body zelené
    else
      col = 'r';  % ostatní červené
    end
    plot(xbod(m), ybod(m), mk, ...
         'MarkerEdgeColor', col, 'MarkerFaceColor', col, ...
         'DisplayName', sprintf('pt%d', m));
  end

  title(sprintf('Sim %d: %s', sim, connResult(sim)));
  xlabel('x [m]'); ylabel('y [m]');
end
