function Text = corona3 (Index_min, Index_max, CoronaDatei, Modus)
% ------------------------------------------------------------------------------
%  SYNTAX
%Text = corona3 (Index_min, Index_max, CoronaDatei, Modus);
% ------------------------------------------------------------------------------
%  DESCRIPTION
%    Dokumentation des Verlaufs der Infektionen mit dem Corona-Virus SARS-CoV-2
%    in Deutschland
%
%    Fallzahl-Daten von https://npgeo-corona-npgeo-de.hub.arcgis.com/datasets/
%                                             dd4580c810204019a7b8eb3e0b329dd6_0
% ------------------------------------------------------------------------------
%  INPUT
%    Index_min (optional): Index_min > 0: Daten ab Index_min auswerten
%                          sonst: Corona-Datei neu formatieren & speichern
%    Index_max (optional): Index_min < Index_max: Daten von Index_min bis
%                                                 Index_max auswerten
%                          sonst: Corona-Datei neu formatieren & speichern
%    CoronaDatei (optional): Name der Datei mit den Corona-Daten
%    Modus (optional):     Modus == 1: Auswertung mit Meldedatum
%                          Modus == 2: Auswertung mit Erkrankungsdatum ohne NC
%                          sonst: Auswertung mit Erkrankungsdatum und Nowcasting
% ------------------------------------------------------------------------------
%  OUTPUT
%    Text: Berechnete Parameter
%
%    Plots mit Fallzahlen und Corona-Pandemie-Parametern
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%  CHANGE TABLE
%  Date       Name       Modification
%  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%  02.08.2020 JMH        Initiale Version basierend auf corona.m und corona2.m
%  03.08.2020 JMH        Weiter aufgeräumt (filtereRohdaten)
%  04.08.2020 JMH        Weiter aufgeräumt (berechneWendepunkt)
%  05.08.2020 JMH        Weiter aufgeräumt (Plot 1)
%  06.08.2020 JMH        Neues Corona-Verzeichnis "/home/jmh/octave/corona"
%  07.08.2020 JMH        Umstellung von zeit.de-Daten auf RKI-Originaldaten
%  09.08.2020 JMH        Weiter aufgeräumt (Tägliche Neuinfektionen und Tote)
%  10.08.2020 JMH        Weiter aufgeräumt (Plot 2)
%  17.08.2020 JMH        Weiter aufgeräumt (Reproduktionszahl R berechnen)
%  17.08.2020 JMH        Berechnung Tägliche Neuinfektionen und Tote korrigiert
%  18.08.2020 JMH        Weiter aufgeräumt (Plot 3)
%  18.08.2020 JMH        berechneWendepunkt repariert
%  19.08.2020 JMH        Berechnung Tägliche Neuinfektionen und Tote korrigiert
%  19.08.2020 JMH        Berechnung Reproduktionszahl R korrigiert
%  19.08.2020 JMH        Plots 2 und 3 korrigiert
%  27.08.2020 CN         Funktion thin_labels ergänzt
%  30.08.2020 JMH        Berechnung Reproduktionszahl R weiter korrigiert
%  30.08.2020 JMH        TagFlag wieder eingeführt
%  30.08.2020 JMH        Referenz-Zeitpunkt der Reproduktionszahl R korrigiert
%  30.08.2020 JMH        Berechnung der 7-Tage-Mittelwerte korrigiert
%  31.08.2020 JMH        Schwelle bei Wendepunkt-Berechnung
%  31.08.2020 JMH        Legenden von Plot 1 besser positioniert
%  31.08.2020 CN         Legende von Plot 2 besser positioniert
%  05.09.2020 CN         "skaliere_Wochen": Achsen-Label verbessert
%  05.09.2020 CN         "thin_labels": keine Aktion, wenn schon ausgedünnt ist
%  05.09.2020 CN         Legende von Plot 2 besser positioniert
%  05.09.2020 JMH        Aufgeräumt
%  12.09.2020 JMH        Fehlermeldungen verbessert
%  15.09.2020 CN         Aktuelle Sterberate korrigiert
%  17.09.2020 JMH        Abstand der Wendepunkte zu den Rändern >= 1 Woche!
%  18.09.2020 CN         "scale_axis" korrigiert
%  18.09.2020 CN         Farben und Legendenposition von Plot 2 verschönert
%  18.09.2020 CN         y-Achse von Plot 3 verschönert
%  18.09.2020 CN         "skaliere_Wochen" korrigiert
%  18.09.2020 JMH        Aufgeräumt
%  24.09.2020 JMH        MAT-Datei mit akkumulierten Ergebnissen von corona3.m
%  26.09.2020 JMH        CoronaDatenRKI.mat in Archiv-Verzeichnis verschoben
%  27.09.2020 JMH        Linienbreiten von Plots 2 und 3 korrigiert
%  12.10.2020 JMH        Verdoppelungszeit im Titel von Plot 3 ausgeben
%  22.10.2020 JMH        Legenden von Plots 2 und 3 besser positioniert
%  24.10.2020 JMH        Funktion thin_labels korrigiert
%  26.10.2020 CN         "berechneWendepunkt" erweitert auf mehrere Wendepunkte
%  27.10.2020 JMH        Legende von Plot 3 besser positioniert
%  27.10.2020 CN         "get_ticklabel" und "set_ticklabel" korrigiert
%  29.10.2020 JMH        Aufgeräumt
%  30.10.2020 JMH        "bestimme_Wendepunkte" korrigiert (Extremwertsuche)
%  30.10.2020 JMH        "Lupe" im unteren Subplot von Plot 1 immer angeschaltet
%  30.10.2020 JMH        "Neuerkrankte" und "Erkrankte" => "Infizierte"
%  30.10.2020 JMH        Legenden verschönert
%  31.10.2020 JMH        Plot 1 repariert
%  02.11.2020 JMH        Plots korrigiert: '.-' => '-'
%  03.11.2020 JMH        Fiesen Fehler bei "skaliere_Wochen" behoben
%  03.11.2020 JMH        Plot 2 repariert
%  21.11.2020 JMH        return => end
%  21.11.2020 JMH        Schwelle bei Wendepunkt-Berechnung wieder auf 1000
%  22.11.2020 JMH        Kranke_MD und Kranke_oNC ergänzt (Einlesen & MAT-Datei)
%  23.11.2020 JMH        TagFlag eliminiert
%  23.11.2020 JMH        Kranke_MD und Kranke_oNC ergänzt (Plots)
%  25.11.2020 JMH        "eval" eliminiert
%  01.01.2021 JMH        x-Achse: Kalenderwochen in 2020/2021
%  04.01.2021 JMH        Oberen Subplot von Plot 1 korrigiert
%  04.01.2021 JMH        Hub-Schwelle für Wendepunkt-Berechnung abgesenkt: 0.28
%  04.01.2021 JMH        Legenden besser positioniert
%  04.01.2021 JMH        Wendepunkt-Beschriftung korrigiert
%  11.01.2021 JMH        Plot 2 aufgehübscht
%  12.01.2021 JMH        berechneWendepunkt repariert: dy_w >= 100
%  13.01.2021 JMH        7-Tage-Inzidenz ergänzt
%  16.01.2021 JMH        Hub-Schwelle für Wendepunkt-Berechnung abgesenkt: 0.26
%  17.01.2021 JMH        Quellen-Angabe (Robert Koch-Institut) ergänzt
%  20.01.2021 JMH        "Datum_Ende" durch "Halbwertszeit" ersetzt
%  24.01.2021 JMH        Faktor für 7-Tage-Inzidenz ergänzt
%  24.01.2021 JMH        "Quelle" => "Rohdaten-Quelle" in Bildüberschriften
%  24.01.2021 CN         Identische x-Achsen bei allen Plots
%  25.01.2021 JMH        7-Tage-Inzidenz geplottet (oberer Subplot von Plot 3)
%  26.01.2021 CN         "octave_label" repariert
%  26.01.2021 JMH        Erster_sinnvoller_Tag = 19.02.2020: Plot 3 aufgehübscht
%  27.01.2021 JMH        Plot 3 noch mehr aufgehübscht
%  27.01.2021 JMH        Berechnung von "Zunahme_krank_7d" korrigiert
%  28.01.2021 JMH        "Verdoppelungszeit" => "Verdopplungszeit"
%  30.01.2021 JMH        Plot 1: Wendepunkt-Beschriftungen korrigiert
%  05.02.2021 JMH        Aktuelle Sterberate korrigiert
%  11.02.2021 JMH        Hub-Schwelle für Wendepunkt-Berechnung abgesenkt: 0.24
%  13.02.2021 JMH        Hub-Schwelle für Wendepunkt-Berechnung abgesenkt: 0.23
%  15.02.2021 CN         "add_print_menu" ergänzt
%  15.02.2021 CN         Plot 3 aufgehübscht
%  15.02.2021 JMH        Aktuelle Sterberate korrigiert
%  15.02.2021 JMH        Plot 2 aufgehübscht
%  17.02.2021 JMH        Hub-Schwelle für Wendepunkt-Berechnung abgesenkt: 0.22
%  02.03.2021 JMH        Hub-Schwelle für Wendepunkt-Berechnung abgesenkt: 0.21
%  04.03.2021 JMH        Dunkelziffer der Infizierten ergänzt
%  13.03.2021 JMH        Plot 1 aufgehübscht
%  14.03.2021 JMH        Plot 1 korrigiert
%  17.03.2021 CN         Plot 1 Verzweigung (Fallzahl) verschlankt
%  17.03.2021 JMH        Aufgeräumt
%  03.04.2021 CN         Lockdown-Daten eingebaut
%  04.04.2021 CN         Robustheit bei Aufruf in Ausschnitten
%  05.04.2021 CN         "scale_axis" korrigiert
%  05.04.2021 JMH        Plot 3 aufgehübscht
%  05.04.2021 CN         "bestlegpos": automatische 'Location' für "legend"
%  08.04.2021 JMH        Aufgeräumt
%  08.04.2021 JMH        Hub-Schwelle für Wendepunkt-Berechnung abgesenkt: 0.19
%  08.04.2021 CN         Unterschiedliche Plot-Größen für Matlab und Octave
%  08.04.2021 JMH        Aufgeräumt
%  08.04.2021 CN         "bestlegpos" optimiert
%  08.04.2021 JMH        1. Impfung gegen Änderitis
%  11.04.2021 JMH        Fehler in Plot 3 repariert
%  24.04.2021 JMH        Ostern und Bundes-Notbremse ergänzt
%  26.04.2021 CN         Text-Clipping verbessert
%  28.04.2021 JMH        Aufgeräumt
%  15.05.2021 CN         Bilder aufgehübscht
%  17.05.2021 JMH        Aufgeräumt
%  09.06.2021 JMH        Positionen einiger Lockdown-Daten korrigiert
%  13.06.2021 CN         Schwelle_AnzKW: 6 => 10
%  13.06.2021 CN         "thin_labels" korrigiert
%  14.06.2021 JMH        Aufgeräumt
%  14.06.2021 JMH        Nochmal aufgeräumt
%  10.08.2021 CN         Verzweigung R-Parameter "gefilterter Eingang oder
%                        nicht" entfernt, "filtereRohdaten" modifiziert
%  11.08.2021 JMH        Aufgeräumt
%  14.09.2021 JMH        Lockdown-Daten aktualisiert
%  13.10.2021 JMH        Lockdown-Daten aktualisiert
%  25.10.2021 CN         Aufgeräumt
%  26.10.2021 JMH        Aufgeräumt
%  14.11.2021 JMH        Lockdown-Daten aktualisiert
%  20.11.2021 CN         "skaliere_Wochen" korrigiert
%  21.11.2021 JMH        Aufgeräumt
%  02.12.2021 JMH        Einwohnerzahl aktualisiert
%  04.12.2021 JMH        x-Achse: Kalenderwochen in 2020, 2021 und/oder 2022
%  19.12.2021 JMH        Lockdown-Daten aktualisiert
%  31.01.2022 JMH        Textausgabe repariert
%  02.02.2022 JMH        Lockdown-Daten repariert
%  20.02.2022 CN         "bestimme_Wendepunkte" korrigiert
%  20.02.2022 JMH        Aufgeräumt
% ------------------------------------------------------------------------------

