selective_stdr.m

%% =====================================================
%  STDR Batch Simulation: Multi-Phase Selectivity Sweep
% =====================================================
clc; clear; close all;
model = 'STDR_3Phase_Selective';
load_system(model);

%% 1. Generate PN excitation
[pnA_vec, pnB_vec, pnC_vec, ~, Fs, ~] = generatePN3phase(1023, 20, 500e3, 'orthogonal');
Ts = 1/Fs; t = (0:length(pnA_vec)-1).' * Ts;
pnA = [t pnA_vec]; pnB = [t pnB_vec]; pnC = [t pnC_vec];
assignin('base','pnA',pnA); assignin('base','pnB',pnB); assignin('base','pnC',pnC);
v = 2e8;

%% 2. Define Experiment Matrix
phasesToFault = {'A', 'B', 'C'};
faultTypes = {'LG', 'SC', 'OC'}; % Varying fault characteristics
distances = [500, 1500, 2500];   % Test different lengths
noisePower = 1e-6;              % Constant low-level noise for realism

GammaMap = containers.Map({'LG','SC','OC'}, {-0.6, -1.0, +1.0});

%% 3. Result Storage
Results = {};
runID = 1;

%% 4. Triple Nested Selectivity Loop
for p = 1:numel(phasesToFault)
    targetPhase = phasesToFault{p};
    
    for f = 1:numel(faultTypes)
        faultType = faultTypes{f};
        Gamma = GammaMap(faultType);
        
        for d = 1:numel(distances)
            faultDistance = distances(d);
            Ndelay = round((2 * faultDistance / v) / Ts);
            
            fprintf('Run %02d | Target Phase: %s | Type: %s | Dist: %d m\n', ...
                    runID, targetPhase, faultType, faultDistance);
            
            %% --- SELECTIVE FAULT INJECTION ---
            % Reset all phases to healthy (Gain = 0)
            set_param([model '/Channel_3Phase/Gamma_A'], 'Gain', '0');
            set_param([model '/Channel_3Phase/Gamma_B'], 'Gain', '0');
            set_param([model '/Channel_3Phase/Gamma_C'], 'Gain', '0');
            set_param([model '/FaultDelayConst'], 'Value', num2str(Ndelay));
            
            % Inject fault only into the target phase
            set_param([model, '/Channel_3Phase/Gamma_', targetPhase], 'Gain', num2str(Gamma));
            
            % Update Noise on all lines
            set_param([model '/Channel_3Phase/Noise_A'], 'Cov', num2str(noisePower));
            set_param([model '/Channel_3Phase/Noise_B'], 'Cov', num2str(noisePower));
            set_param([model '/Channel_3Phase/Noise_C'], 'Cov', num2str(noisePower));
            
            % Run Simulation
            out = sim(model);
            
            % Extract results for all phases
            dA = out.dist_A.Data(end); dB = out.dist_B.Data(end); dC = out.dist_C.Data(end);
            
            % Store results
            Results(runID,:) = {targetPhase, faultType, faultDistance, dA, dB, dC};
            runID = runID + 1;
        end
    end
end

%% 5. Convert to Table
ResultsTable = cell2table(Results, 'VariableNames', ...
    {'FaultedPhase','FaultType','TrueDist','EstA','EstB','EstC'});
writetable(ResultsTable, 'STDR_Full_Selectivity_Study.csv');

%% 6. Visualization: Performance per Phase
figure('Color','w','Name','Cross-Phase Selectivity','Position', [100 100 1000 400]);

subplot(1,3,1); % Focus on Phase A detection
gscatter(ResultsTable.TrueDist, ResultsTable.EstA, ResultsTable.FaultedPhase);
title('Detector A Response'); xlabel('True Dist (m)'); ylabel('Est Dist (m)'); grid on;

subplot(1,3,2); % Focus on Phase B detection
gscatter(ResultsTable.TrueDist, ResultsTable.EstB, ResultsTable.FaultedPhase);
title('Detector B Response'); xlabel('True Dist (m)'); grid on;

subplot(1,3,3); % Focus on Phase C detection
gscatter(ResultsTable.TrueDist, ResultsTable.EstC, ResultsTable.FaultedPhase);
title('Detector C Response'); xlabel('True Dist (m)'); grid on;

saveas(gcf, 'STDR_CrossPhase_Selectivity.png');
disp('✔ Comprehensive Selectivity Sweep Complete.');