
% Matlab script that compares serpent output files and prints statistically
% significant differences.

% Usage: read two or more res files, set variable ref to reference case
%        index and run script.

printf("\nComparing Serpent output files. The relative difference in each ");
printf("parameter is compared to the 2*sigma confidence interval of\n");
printf("the parameter with the largest error. Reference case = %d.\n\n", ref);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define a function to compare values:    

function count = Compare (param, pname, ref)

% Reset counters:

count = [0 0];

% Get size:

[nc nv] = size(param);

% Loop over values:

for j = 1:2:nv;

% Loop over cases:

for i = 1:nc;

if (i != ref);

if (param(ref,j) != 0.0)

% Calculate relative difference:

d = param(i,j)/param(ref,j) - 1.0;

% Calculate combined error

e = 2*max(param(i, j + 1), param(ref, j + 1));

% Add to total counter:

count(1) = count(1) + 1;

% Compare:

if (e > 0);
if (param(i,j)*param(ref,j) > 0.0);
if (abs(d) > e);

% Print:

printf("%s (%3d, %3d): ", pname, i, j);

printf("%12.5E (%1.5f) vs. %12.5E (%1.5f): diff = %9.5f (%4.1f x 2*sigma)\n", param(ref,j), param(ref,j+1), param(i,j), param(i,j+1), d, abs(d)/e); 

% Add to error counter:

count(2) = count(2) + 1;

end;
end;
end;

end;

end;

end;

end;
    
endfunction


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Reset counter:

count = [0 0];

% Slowing-down and thermalization:

count = count + Compare(COL_SLOW                  , 'COL_SLOW                  ', ref);
count = count + Compare(COL_THERM                 , 'COL_THERM                 ', ref);
count = count + Compare(COL_TOT                   , 'COL_TOT                   ', ref);
count = count + Compare(SLOW_TIME                 , 'SLOW_TIME                 ', ref);
count = count + Compare(THERM_TIME                , 'THERM_TIME                ', ref);
count = count + Compare(SLOW_DIST                 , 'SLOW_DIST                 ', ref);
count = count + Compare(THERM_DIST                , 'THERM_DIST                ', ref);
count = count + Compare(THERM_FRAC                , 'THERM_FRAC                ', ref);

% Criticality eigenvalues:

count = count + Compare(ANA_KEFF                  , 'ANA_KEFF                  ', ref);
count = count + Compare(IMP_KEFF                  , 'IMP_KEFF                  ', ref);
count = count + Compare(COL_KEFF                  , 'COL_KEFF                  ', ref);
count = count + Compare(ABS_KEFF                  , 'ABS_KEFF                  ', ref);
count = count + Compare(ABS_KINF                  , 'ABS_KINF                  ', ref);
count = count + Compare(ABS_GC_KINF               , 'ABS_GC_KINF               ', ref);
count = count + Compare(IMPL_ALPHA_EIG            , 'IMPL_ALPHA_EIG            ', ref);
count = count + Compare(FIXED_ALPHA_EIG           , 'FIXED_ALPHA_EIG           ', ref);

% Normalization:

count = count + Compare(TOT_POWER                 , 'TOT_POWER                 ', ref);
count = count + Compare(TOT_GENRATE               , 'TOT_GENRATE               ', ref);
count = count + Compare(TOT_FISSRATE              , 'TOT_FISSRATE              ', ref);
count = count + Compare(TOT_ABSRATE               , 'TOT_ABSRATE               ', ref);
count = count + Compare(TOT_LEAKRATE              , 'TOT_LEAKRATE              ', ref);
count = count + Compare(TOT_LOSSRATE              , 'TOT_LOSSRATE              ', ref);
count = count + Compare(TOT_FLUX                  , 'TOT_FLUX                  ', ref);
count = count + Compare(TOT_RR                    , 'TOT_RR                    ', ref);
count = count + Compare(TOT_SOLU_ABSRATE          , 'TOT_SOLU_ABSRATE          ', ref);
count = count + Compare(TOT_POWDENS               , 'TOT_POWDENS               ', ref);

% Point-kinetic parameters:

count = count + Compare(ANA_PROMPT_LIFETIME       , 'ANA_PROMPT_LIFETIME       ', ref);
count = count + Compare(IMPL_PROMPT_LIFETIME      , 'IMPL_PROMPT_LIFETIME      ', ref);
count = count + Compare(ANA_REPROD_TIME           , 'ANA_REPROD_TIME           ', ref);
count = count + Compare(IMPL_REPROD_TIME          , 'IMPL_REPROD_TIME          ', ref);
count = count + Compare(DELAYED_EMTIME            , 'DELAYED_EMTIME            ', ref);