%
% Faktor für 7-Tage-Inzidenz bestimmen
%
Einwohner      = 83129285;
Inzidenzfaktor = 7 * 100000 / Einwohner;

%
% Die Zeitachse ist in Tagen skaliert: 01.01.2020 = Tag Nr. 1
%
Tag_0 = datenum (2019, 12, 31);

%
% Lockdown-Daten von
% https://de.wikipedia.org/wiki/COVID-19-Pandemie_in_Deutschland
%
Text_lock  = {'1. Lockdown',           ...
              'Lockerungen',           ...
              'Lockdown light',        ...
              '2. Lockdown',           ...
              'Bundesnotbremse',       ...
              'Ende Bundesnotbremse',  ...
              '3G-Regel',              ...
              'Tests kostenpflichtig', ...
              'Tests kostenlos',       ...
              '3G (Arbeit und OEPNV)'};

Lockdown   = [2020,  3, 22, 2;   ... % 22.03.2020 1. Lockdown
              2020,  5,  6, 1;   ... % 06.05.2020 Lockerungen
              2020, 11,  2, 1.5; ... % 02.11.2020 Lockdown light
              2020, 12, 16, 0.5; ... % 16.12.2020 2. Lockdown
              2021,  4, 23, 5.5; ... % 23.04.2021 Bundesnotbremse
              2021,  6, 30, 4.5; ... % 30.06.2021 Ende Bundesnotbremse
              2021,  8, 23, 1.5; ... % 23.08.2021 3G-Regel
              2021, 10, 11, 0.5; ... % 11.10.2021 Tests kostenpflichtig
              2021, 11, 13, 1.5; ... % 13.11.2021 Tests kostenlos
              2021, 11, 24, 2.5];    % 24.11.2021 3G (Arbeit und ÖPNV)

deltaHoehe = 0.15 * Lockdown (:, 4);
nLock      = numel (deltaHoehe);

t_lock     = (berechneZeit (Lockdown (:, 1), Lockdown (:, 2), ...
                            Lockdown (:, 3), 0, 0, Tag_0)).';

Datum_lock = sprintf ('%02d.%02d.%4d', (Lockdown (:, [3, 2, 1])).');
Datum_lock = (reshape (Datum_lock, 10, (numel (Datum_lock) / 10))).';

%
% Rückgabewert initialisieren
%
Text = [];

%
% Pfad bestimmen
%
if exist ('/home/jmh', 'dir')
  Pfad = '/home/jmh/octave/corona/';
else
  Pfad = '';
end

%
% Benötigte Toolbox laden
%
persistent pkg_loaded

if isempty (pkg_loaded)
  if ~isempty (strfind (license, 'GNU'))
    pkg load signal
  end
  pkg_loaded = 1;
end

%
% ASCII-Datei mit Daten einlesen
%
if nargin < 3
  CoronaDatei = 'CoronaDatenRKI.txt';
end

[~, Rumpfname] = fileparts (CoronaDatei);
CoronaDatei    = [Rumpfname, '.txt'];

Daten          = load ([Pfad, CoronaDatei]);

%
% Bei Bedarf Daten beschneiden
%
if nargin > 0

  if nargin < 2
    Index_max = size (Daten, 1);
  end

  Index_max = min (Index_max, size (Daten, 1));

  if (0 < Index_min) && (Index_min < Index_max)

    %
    % Daten nur von Index_min bis Index_max auswerten
    %
    Daten = Daten (Index_min:Index_max, :);

  else

    %
    % Gefilterte, äquidistant gesampelte ASCII-Datei schreiben:
    % Zeitinkrement wählen (1/12 => 12x in 24h, 0.11 gewählt,
    % damit die Zeiten etwas "abwechslungsreicher" erscheinen.)
    %
    ts     = 0.11;

    Zeit_d = (berechneZeit (Daten (:, 1), Daten (:, 2), Daten (:, 3), ...
                            Daten (:, 4), Daten (:, 5), Tag_0)).';

    % Filterung durchführen
    [Zeit_d_sn, Kranke_fn]     = filtereRohdaten (Zeit_d, Daten (:,  6), ts);
    [~,         Gesund_fn]     = filtereRohdaten (Zeit_d, Daten (:,  7), ts);
    [~,         Tote_fn]       = filtereRohdaten (Zeit_d, Daten (:,  8), ts);
    [~,         Kranke_MD_fn]  = filtereRohdaten (Zeit_d, Daten (:,  9), ts);
    [~,         Kranke_oNC_fn] = filtereRohdaten (Zeit_d, Daten (:, 10), ts);

    Daten_flt       = [berechneDatum(Zeit_d_sn, Tag_0), ...
                       Kranke_fn,                       ...
                       Gesund_fn,                       ...
                       Tote_fn,                         ...
                       Kranke_MD_fn,                    ...
                       Kranke_oNC_fn];

    % Letzten Datensatz am Ende der gefilterten Daten wieder dranhängen
    Daten_flt       = [Daten_flt; Daten(end, :)];

    CoronaDatei_flt = [Rumpfname, '_flt.txt'];
    Datei           = fopen ([Pfad, CoronaDatei_flt], 'wt');

    Format = '%4d %02d %02d %02d %02d %11.2f %11.2f %11.2f %11.2f %11.2f\n';

    for Index = 1:size (Daten_flt, 1)
      fprintf (Datei, Format, Daten_flt (Index, :));
    end

    fclose (Datei);

    return;

  end

end

%
% Daten für die Plots extrahieren
%
Zeit_d     = berechneZeit (Daten (:, 1), Daten (:, 2), Daten (:, 3), ...
                           Daten (:, 4), Daten (:, 5), Tag_0);

Kranke_MD  = (Daten (:,  9)).';
Kranke_oNC = (Daten (:, 10)).';
Kranke     = (Daten (:,  6)).';
Gesund     = (Daten (:,  7)).';
Tote       = (Daten (:,  8)).';

if nargin < 4
  Modus = 3;
end

if     Modus == 1
  Kranke = Kranke_MD;
elseif Modus == 2
  Kranke = Kranke_oNC;
end

Momentan   = Kranke - (Gesund + Tote);

%
% Textausgabe erstellen
%
Stand = sprintf ('Stand %02d.%02d.%4d %02d:%02d', Daten (end, [3, 2, 1, 4, 5]));

if     Modus == 1
  Stand     = [Stand, ', Auswertung mit Meldedatum'];
  StandText = Stand;
elseif Modus == 2
  Stand     = [Stand, ', Erkrankungsdatum ohne Nowcasting'];
  StandText = Stand;
else
  StandText = [Stand, ', Quelle: Robert Koch-Institut'];
  Stand     = [Stand, ', Rohdaten-Quelle: Robert Koch-Institut'];
end

KrankeEnd   = round (Kranke   (end));
MomentanEnd = round (Momentan (end));
ToteEnd     = round (Tote     (end));

Text = sprintf ('COVID-19-Pandemie, Ausbreitung in Deutschland, %s', StandText);
Text = char (Text, sprintf ('Summe Infizierte    = %8d Personen', KrankeEnd));
Text = char (Text, sprintf ('Momentan Kranke     = %8d Personen', MomentanEnd));
Text = char (Text, sprintf ('Gestorbene          = %8d Personen', ToteEnd));

%
% Wendepunkte berechnen
%
[w_krank,    yw_krank, dyw_krank] = berechneWendepunkt (Zeit_d, Kranke,   0);
[w_momentan, yw_momentan]         = berechneWendepunkt (Zeit_d, Momentan, 1);
[w_tot,      yw_tot,   dyw_tot]   = berechneWendepunkt (Zeit_d, Tote,     0);

w_Datum_krank       = berechneDatum (w_krank,    Tag_0);
w_Datum_momentan    = berechneDatum (w_momentan, Tag_0);
w_Datum_tot         = berechneDatum (w_tot,      Tag_0);

idx                 = [3, 2, 1];
Wendedatum_krank    = sprintf ('%02d.%02d.%4d', (w_Datum_krank    (:, idx)).');
Wendedatum_momentan = sprintf ('%02d.%02d.%4d', (w_Datum_momentan (:, idx)).');
Wendedatum_tot      = sprintf ('%02d.%02d.%4d', (w_Datum_tot      (:, idx)).');

if numel (w_krank)    > 1
  Anzahl              = numel (w_krank);
  Wendedatum_krank    = (reshape (Wendedatum_krank,    10, Anzahl)).';
end
if numel (w_momentan) > 1
  Anzahl              = numel (w_momentan);
  Wendedatum_momentan = (reshape (Wendedatum_momentan, 10, Anzahl)).';
end
if numel (w_tot)      > 1
  Anzahl              = numel (w_tot);
  Wendedatum_tot      = (reshape (Wendedatum_tot,      10, Anzahl)).';
end

%
% Runden auf ganzzahlige Werte
%
Kranke   = round (Kranke);
Momentan = round (Momentan);
Tote     = round (Tote);

%
% Tägliche Neuinfektionen und Tote bestimmen
%
if std (diff (Zeit_d)) < 0.005
  Daten_sind_aequidistant = true;
else
  Daten_sind_aequidistant = false;
end

Tag_min       = ceil  (Zeit_d (1) + 1 / (24 * 60));
Tag_max       = floor (Zeit_d (end));

Tagesstempel  = zeros (1, (Tag_max - Tag_min + 1));
Zunahme_krank = Tagesstempel;
Zunahme_tot   = Tagesstempel;

Index         = 0;

for Tag = Tag_min:Tag_max

  Index = Index + 1;

  if Daten_sind_aequidistant
    Index_gestern = find ((Zeit_d <  Tag),        1, 'last');
    Index_heute   = find ((Zeit_d < (Tag + 1)),   1, 'last');
  else
    Index_gestern = find ((Zeit_d < (Tag + 0.2)), 1, 'last');
    Index_heute   = find ((Zeit_d < (Tag + 1.2)), 1, 'last');
  end

  Tagesstempel  (Index) = Tag;
  Zunahme_krank (Index) = Kranke (Index_heute) - Kranke (Index_gestern);
  Zunahme_tot   (Index) = Tote   (Index_heute) - Tote   (Index_gestern);

end

Tage_krank_7d    = zeros (1, (numel (Zunahme_krank) - 6));
Zunahme_krank_7d = Tage_krank_7d;

for Index = 1:numel (Zunahme_krank_7d)
  Tage_krank_7d    (Index) = Tagesstempel (Index + 3);
  Zunahme_krank_7d (Index) = mean (Zunahme_krank (Index:(Index + 6)));
end

if isempty (Zunahme_krank_7d)
  Tage_krank_7d    = NaN;
  Zunahme_krank_7d = NaN;
end

Tage_tot_7d      = zeros (1, (numel (Zunahme_tot)   - 7));
Zunahme_tot_7d   = Tage_tot_7d;

for Index = 1:numel (Zunahme_tot_7d)
  Tage_tot_7d      (Index) = Tagesstempel (Index + 3);
  Zunahme_tot_7d   (Index) = mean (Zunahme_tot   (Index:(Index + 6)));
end

if isempty (Zunahme_tot_7d)
  Tage_tot_7d      = NaN;
  Zunahme_tot_7d   = NaN;
end

if ~isnan (Zunahme_krank_7d (end))

  Unit = 'Infizierte/Woche/100.000 Einwohner';
  Text = char (Text, sprintf ('7-Tage-Inzidenz     = %10.1f %s', ...
                              (Zunahme_krank_7d (end) * Inzidenzfaktor), Unit));

  Text = char (Text, sprintf ('7-Tage-Mittelwert K = %10.1f Infizierte/Tag', ...
                              Zunahme_krank_7d (end)));

end

if ~isnan (Zunahme_tot_7d   (end))
  Text = char (Text, sprintf ('7-Tage-Mittelwert T = %10.1f Gestorbene/Tag', ...
                              Zunahme_tot_7d   (end)));
end

Zaehler = Zunahme_tot_7d (end);

%
% Die Maxima von Zunahme_tot_7d und Zunahme_krank_7d sind um 14 Tage verschoben.
% Außerdem ist der Vektor Zunahme_krank_7d um einen Tag länger als der Vektor
% Zunahme_tot_7d: Daher muß man bei Zunahme_krank_7d um 15 Tage zurückgehen.
%
if numel (Zunahme_krank_7d) > 15
  Nenner = Zunahme_krank_7d (end - 15);
else
  Nenner = Zunahme_krank_7d (1);
end

if Nenner > 0
  Sterberate = 100 * Zaehler / Nenner;
else
  Sterberate = NaN;
end

if ~isnan (Sterberate)
  Text = char (Text, sprintf ('Aktuelle Sterberate = %11.2f%%', Sterberate));
end

Sterberate_ref = 0.8;

if Sterberate > Sterberate_ref
  Dunkelziffer = round (100 * (Sterberate / Sterberate_ref - 1));
else
  Dunkelziffer = NaN;
end

if ~isnan (Dunkelziffer)
  Unit = 'der gemeldeten Infizierten zusätzlich';
  Unit = sprintf ('%s, wenn Sterberate_ref = %4.2f%%', Unit, Sterberate_ref);
  Text = char (Text, sprintf ('Dunkelziffer        = %8d%% %s', ...
               Dunkelziffer, Unit));
end

KeineZunahme     = find (Zunahme_krank_7d <= 0);
Tage_krank_7d    (KeineZunahme) = [];
Zunahme_krank_7d (KeineZunahme) = [];
Tage_krank_7d    = [NaN, Tage_krank_7d];
Zunahme_krank_7d = [NaN, Zunahme_krank_7d];

KeineZunahme     = find (Zunahme_tot_7d   <= 0);
Tage_tot_7d      (KeineZunahme) = [];
Zunahme_tot_7d   (KeineZunahme) = [];
Tage_tot_7d      = [NaN, Tage_tot_7d];
Zunahme_tot_7d   = [NaN, Zunahme_tot_7d];

Zunahme_krank_orig = Zunahme_krank;

if Zunahme_krank (end) < 1
  Zunahme_krank (end) = NaN;
end

KeineZunahme     = find (Zunahme_krank    <  1);
Tage_krank       = Tagesstempel;
Tage_krank       (KeineZunahme) = [];
Zunahme_krank    (KeineZunahme) = [];
Tage_krank       = [NaN, Tage_krank];
Zunahme_krank    = [NaN, Zunahme_krank];

if Zunahme_tot   (end) < 1
  Zunahme_tot   (end) = NaN;
end

KeineZunahme     = find (Zunahme_tot      <  1);
Tage_tot         = Tagesstempel;
Tage_tot         (KeineZunahme) = [];
Zunahme_tot      (KeineZunahme) = [];
Tage_tot         = [NaN, Tage_tot];
Zunahme_tot      = [NaN, Zunahme_tot];

%
% Reproduktionszahl R aus den Messwerten berechnen
%
Nsumme = 7;
Start  = 6; % Von den 11 Messwerten ist der Wert Nr. 6 der mittlere (Referenz)
Offset = 5;

meanWoche = zeros (1, (numel (Zunahme_krank_orig) - (Nsumme - 1)));
stdWoche  = meanWoche;

for Index = 1:numel (meanWoche)
  meanWoche (Index) = mean (Zunahme_krank_orig (Index:(Index + (Nsumme - 1))));
  stdWoche  (Index) = std  (Zunahme_krank_orig (Index:(Index + (Nsumme - 1))));
end

R    = NaN * ones (1, (numel (meanWoche) - 4));
stdR = R;

for Index = 1:numel (R)
  if meanWoche (Index) > 0

    R    (Index) = meanWoche (Index + 4) / meanWoche (Index);

    % stdR^2 = (dR/dZ * stdZ)^2 + (dR/dN * stdN)^2
    stdR (Index) = sqrt ((stdWoche (Index + 4) / meanWoche (Index))^2 + ...
                         (stdWoche (Index) * meanWoche (Index + 4) / ...
                          (meanWoche (Index))^2)^2);

  end
end

Erster_sinnvoller_Tag = 50; % 19.02.2020

if (berechneZeit (Daten (end, 1), Daten (end, 2), Daten (end, 3), ...
                  Daten (end, 4), Daten (end, 5), Tag_0) < ...
    (Erster_sinnvoller_Tag + Offset)) || isempty (R)

  R_ist_sinnvoll = 0;

  Tagesstempel   = Tagesstempel (1);
  R              = NaN;
  stdR           = NaN;

else

  R_ist_sinnvoll = 1;

  Tagesstempel   = Tagesstempel (Start:(end - Offset));
  Sinnvoll       = find (Tagesstempel >= Erster_sinnvoller_Tag);

  Tagesstempel   = Tagesstempel (Sinnvoll);
  R              = R            (Sinnvoll);
  stdR           = stdR         (Sinnvoll);

end

if (R_ist_sinnvoll > 0) && (numel (R) > 1)

  if     R (end - 1) < 1
%    Dauer            = -round (4 * log (Momentan (end)) / log (R (end - 1)));
%    Enddatum         = berechneDatum ((floor (Zeit_d (end)) + Dauer), Tag_0);
%    Datum_Ende       = sprintf ('%02d.%02d.%4d', Enddatum ([3, 2, 1]));
    Halbwertszeit    = -4 * log (2) / log (R (end - 1));
  elseif R (end - 1) > 1
    Verdopplungszeit =  4 * log (2) / log (R (end - 1));
  end

  Text = char (Text, ...
               sprintf ('Reproduktionszahl R = %11.2f+/-%.2f (%.2f+/-%.2f)', ...
                        R (end - 1), stdR (end - 1), R (end), stdR (end)));

%  if     exist ('Datum_Ende',       'var')
%    Text = char (Text, sprintf ('Pandemie-Ende aus R = %s', Datum_Ende));
  if     exist ('Halbwertszeit',    'var')
    Text = char (Text, sprintf ('Halbwertszeit       = %10.1f Tage', ...
                                Halbwertszeit));
  elseif exist ('Verdopplungszeit', 'var')
    Text = char (Text, sprintf ('Verdopplungszeit    = %10.1f Tage', ...
                                Verdopplungszeit));
  end

end

%
% Plot-Farben definieren
%
dunkelrot   = [0.5, 0,   0];
dunkelgruen = [0,   0.4, 0];

%
% 1. Plot: Verlauf der Infektionen mit dem Corona-Virus in Deutschland
%
figure (1);
clf;

if exist ('/home/jmh', 'dir')
  set (1, 'Position',[7, 664, 560, 420]);
else
  cn_size = [640, 480];
  if isempty (strfind (license, 'GNU'))
    cn_size = [800, 600];
  end
  set (1, 'Position',[5, 55, cn_size]);
  add_print_menu (1);
end

ax1  = subplot (211);
ax2  = subplot (212);
pos1 = get (ax1, 'Position');
pos2 = get (ax2, 'Position');
dy1  = 0.14;
dy2  = 0.02;
set (ax1, 'Position',(pos1 + [0, -dy1, 0,  dy1]));
set (ax2, 'Position',(pos2 + [0,  0,   0, -dy2]));

Y_Label = 'Anzahl Personen';
Divisor = 1;
if max (Kranke) > 7e5
  Y_Label = 'Anzahl [1000 Personen]';
  Divisor = 1000;
end

%
% Horizontaler Abstand der Beschriftungen von den senkrechten Linien
%
Zwei    = ones (1, 2);
Abstand = (Zeit_d (end) - Zeit_d (1)) / 7 / 250;

x_min   = KW (floor (Zeit_d (1)));
x_max   = KW (ceil  (Zeit_d (end)));
y_min   = 0;
y_max1  = ceil (max ((Kranke (end) / Divisor), 1) / 100) * 100;
y_max2  = ceil (max ((Tote   (end) / Divisor), 1) /  10) *  10;

axes (ax1);
hold on;
hp2 = plot (KW (Zeit_d), (Momentan / Divisor), '-b', 'LineWidth',2);
hp1 = plot (KW (Zeit_d), (Kranke   / Divisor), '-r', 'LineWidth',2);
ht1 = [];
if any (~isnan (w_momentan))
  Null = zeros (size (w_momentan));
  plot   ((KW (w_momentan) * Zwei).',  [Null, (yw_momentan / Divisor)].', ...
          '--b');
  ht1  = ...
    text ((KW (w_momentan) + Abstand), (2.0 * (yw_momentan / Divisor)), ...
          Wendedatum_momentan, ...
          'FontSize',9, 'Color','b',       'Clipping','on');
end
ht2 = [];
if any (~isnan (w_krank))
  Null = zeros (size (w_krank));
  plot   ((KW (w_krank)    * Zwei).',  [Null, (yw_krank    / Divisor)].', ...
          '--r');
  ht2  = ...
    text ((KW (w_krank)    + Abstand), (0.7 * (yw_krank    / Divisor)), ...
          Wendedatum_krank,    ...
          'FontSize',9, 'Color',dunkelrot, 'Clipping','on');
end
hold off;

axis ([x_min, x_max, y_min, y_max1]);
scale_axis (10);
skaliere_Wochen (gca, x_min, x_max);
set (ax1, 'xTickLabel',[]);
ht = [ht1; ht2];
limit_text_to_axis (ht (ishandle (ht)), ax1);
grid ('on');
box  ('on');

axes (ax2);
hold on;
hp3 = plot (KW (Zeit_d), (Tote     / Divisor), '-k', 'LineWidth',2);
hp4 = plot (KW (Zeit_d), (Momentan / Divisor), '-b', 'LineWidth',2);
ht1 = [];
if any (~isnan (w_tot))
  Null = zeros (size (w_tot));
  plot   ((KW (w_tot)      * Zwei).',  [Null, (yw_tot      / Divisor)].', ...
          '--k');
  ht1  = ...
    text ((KW (w_tot)      + Abstand), (0.8 * (yw_tot      / Divisor)), ...
          Wendedatum_tot,      ...
          'FontSize',9, 'Color','k',       'Clipping','on');
end
ht2 = [];
if any (~isnan (w_momentan))
  Null = zeros (size (w_momentan));
  plot   ((KW (w_momentan) * Zwei).',  [Null, (yw_momentan / Divisor)].', ...
          '--b');
  ht2  = ...
    text ((KW (w_momentan) + Abstand), (0.7 * (yw_momentan / Divisor)), ...
          Wendedatum_momentan, ...
          'FontSize',9, 'Color','b',       'Clipping','on');
end
hold off;

axis ([x_min, x_max, y_min, y_max2]);
scale_axis (7);
skaliere_Wochen (gca, x_min, x_max);
thin_labels     (gca, 'x', 14);
ht = [ht1; ht2];
limit_text_to_axis (ht (ishandle (ht)), ax2);
grid ('on');
box  ('on');

axes (ax1);
LegPos = bestlegpos ([hp1, hp2], ax1);
hl     = legend ([hp1, hp2], ...
                 {sprintf('Infizierte (akkumuliert): %d', Kranke(end)); ...
                  sprintf('momentan Kranke: %d',          Momentan(end))}, ...
                 'Location',LegPos);
stretch_legend (hl, 1.07, LegPos);

Titel = {'Verlauf der Infektionen mit SARS-CoV-2 in Deutschland'; ...
         sprintf('(%s)', Stand)};

title  (Titel);
adjust_title_position (gca, [0.5, 1.02, 0], 'GNU');
xlabel ('');
ylabel (Y_Label);

octave_label (ax1, 'y');

axes (ax2);
LegPos = bestlegpos (hp3, ax2);
hl     = legend ([hp3, hp4], ...
                 {sprintf('an COVID-19 Gestorbene: %d',   Tote(end));   ...
                  sprintf('momentan Kranke: %d',          Momentan(end))}, ...
                 'Location',LegPos);
stretch_legend (hl, 1.07, LegPos);

%xlabel ('Zeit [Kalenderwoche im Jahr 2020]');
ylabel (Y_Label);

octave_label   (ax2, 'y');
cut_last_label (ax2, 'y');

drawnow;

%
% 2. Plot: Tägliche Neuinfektionen und Tote
%
figure (2);
clf;

if exist ('/home/jmh', 'dir')
  set (2, 'Position',[521, 664, 560, 420]);
else
  set (2, 'Position',[230, 55, cn_size]);
  add_print_menu (2);
end

x_min = KW (floor (Zeit_d (1)));
x_max = KW (ceil  (Zeit_d (end)));
y_max = 10^ceil (log10 (max ([Zunahme_krank(:); Zunahme_tot(:); 10])));

hold off;
hp1 = semilogy (KW (Tage_krank (1:(end - 1))), Zunahme_krank (1:(end - 1)), ...
                '*', 'Color',dunkelrot,   'LineWidth',1);
hold on;
hp2 = semilogy (KW (Tage_tot   (1:(end - 1))), Zunahme_tot   (1:(end - 1)), ...
                '*', 'Color',dunkelgruen, 'LineWidth',1);
hp3 = semilogy (KW (Tage_krank_7d),            Zunahme_krank_7d,            ...
                '-', 'Color','r',         'LineWidth',2);
hp4 = semilogy (KW (Tage_tot_7d),              Zunahme_tot_7d,              ...
                '-', 'Color','k',         'LineWidth',2);
semilogy (KW (Tage_krank (end)), Zunahme_krank (end), '*r', 'LineWidth',1);
semilogy (KW (Tage_tot   (end)), Zunahme_tot   (end), '*k', 'LineWidth',1);
ht1 = [];
if any (~isnan (w_tot))
  Eins = ones (size (w_tot));
  semilogy    ((KW (w_tot)   * Zwei).',  [Eins, dyw_tot].',   '--k');
  ht1  = text ((KW (w_tot)   + Abstand), (0.3 * dyw_tot),   ...
               Wendedatum_tot,   ...
               'FontSize',9, 'Color','k',       'Clipping','on');
end
ht2 = [];
if any (~isnan (w_krank))
  Eins = ones (size (w_krank));
  semilogy    ((KW (w_krank) * Zwei).',  [Eins, dyw_krank].', '--r');
  ht2  = text ((KW (w_krank) + Abstand), (0.1 * dyw_krank), ...
               Wendedatum_krank, ...
               'FontSize',9, 'Color',dunkelrot, 'Clipping','on');
end
hold off;

axis ([x_min, x_max, 1, 10^log10(y_max)]);
set  (gca, 'yTick',10.^(0:log10 (y_max)));
skaliere_Wochen (gca, x_min, x_max);
thin_labels     (gca, 'x', 14);
ht = [ht1; ht2];
limit_text_to_axis (ht (ishandle (ht)), gca);
grid ('on');
box  ('on');

LegPos = bestlegpos ([hp1, hp2], gca);
hl     = legend ([hp1, hp2, hp3, hp4], ...
                 {'Infizierte/Tag';                     ...
                  'Gestorbene/Tag';                     ...
                  '7-Tage-Mittelwert "Infizierte/Tag"'; ...
                  '7-Tage-Mittelwert "Gestorbene/Tag"'}, 'Location',LegPos);
stretch_legend (hl, 1.07, LegPos);

Titel = {'COVID-19: Infizierte/Tag und Gestorbene/Tag in Deutschland'; ...
         sprintf('(%s)', Stand)};

title  (Titel);
adjust_title_position (gca, [0.5, 1.01, 0], 'GNU');
%xlabel ('Zeit [Kalenderwoche im Jahr 2020]');
ylabel ('Anzahl Personen in Deutschland');

drawnow;

%
% 3. Plot: 7-Tage-Inzidenz und Reproduktionszahl R
%
figure (3);
clf;

if exist ('/home/jmh', 'dir')
  set (3, 'Position',[1035, 664, 560, 420]);
else
  set (3, 'Position',[480, 55, cn_size]);
  add_print_menu (3);
end

ax1  = subplot (211);
ax2  = subplot (212);
pos1 = get (ax1, 'Position');
pos2 = get (ax2, 'Position');
dy   = 0.06;
set (ax1, 'Position',(pos1 + [0, -dy, 0, dy]));
set (ax2, 'Position',(pos2 + [0,  0,  0, dy]));

Inzidenz    = Inzidenzfaktor * Zunahme_krank_7d;
yw_Inzidenz = Inzidenzfaktor * dyw_krank;

x_min       = KW (floor (Zeit_d (1)));
x_max       = KW (ceil  (Zeit_d (end)));
y_min       = 0;
y_max1      = ceil (max ([Inzidenz, 1]) / 10) * 10;
y_max2      = 2.1;

axes (ax1);
hold on;
hp = plot (KW (Tage_krank_7d), Inzidenz, '-', 'Color',dunkelrot, 'LineWidth',2);
Eins = ones (size (nLock));
plot        ((KW (t_lock) * Zwei).',  (Eins * [0, y_max1]).', '-k');
ht1  = text ((KW (t_lock) + Abstand), ((0.05 + deltaHoehe) * y_max1), ...
             Datum_lock, 'FontSize',9, 'Color','k', 'Clipping','on');
ht2  = text ((KW (t_lock) + Abstand),         (deltaHoehe  * y_max1), ...
             Text_lock,  'FontSize',9, 'Color','k', 'Clipping','on');
if any (~isnan (w_krank))
  Null = zeros (size (w_krank));
  plot ((KW (w_krank) * Zwei).',  [Null, yw_Inzidenz].', '--r');
end
hold off;

axis ([x_min, x_max, y_min, y_max1]);
scale_axis (10);
skaliere_Wochen (gca, x_min, x_max);
set (ax1, 'xTickLabel',[]);
ht = [ht1; ht2];
limit_text_to_axis (ht (ishandle (ht)), ax1);
grid ('on');
box  ('on');

axes (ax2);
Ro = R + stdR;
Ru = R - stdR;
hold on;
hp1 = plot ([x_min, x_max],    [1, 1], '-g', 'LineWidth',2);
hp3 = plot (KW (Tagesstempel), Ro,     '-b', 'LineWidth',1);
hp4 = plot (KW (Tagesstempel), Ru,     '-b', 'LineWidth',1);
hp2 = plot (KW (Tagesstempel), R,      '-r', 'LineWidth',2);
ht  = [];
if any (~isnan (w_krank))
  Eins = ones (size (w_krank));
  plot        ((KW (w_krank) * Zwei).',  (Eins * [0, y_max]).', '--r');
  ht   = text ((KW (w_krank) + Abstand), (Eins * 0.25), Wendedatum_krank, ...
               'FontSize',9, 'Color',dunkelrot, 'Clipping','on');
end
hold off;

axis ([x_min, x_max, y_min, y_max2]);
scale_axis (10);
skaliere_Wochen (gca, x_min, x_max);
thin_labels     (gca, 'x', 14);
limit_text_to_axis (ht (ishandle (ht)), ax2);
grid ('on');
box  ('on');

axes (ax1);
LegPos = bestlegpos (hp, ax1);
hl     = legend (hp, ...
                 {['7-Tage-Mittelwert', char(10), ...
                   '"Infizierte/Woche/100.000 Einwohner"']}, ...
                 'Location',LegPos);
stretch_legend (hl, 1.07, LegPos);

%if     exist ('Datum_Ende',       'var')
%  Titel = {['SARS-CoV-2: 7-Tage-Inzidenz und Reproduktionszahl R', ...
%            ' (Generationsdauer = 4 Tage) in Deutschland']; ...
%           sprintf('(Ende der Pandemie = %s, %s)', Datum_Ende, Stand)};
if     exist ('Halbwertszeit',    'var')
  Titel = {['SARS-CoV-2: 7-Tage-Inzidenz und Reproduktionszahl R', ...
            ' (Generationsdauer = 4 Tage) in Deutschland']; ...
           sprintf('(Halbwertszeit = %.1f Tage, %s)', Halbwertszeit, Stand)};
elseif exist ('Verdopplungszeit', 'var')
  Titel = {['SARS-CoV-2: 7-Tage-Inzidenz und Reproduktionszahl R', ...
            ' (Generationsdauer = 4 Tage) in Deutschland']; ...
           sprintf('(Verdopplungszeit = %.1f Tage, %s)', ...
                   Verdopplungszeit, Stand)};
else
  Titel = {['SARS-CoV-2: 7-Tage-Inzidenz und Reproduktionszahl R', ...
            ' (Generationsdauer = 4 Tage) in Deutschland']; ...
           sprintf('(%s)', Stand)};
end

title  (Titel);
adjust_title_position (gca, [0.5, 1.02, 0], 'GNU');
xlabel ('');
ylabel ('Inzidenz-Wert');

octave_label (ax1, 'y');

axes (ax2);
LegPos = bestlegpos (hp2, ax2);
if isnan (R)
  hl = legend (hp1, 'R = 1 => lineares Wachstum', 'Location',LegPos);
else
  hl = legend ([hp1, hp2, hp3], ...
               {'R = 1 => lineares Wachstum'; ...
                'Reproduktionszahl R';        ...
                'Standardabweichung von R'}, 'Location',LegPos);
end
stretch_legend (hl, 1.07, LegPos);

%xlabel ('Zeit [Kalenderwoche im Jahr 2020]');
ylabel ('Reproduktionszahl R');

octave_label   (ax2, 'y');
cut_last_label (ax2, 'y');

drawnow;

%
% MAT-Datei mit akkumulierten Ergebnissen von corona3.m
%
if exist ('/home/jmh', 'dir') && (nargin < 1)

  Pfad      = '/bckdsk/tmp/RKI-Download/mat/';
  Dateiname = [Pfad, 'CoronaDatenRKI.mat'];

  if exist (Dateiname, 'file') == 2
    load (Dateiname);
  end

  IndexTag = floor (max (Zeit_d)) - 232; % 20.08.2020 = Tag Nr. 1

  CoronaDatenRKI (IndexTag).Text             = Text;

  CoronaDatenRKI (IndexTag).Zeit_d           = Zeit_d;
  CoronaDatenRKI (IndexTag).Kranke_MD        = Kranke_MD;
  CoronaDatenRKI (IndexTag).Kranke_oNC       = Kranke_oNC;
  CoronaDatenRKI (IndexTag).Kranke           = Kranke;
  CoronaDatenRKI (IndexTag).Momentan         = Momentan;
  CoronaDatenRKI (IndexTag).Tote             = Tote;

  CoronaDatenRKI (IndexTag).Tage_krank       = Tage_krank;
  CoronaDatenRKI (IndexTag).Zunahme_krank    = Zunahme_krank;
  CoronaDatenRKI (IndexTag).Tage_tot         = Tage_tot;
  CoronaDatenRKI (IndexTag).Zunahme_tot      = Zunahme_tot;

  CoronaDatenRKI (IndexTag).Tage_krank_7d    = Tage_krank_7d;
  CoronaDatenRKI (IndexTag).Zunahme_krank_7d = Zunahme_krank_7d;
  CoronaDatenRKI (IndexTag).Tage_tot_7d      = Tage_tot_7d;
  CoronaDatenRKI (IndexTag).Zunahme_tot_7d   = Zunahme_tot_7d;

  CoronaDatenRKI (IndexTag).Tagesstempel     = Tagesstempel;
  CoronaDatenRKI (IndexTag).R                = R;
  CoronaDatenRKI (IndexTag).stdR             = stdR;

  save (Dateiname, 'CoronaDatenRKI');

end

end

function Zeit_d = berechneZeit (Jahr, Monat, Tag, Stunde, Minute, Tag_0)

Zeit_d = (datenum (Jahr, Monat, Tag) + (Stunde + Minute / 60) / 24 - Tag_0).';

end

function DatumUhrzeit = berechneDatum (Zeit_d, Tag_0)

DatumUhrzeit = zeros (numel (Zeit_d), 5);

for Index = 1:numel (Zeit_d)

  % Runden auf ganze Minuten:
  Vektor = datevec (round (24 * 60 * (Zeit_d (Index) + Tag_0)) / (24 * 60));

  DatumUhrzeit (Index, :) = Vektor (1:5);

end

end

function Kalenderwoche = KW (Tag)

Kalenderwoche = (Tag + 8) / 7;

end

function [xData_sn, yData_fn, yData_fo] = filtereRohdaten (xData, yData, ts)
% Das Unterprogramm macht:
% a) Resampling der Daten auf dem Wertebereich von "xData" mit Intervall 0,08
% b) Tiefpassfilterung der Daten
% c) Dezimation der gefilterten Daten sowie
%    Resampling entlang der Original-Abtastwerte von "xData"
% Eingang:
%   xData = x-Werte (Zeitwerte), i.a. nicht äquidistant
%   yData = y-Werte (Messwerte) dazu
%   ts    = Sampling-Intervall in der Skalierung von "xData"
% Ausgang:
%   xData_sn = x-Werte im Wertebereich von "xData", äquidist. mit Inkrement "ts"
%   yData_fn = gefilterte "yData"-Werte, äquidistant mit "ts" abgetastet
%   yData_fo = gefilterte "yData"-Werte, abgetastet mit "xData"

if nargin < 3
  ts = 0.11;
end

% Filter-Koeffizienten für Tiefpass, erzeugt mit Matlab "fdatool":
% Lowpass, FIR, Window, Order 64, Scale Passb., Kaiser (0.5), wc 0.01 (norm.)
b2      = [0.013130554374, 0.013328349337, 0.013521939100, 0.013711139006, ...
           0.013895767916, 0.014075648442, 0.014250607172, 0.014420474893, ...
           0.014585086814, 0.014744282772, 0.014897907442, 0.015045810539, ...
           0.015187847011, 0.015323877229, 0.015453767167, 0.015577388580, ...
           0.015694619166, 0.015805342730, 0.015909449334, 0.016006835441, ...
           0.016097404051, 0.016181064830, 0.016257734224, 0.016327335573, ...
           0.016389799210, 0.016445062553, 0.016493070186, 0.016533773934, ...
           0.016567132926, 0.016593113649, 0.016611689989, 0.016622843271, ...
           0.016626562278];

b       = [b2, b2((end - 1):-1:1)].';
nfilt   = numel (b);

% Work-Around um Kopie des letzten Wertes
del_x = diff (xData ((end - 2):end));
del_y = diff (yData ((end - 2):end));
if del_y (end) == 0
  del_y (end) = del_y (end - 1) * del_x (end) / del_x (end - 1);
  yData (end) = yData (end) + del_y (end);
end

% Resampling der Daten auf dem Wertebereich von "xData" mit Intervall 0,08
delta_x = 0.08; % 0.02 = kaum Filterwirkung, 0.1 = starke Filterwirkung
x_neu   = (xData (1):delta_x:xData (end)).';

if numel (x_neu) > 1
  y_neu = interp1 (xData, yData, x_neu, 'pchip');
else
  y_neu = yData (1);
end

if numel (y_neu) < (3 * (nfilt - 1) + 1)
  for Index = 1:(3 * (nfilt - 1) + 1 - numel (y_neu))
    y_neu = [y_neu; y_neu(end)];
  end
end

% Tiefpassfilterung der Daten
yf = filtfilt (b, 1, y_neu);
yf = yf (1:numel (x_neu));

if numel (x_neu) > 1

  % Dezimation der gefilterten Daten
  xData_sn = (x_neu (1):ts:(x_neu (end) - ts / 2)).';
  yData_fn = interp1 (x_neu, yf, xData_sn, 'linear');

  % Resampling entlang der Original-Abtastwerte von "xData"
  yData_fo = interp1 (x_neu, yf, xData,    'linear');

else

  xData_sn = xData (1);
  yData_fn = yData (1);
  yData_fo = yData (1);

end

end

function [x_w, y_w, dy_w] = berechneWendepunkt (xData, Messung, dFlag, Hub_min)
% Das Unterprogramm berechnet zu den Daten "Messung" die Wendepunkte.
% Voraussetzung: d/dt(Messung) hat konvexe Verlaufsabschnitte (=> suche Maxima)
% Spezialfall: Eingang ist schon Differenzkurve (=> dFlag == 1)
% Eingang:
%   xData   = Sampling-Punkte der Messung
%   Messung = Messung
%   dFlag   = Spezialfall "Eingang ist schon Differenzkurve"
%   Hub_min = Schwelle für rel. Unterschied zw. Maximum und benachbarten Minima
%                 der d/dx-Kurve zur Erkennung als Wendepunkt
% Ausgang:
%   x_w  = x-Achsen-Werte für die Wendepunkte
%   y_w  = Messwerte an den Wendepunkten
%   dy_w = d/dx(Messwert) an den Wendepunkten

if nargin < 4
  Hub_min = 0.19;
end

x_w  = NaN;
y_w  = NaN;
dy_w = NaN;

if (std (diff (xData)) > 0.005) || (median (diff (xData)) > 0.99)
  % Ungefilterte Messung muß gefiltert werden!
  [xData, Messung] = filtereRohdaten (xData, Messung);
end

if dFlag == 1 % Eingang ist schon Ableitung => direkte Maximasuche

  [x_w, y_w, dy_w] = bestimme_Wendepunkte (xData, Messung, Messung, Hub_min);

else % Maxima der gefilterten Ableitung von Messung gesucht

  % Filter für Tiefpass auf diff(Messung), erzeugt mit Matlab "fdatool":
  % Lowpass, FIR, Window, Order 64, Scale Passb., Hamming, wc 0.01 (normalized)
  b2    = [0.001993146074, 0.002071215513, 0.002262178644, 0.002567823732, ...
           0.002988588309, 0.003523532103, 0.004170327174, 0.004925265581, ...
           0.005783284650, 0.006738009699, 0.007781813837, 0.008905894214, ...
           0.010100363868, 0.011354358137, 0.012656154344, 0.013993303335, ...
           0.015352771243, 0.016721089727, 0.018084512794, 0.019429178237, ...
           0.020741271636, 0.022007190808, 0.023213708627, 0.024348132090, ...
           0.025398455582, 0.026353506350, 0.027203080302, 0.027938066343, ...
           0.028550557648, 0.029033948396, 0.029383014713, 0.029593978752, ...
           0.029664555079];

  b     = [b2, b2((end - 1):-1:1)].';
  nfilt = numel (b);

  dxData = mean (diff (xData));
  dM     = diff (Messung) / dxData;

  if numel (dM) < (3 * (nfilt - 1) + 1)
    return;
  end

  dF = filtfilt (b, 1, dM);
  [x_w, y_w, dy_w] = bestimme_Wendepunkte (xData, Messung, dF, Hub_min);

end

% Sichere Bestimmung mit >= 1 Woche Abstand zu den Rändern und
% bei y_w >= 1000 und dy_w >= 100!
for k = 1:numel (x_w)
  if ((x_w (k) - min (xData)) < 7) || ((max (xData) - x_w (k)) < 7) || ...
     (y_w (k) < 1000) || (dy_w (k) < 100)
    x_w  (k) = NaN;
    y_w  (k) = NaN;
    dy_w (k) = NaN;
  end
end

% behalte nur ~isnan, auf jeden Fall aber eines:
if numel (x_w) > 0
  maske = 1;
  if any (~isnan (x_w))
    maske = ~isnan (x_w);
  end
  x_w  = x_w  (maske);
  y_w  = y_w  (maske);
  dy_w = dy_w (maske);
else
  x_w  = NaN;
  y_w  = NaN;
  dy_w = NaN;
end

end

function [x_w, y_w, dy_w] = bestimme_Wendepunkte (xData, yData, dyData, Hub_min)
% Unterprogramm zu "berechneWendepunkt"
% Eingang:
%   xData   = Sampling-Punkte der Messung
%   yData   = Messung
%   dyData  = d/dx(Messwert)
%   Hub_min = Schwelle für rel. Unterschied zw. Maximum und benachbarten Minima,
%                 d.h. Schwelle für "(my - mn)/my"
%                 mit  my = Messung(Wendepunkt)
%                      mn = größeres der benachbarten relativen Minima
% Ausgang:
%   x_w  = x-Achsen-Werte für die Wendepunkte
%   y_w  = Messwerte an den Wendepunkten
%   dy_w = d/dx(Messwert) an den Wendepunkten

if nargin < 4
  Hub_min = 0.3;
end

% Spaltenvektoren erzwingen:
xData  = xData  (:);
yData  = yData  (:);
dyData = dyData (:);

index_w = find ((dyData (2:(end - 1)) >= dyData (1:(end - 2))) & ...
                (dyData (2:(end - 1)) >  dyData (3:end)));
AnzW    = numel (index_w);
if AnzW > 0 % ein Wendepunkt oder mehrere Wendepunkte
  index_n = zeros ((AnzW - 1), 1);
  for k = 1:(AnzW - 1)
    index_bereich  = index_w (k):index_w (k + 1);
    dyData_Bereich = dyData (index_bereich);
    idx            = ...
      find ((dyData_Bereich (2:(end - 1)) <  dyData_Bereich (1:(end - 2))) & ...
            (dyData_Bereich (2:(end - 1)) <= dyData_Bereich (3:end))) + ...
      index_w (k) - 1;
    if ~isempty (idx)
      index_n (k) = idx (1);
    else
      index_n (k) = index_w (k);
    end
  end
  [~, index_a] = min (dyData (1:index_w (1)));
  [~, index_e] = min (dyData ((index_w (end) + 1):end));
  index_e      = index_e + index_w (end);
  index_n      = [index_a; index_n; index_e];

  while true
    AnzW = numel (index_w);
    if AnzW == 0
      break;
    end
    hub1 = zeros (AnzW, 1);
    hub2 = zeros (AnzW, 1);
    for k = 1:AnzW
      hub1 (k) = (dyData (index_w (k)) - dyData (index_n (k)))     / ...
                  dyData (index_w (k));
      hub2 (k) = (dyData (index_w (k)) - dyData (index_n (k + 1))) / ...
                  dyData (index_w (k));
    end
    [minhub1, idxmin1] = min (hub1); %  linke Flanke
    [minhub2, idxmin2] = min (hub2); % rechte Flanke
    if     (minhub1 < minhub2) && (minhub1 < Hub_min)
      % links löschen
      index_w (idxmin1) = [];
      index_n (idxmin1) = [];
    elseif (minhub2 < minhub1) && (minhub2 < Hub_min)
      % rechts löschen
      index_w (idxmin2)     = [];
      index_n (idxmin2 + 1) = [];
    else
      break;
    end
  end

  dy_w  = dyData (index_w);

else

  index_w = [];
  dy_w  = zeros (size (index_w));

end

x_w = zeros (size (index_w));
y_w = zeros (size (index_w));
for k = 1:numel (index_w)
  x_w (k) = mean (xData (index_w (k) + [0, 1]));
  y_w (k) = mean (yData (index_w (k) + [0, 1]));
end

end

function LegPos = bestlegpos (hlines, haxes)

xlim = get(haxes, 'xLim');
ylim = get(haxes, 'yLim');
dxlim = diff(xlim);

logflag = false;
if strcmp(get(haxes,'yScale'), 'log')
  logflag = true;
  ylim = log10(ylim);
end

anzobenlinks   = 0;
anzuntenlinks  = 0;
anzobenrechts  = 0;
anzuntenrechts = 0;
anzmittelinks  = 0;
anzmitterechts = 0;
anzobenmitte   = 0;
anzuntenmitte  = 0;

margin = 0.025;
grenze = 0.35;

ylim_unten = grenze;
ylim_oben  = 1-grenze;
ylim_mitte = [grenze (1-grenze)];

xlim_links  = xlim(1) + [margin grenze]*dxlim;
xlim_rechts = xlim(2) - [grenze margin]*dxlim;
xlim_mitte  = xlim(1) + [grenze (1-grenze)]*dxlim;

for k=1:numel(hlines)
  xdata = get(hlines(k),'xData');
  ydata = get(hlines(k),'yData');
  if logflag
    ydata = log10(ydata);
  end
  ydata = (ydata - ylim(1))/diff(ylim);
  ydata(ydata<0) = NaN;
  ydata(ydata>1) = NaN;

  links  = (xdata > xlim_links(1) ) & (xdata <= xlim_links(2) );
  rechts = (xdata > xlim_rechts(1)) & (xdata <= xlim_rechts(2));
  mitte  = (xdata > xlim_mitte(1) ) & (xdata <  xlim_mitte(2) );

  anzobenlinks   = anzobenlinks   + sum(ydata(links)  > ylim_oben );
  anzuntenlinks  = anzuntenlinks  + sum(ydata(links)  < ylim_unten);
  anzobenrechts  = anzobenrechts  + sum(ydata(rechts) > ylim_oben );
  anzuntenrechts = anzuntenrechts + sum(ydata(rechts) < ylim_unten);
  anzmittelinks  = anzmittelinks  + sum(ydata(links)  > ylim_mitte(1) & ...
                   ydata(links)  < ylim_mitte(2));
  anzmitterechts = anzmitterechts + sum(ydata(rechts) > ylim_mitte(1) & ...
                   ydata(rechts) < ylim_mitte(2));
  anzobenmitte  = anzobenmitte  + sum(ydata(mitte)   > ylim_oben  );
  anzuntenmitte = anzuntenmitte + sum(ydata(mitte)   < ylim_unten );
end

anzahl = [anzobenlinks , anzobenrechts , ...
          anzmittelinks, anzmitterechts, ...
          anzuntenlinks, anzuntenrechts, ...
          anzobenmitte,  anzuntenmitte];

[~, idx] = min(anzahl);

LegPosNames = {'NorthWest', 'NorthEast', 'West', 'East', ...
               'SouthWest', 'SouthEast','North','South'};
LegPos = LegPosNames{idx};

end

function stretch_legend (hl, xboxstretch, location)
% Vermeidet das Hinausragen des Textes aus dem Plot-Legende-Rahmen,
% wenn das Bild gespeichert wird.

set (hl, 'FontSize',9);

if ~isfield (get (hl), 'outerposition')
  return; % Matlab hat das nicht / braucht keine Aktion hier.
end

op = get (hl, 'outerposition');
if     strfind (lower (location), 'east')
  op (1) = op (1) - (xboxstretch - 1) * op (3);
elseif strfind (lower (location), 'west')
  % nothing else
else
  op (1) = op (1) - (xboxstretch - 1) * op (3) / 2;
end
op (3) = xboxstretch * op (3);
set (hl, 'outerposition',op);

end

function new_axis = scale_axis (nticks, ha, axname)

if nargin < 2
  ha = gca;
end
if nargin < 3
  axname = 'y';
end

if     strcmpi (axname (1), 'x')
  lim = get (ha, 'xLim');
elseif strcmpi (axname (1), 'y')
  lim = get (ha, 'yLim');
else
  error ('%s: Axis specification = %s', mfilename, axname);
end
scale_region = 10.^(-1:floor (log10 (lim (2))));
scale_values = [1, 2, 2.5, 4, 5];
scale_values = (scale_region.' * scale_values).';
scale_values = (scale_values (:)).';
new_lim      = [(floor(lim(1) ./ scale_values) .* scale_values); ...
                (ceil(lim(2)  ./ scale_values) .* scale_values)];
delticks     = diff (new_lim) ./ scale_values;
[~, idx]     = min (abs (delticks - nticks));
delticks     = scale_values (idx);
new_lim      = [floor(lim(1) / delticks), ceil(lim(2)  / delticks)] * delticks;
new_axis     = new_lim (1):delticks:new_lim (2);
new_axis (new_axis < lim (1)) = [];
new_axis (new_axis > lim (2)) = [];
if strcmpi (axname (1), 'y')
  set (ha, 'yTick',new_axis);
else
  set (ha, 'xTick',new_axis);
end

end

function skaliere_Wochen (ha, x_min, x_max, Schwelle_AnzKW)

if nargin < 4
  Schwelle_AnzKW = [10 80];
end

if     (x_max - x_min) > Schwelle_AnzKW (2)
  Skalenwerte = ceil (x_min):4:floor (x_max);
  set (ha, 'xTick',Skalenwerte);
  naechstesJahr                    = find (Skalenwerte >  53);
  uebernaechstesJahr               = find (Skalenwerte > 105);
  Skalenwerte (naechstesJahr)      = Skalenwerte (naechstesJahr)      - 53;
  Skalenwerte (uebernaechstesJahr) = Skalenwerte (uebernaechstesJahr) - 52;
  set (ha, 'xTickLabel',str2cell (num2str (Skalenwerte.')));
elseif (x_max - x_min) > Schwelle_AnzKW (1)
  Skalenwerte = ceil (x_min):floor (x_max);
  set (ha, 'xTick',Skalenwerte);
  naechstesJahr                    = find (Skalenwerte >  53);
  uebernaechstesJahr               = find (Skalenwerte > 105);
  Skalenwerte (naechstesJahr)      = Skalenwerte (naechstesJahr)      - 53;
  Skalenwerte (uebernaechstesJahr) = Skalenwerte (uebernaechstesJahr) - 52;
  set (ha, 'xTickLabel',str2cell (num2str (Skalenwerte.')));
else
  Skalenwerte = (ceil (x_min * 7) / 7):(1 / 7):(floor (x_max * 7) / 7);
  set (ha, 'xTick',Skalenwerte);
  naechstesJahr                    = find (Skalenwerte >  53);
  uebernaechstesJahr               = find (Skalenwerte > 105);
  Skalenwerte (naechstesJahr)      = Skalenwerte (naechstesJahr)      - 53;
  Skalenwerte (uebernaechstesJahr) = Skalenwerte (uebernaechstesJahr) - 52;
  set (ha, 'xTickLabel',str2cell (num2str (Skalenwerte.')));
  if (x_max - x_min) > 0.9
    xtl  = get_ticklabel (ha, 'x');
    nums = str2num (char (xtl));
    xtl (nums ~= round (nums)) = {' '};
    set_ticklabel (ha, 'x', xtl);
  end
end

if floor (x_max) <= 53
  xlabel ('Zeit [Kalenderwoche im Jahr 2020]');
else
  if floor (x_max) <= 105
    if ceil (x_min) <= 53
      xlabel ('Zeit [Kalenderwoche in 2020/2021]');
    else
      xlabel ('Zeit [Kalenderwoche im Jahr 2021]');
    end
  else
    if     ceil (x_min) <=  53
      xlabel ('Zeit [Kalenderwoche in 2020-2022]');
    elseif ceil (x_min) <= 105
      xlabel ('Zeit [Kalenderwoche in 2021/2022]');
    else
      xlabel ('Zeit [Kalenderwoche im Jahr 2022]');
    end
  end
end

end

function thin_labels (ha, axname, numLabels)

if     strcmpi (axname (1), 'x')
  x = 'X';
elseif strcmpi (axname (1), 'y')
  x = 'Y';
elseif strcmpi (axname (1), 'z')
  x = 'Z';
else
  error ('%s: Axis specification = %s', mfilename, axname);
end

Labels    = get_ticklabel (ha, axname);
AnzLabels = numel (Labels);
for k = 1:AnzLabels
  if Labels {k} == ' '
    return; % ist schon ausgedünnt!
  end
end

factor    = max (1, round (AnzLabels / numLabels));
lab2num   = floor (str2num (char (Labels {1})));
start     = (factor - mod (lab2num, factor));

maske     = false (1, AnzLabels);
maske (start:factor:end) = true;
for k = 1:AnzLabels
  if ~maske (k)
    Labels {k} = ' ';
  end
end
set_ticklabel (ha, x, Labels);

end

function octave_label (ha, axname)
% Ändert die Label vom %g-Format in das %d-Format für große Zahlen:
% Ansonsten kollidiert der Exponent mit der Bildüberschrift
% (eine Matlab Core-Inkompetenz).

if ~isempty (strfind (license, 'GNU'))
  return;
end

if     strcmpi (axname (1), 'x')
  x = 'X';
elseif strcmpi (axname (1), 'y')
  x = 'Y';
elseif strcmpi (axname (1), 'z')
  x = 'Z';
else
  error ('%s: Axis specification = %s', mfilename, axname);
end

Axis       = get  (ha, [x, 'Axis']);
TickValues = get  (Axis, 'TickValues');

if ~all (isint (TickValues))
  return;
end

TickLabel  = cell (size (TickValues));
for k = 1:numel (TickValues)
  TickLabel {k} = sprintf ('%d', TickValues (k));
end

set_ticklabel (ha, x, TickLabel);
% Danach nicht mehr umskalieren, sonst kommt Chaos heraus!

end

function cut_last_label (ha, axname)

if     strcmpi (axname (1), 'x')
  x = 'X';
elseif strcmpi (axname (1), 'y')
  x = 'Y';
elseif strcmpi (axname (1), 'z')
  x = 'Z';
else
  error ('%s: Axis specification = %s', mfilename, axname);
end

Ticks  = get (ha, [x, 'Tick']);
Limits = get (ha, [x, 'Lim']);
Labels = get (ha, [x, 'TickLabel']);
if Ticks (end) > (Limits (2) - (Ticks (end) - Ticks (end - 1)) / 4)
  Labels {end} = [];
  set (ha, [x, 'TickLabel'],Labels);
end

end

function ticklabel = get_ticklabel (ha, axname)
% Diese Funktion ist für Matlab nötig, da dort mit "get (ha, 'xTickLabel')"
% nur die gerade sichtbaren Label zurückgegeben werden, anders als in Octave.
% Diese Funktion gibt also immer alle definierten Label zurück:
%    "get_ticklabel (ha, 'x')" oder "get_ticklabel (ha)" die Label der x-Achse,
%    "get_ticklabel (ha, 'y')" die Label der y-Achse und
%    "get_ticklabel (ha, 'z')" die Label der z-Achse.
%    "ha" ist der axis handle, z.B. "gca".

if nargin < 1
  ha = gca;
end
if nargin < 2
  axname = 'x';
end

if     strcmpi (axname (1), 'x')
  x = 'X';
elseif strcmpi (axname (1), 'y')
  x = 'Y';
elseif strcmpi (axname (1), 'z')
  x = 'Z';
else
  error ('%s: Axis specification = %s', mfilename, axname);
end

lim       = get (ha, [x, 'Lim']);
ticks     = get (ha, [x, 'Tick']);
maske     = (ticks >= lim (1)) & (ticks <= lim (2));
set (ha, [x, 'Tick'],ticks (maske));
ticklabel = get (ha, [x, 'TickLabel']);

end

function set_ticklabel (ha, axname, ticklabel)
% Diese Funktion ist für Matlab nötig, da dort mit "set (ha, 'xTickLabel')"
% nur die gerade sichtbaren Label gesetzt werden, anders als in Octave.
% Diese Funktion setzt also immer alle definierten Label:
%    "set_ticklabel (ha, 'x', ticklabel)" die Label der x-Achse,
%    "set_ticklabel (ha, 'y', ticklabel)" die Label der y-Achse und
%    "set_ticklabel (ha, 'z', ticklabel)" die Label der z-Achse.
%    "ha" ist der axis handle, z.B. "gca".

if     strcmpi (axname (1), 'x')
  x = 'X';
elseif strcmpi (axname (1), 'y')
  x = 'Y';
elseif strcmpi (axname (1), 'z')
  x = 'Z';
else
  error ('%s: Axis specification = %s', mfilename, axname);
end

lim   = get (ha, [x, 'Lim']);
ticks = get (ha, [x, 'Tick']);
maske = (ticks >= lim (1)) & (ticks <= lim (2));
set (ha, [x, 'Tick'],ticks (maske));
set (ha, [x, 'TickLabel'],ticklabel);

end

function adjust_title_position (ha, normpos, liccond)

if ~isempty (strfind (license, liccond))
  t = get (ha, 'Title');
  set (t, 'Units','normalized');
  set (t, 'Position',normpos);
end

end

function limit_text_to_axis (hText, hAxes, MarginFactor)
% schneidet den Text ab, der über den rechten Bildrand hinausgeht
% input: hText = Text-Handle
% input: hAxes = Axis-Handle

if nargin < 3
  MarginFactor = 0.01;
end

xl      = get (hAxes,'XLim');
xmargin = xl (2) + MarginFactor * diff (xl);

for k = 1:numel (hText)
  while true
    ext = get (hText (k),'Extent');
    txt = get (hText (k),'String');
    if (ext (1) + ext (3)) <= xmargin
      break;
    end
    if numel (txt) < 1
      break;
    end
    set (hText (k),'String', txt (1:(end - 1)))
  end
end

end

function c = str2cell (s)

nr = size (s, 1);
c  = cell (nr, 1);
for k = 1:nr
  c {k} = strtrim (s (k, :));
end

end