% Six-factor formula:

count = count + Compare(SIX_FF_ETA                , 'SIX_FF_ETA                ', ref);
count = count + Compare(SIX_FF_F                  , 'SIX_FF_F                  ', ref);
count = count + Compare(SIX_FF_P                  , 'SIX_FF_P                  ', ref);
count = count + Compare(SIX_FF_EPSILON            , 'SIX_FF_EPSILON            ', ref);
count = count + Compare(SIX_FF_LF                 , 'SIX_FF_LF                 ', ref);
count = count + Compare(SIX_FF_LT                 , 'SIX_FF_LT                 ', ref);
count = count + Compare(SIX_FF_KINF               , 'SIX_FF_KINF               ', ref);
count = count + Compare(SIX_FF_KEFF               , 'SIX_FF_KEFF               ', ref);

% Delayed neutron parameters:

count = count + Compare(BETA_EFF                  , 'BETA_EFF                  ', ref);
count = count + Compare(BETA_ZERO                 , 'BETA_ZERO                 ', ref);
count = count + Compare(DECAY_CONSTANT            , 'DECAY_CONSTANT            ', ref);

% Few-group cross sections:

count = count + Compare(FLUX                      , 'FLUX                      ', ref);
count = count + Compare(LEAK                      , 'LEAK                      ', ref);
count = count + Compare(TOTXS                     , 'TOTXS                     ', ref);
count = count + Compare(FISSXS                    , 'FISSXS                    ', ref);
count = count + Compare(CAPTXS                    , 'CAPTXS                    ', ref);
count = count + Compare(ABSXS                     , 'ABSXS                     ', ref);
count = count + Compare(ELAXS                     , 'ELAXS                     ', ref);
count = count + Compare(INELAXS                   , 'INELAXS                   ', ref);
count = count + Compare(SCATTXS                   , 'SCATTXS                   ', ref);
count = count + Compare(N2NXS                     , 'N2NXS                     ', ref);
count = count + Compare(REMXS                     , 'REMXS                     ', ref);
count = count + Compare(NUBAR                     , 'NUBAR                     ', ref);
count = count + Compare(NSF                       , 'NSF                       ', ref);
count = count + Compare(RECIPVEL                  , 'RECIPVEL                  ', ref);
count = count + Compare(FISSE                     , 'FISSE                     ', ref);

% Fission spectra:

count = count + Compare(CHI                       , 'CHI                       ', ref);
count = count + Compare(CHIP                      , 'CHIP                      ', ref);
count = count + Compare(CHID                      , 'CHID                      ', ref);

% Group-transfer probabilities and cross sections:

count = count + Compare(GTRANSFP                  , 'GTRANSFP                  ', ref);
count = count + Compare(GTRANSFXS                 , 'GTRANSFXS                 ', ref);

% Diffusion parameters:

count = count + Compare(DIFFAREA                  , 'DIFFAREA                  ', ref);
count = count + Compare(DIFFCOEF                  , 'DIFFCOEF                  ', ref);
count = count + Compare(TRANSPXS                  , 'TRANSPXS                  ', ref);
count = count + Compare(MUBAR                     , 'MUBAR                     ', ref);

% PN scattering cross sections:

count = count + Compare(SCATT0                    , 'SCATT0                    ', ref);
count = count + Compare(SCATT1                    , 'SCATT1                    ', ref);
count = count + Compare(SCATT2                    , 'SCATT2                    ', ref);
count = count + Compare(SCATT3                    , 'SCATT3                    ', ref);
count = count + Compare(SCATT4                    , 'SCATT4                    ', ref);
count = count + Compare(SCATT5                    , 'SCATT5                    ', ref);

% P1 diffusion parameters:

count = count + Compare(P1_TRANSPXS               , 'P1_TRANSPXS               ', ref);
count = count + Compare(P1_DIFFCOEF               , 'P1_DIFFCOEF               ', ref);
count = count + Compare(P1_MUBAR                  , 'P1_MUBAR                  ', ref);

% Assembly discontinuity factors:

count = count + Compare(ADFS                      , 'ADFS                      ', ref);
count = count + Compare(ADFC                      , 'ADFC                      ', ref);

% Power distributions in lattices:

count = count + Compare(POWDISTR10                , 'POWDISTR10                ', ref);

printf("\nComparison failed %d / %d times (%1.1f%%).\n\n", count(2), count(1), 100*(count(2)/count(1)));

