gpxconv

#!/usr/bin/perl
################################################################################
# GPX track converter: post-process routes and tracks produced by GPS loggers etc.
#
# Reads from file(s) given as argument or STDIN and writes to STDOUT or a file.
# Concatenates segments of multiple given tracks or routes (sequential composition),
# using (header) metadata of the first track, preserving segments, and collecting waypoints.
# Optionally augments main input by merging data from alternative input (from parallel tracks),
# ignoring metadata and segment information of the alternative input.
#
# Reports statistics including total length and time, moving time, average moving speed and max speed,
# min and max altitude, total ascent and descent, as well as max ascent and descent rate.
# Reports missing data and by default ignores points containing implausible data w.r.t.
# altitude, ascent/descent rate, speed, acceleration/deceleration, and direction change.
# Between segments and on start/end of merged sections, speed and ascent/descent rate are ignored.
# Corrects implausible elevation gain/loss after gaps, typically on exit of tunnel or building
# Optionally filters out points before or after given time limits.
# Optionally filters out points with an element value below or above given limits.
# Optionally reduces output size by removing points that are rather close-by or on rather straight line.
# Optionally prunes comments, extensions, or extension elements with value below or above given limits.
# By default carries over missing altitude and time data between segments.
# By default fills in missing altitude and time data by interpolation within segments.
# By default inserts interpolated points in long time gaps (default 1800 seconds sufficient for exiftool).
# Optionally smoothens tracks.
# Optionally corrects elevation w.r.t. geoid height (= ellipsoidal - orthometric height).
# Optionally produces additional statistics on a per-segment or per-day basis.
# Optionally analyzes climbing phases with extra ascent/descent and speed statistics.
# Optionally calculates approximate total energy spent by biking along the resulting track.
# Prints information and warnings (unless disabled), as well as any errors, to STDERR.
#

use constant TOOL_NAME    => "GPXConv";
use constant TOOL_VERSION => "2.9.3";
my $TOOL_ID = TOOL_NAME." v".TOOL_VERSION;
use constant MIN_INPUT_TIMEDIFF => 2;

my $usage =
"Usage:  gpxconv <option(s) and infile(s)> [> <outfile>] [2> <logfile>]

Command-line options:
  -h | -help | --help          - print these hints and exit
  -version                     - print tool name and version and exit
  -swim, -walk, -cycle, -drive, -fly - type of recorded activity, default: drive
  -no_sanitize                 - do not sanitize trackpoints with implausible data
  -no_insert                   - do not insert trackpoints on long time gaps
  -smooth                      - smoothen tracks
  -phases                      - analyze and provide statistics for ascent/descent phases
  -segs [<n>..[<m>][(+|-)<d>]] - add statistics per segment n..m, may adapt indexes +/-d
  -days [<n>..[<m>][(+|-)<d>]] - add statistics per day n..m, may adapt indexes +/-d
  -tz <hours>                  - offset from UTC (timezone+DST) for -days, default taken from system
  -split <basename>            - produce GPX files per day or segment, with name <basename>_<i>.gpx
  -lean_stat_wpts              - only one stat. wpt per segment/day with data except ascent/descent
  -info_wpts                   - provide info points also in the form of waypoints
  -warn_wpts                   - provide warning points also in the form of waypoints
  -cmt_wpts <pattern>          - add info waypoints for trackpoints with matching 'cmt' element
  -src_wpts <pattern>          - add warning waypoints for trackpoints with matching 'src' element
  -merge <file>                - augment tracks in infile(s) with data from given alternative file
  -weight <value>              - calculate biking energy for given weight in kg
  -begin <time>                - ignore points before given time
  -end <time>                  - ignore points after given time
  -resolution <seconds>        - minimum time difference between trackpoints on input, default: ".MIN_INPUT_TIMEDIFF."
  -min <name> <limit>          - filter out points with element value below limit
  -max <name> <limit>          - filter out points with element value above limit
  -min_ext <name> <limit>      - prune extension elements with value below limit
  -max_ext <name> <limit>      - prune extension elements with value above limit
  -min_dist <meters>           - output only points that have at least the given distance
  -min_angle <degree>          - output only points after at least the given deviation angle
  -prune_wpts                  - remove waypoints
  -prune_cmts                  - remove comments
  -prune_exts                  - remove extensions
  -ni                          - do not output information
  -nw                          - do not output warnings
  -o <outfile>                 - output the resulting track to <outfile>
  -debug                       - enable internal consistency checks and debug output
Other options can be set by changing values of the configuration variables.

(c) 2012-2025 David von Oheimb - License: MIT - ".TOOL_NAME." version ".TOOL_VERSION;
################################################################################

use strict;
use warnings;
# use feature 'signatures'; # https://perldoc.perl.org/perlsub#Signatures

### configuration variables

my $activity = "-drive"; # default

# basic output control
use constant PRUNE_WPTS           => 0; # by default do not ignore waypoints
use constant PRUNE_AUTOMATIC_WPTS => 1; # ignore waypoints containing <name><!\[CDATA\[((day )?\d+)?\]\]></name> as produced by OruxMaps
use constant PRUNE_CMTS           => 0; # by default do not remove comments
use constant PRUNE_DESCS          => 0; # do not remove descriptions
use constant PRUNE_LINKS          => 0; # do not remove links
use constant PRUNE_EXTENSIONS     => 0; # by default do not remove extensions
use constant PRUNE_RECORDED_SPEED => 1; # remove existing speed extensions with recorded speed (<gpxtpx:speed>)
use constant INCLUDE_SPEED        => 1; # include our simple speed extensions with recorded or computed speed (<speed>)
use constant DISABLE_TRKPT_CMT    => 0; # do not add 'cmt' entries to trackpoints, e.g., comments on ignored points
use constant DISABLE_TRKPT_SRC    => 0; # do not add 'src' entries to trackpoints, e.g., on time gaps and interpolated data
use constant DISABLE_TRKPT_FIX    => 0; # do not add 'fix' entries to trackpoints, on points inserted with interpolated data

# control filling in missing altitude and/or time
use constant CARRY_OVER_ELE    => 1; # copy missing altitude from end of previous segment or begin of following segment # better not: last available value (prev_avail_ele)
use constant CARRY_OVER_TIME   => 1; # copy missing time     from end of previous segment or begin of following segment
use constant FILL_MISSING_ELE  => 1; # interpolate missing altitude
use constant FILL_MISSING_TIME => 1; # interpolate missing time

# control use of recorded vs. computed speed
use constant MAX_TIMEDIFF_RECORDED_SPEED => 2; # maximum number of seconds for which recorded speed is used

# analysis control
use constant MIN_SPEED_MOVING =>  1; # threshold for detecting movement (in km/h)
use constant ANALYZE_PHASES   =>  0; # by default do not analyze ascent/descent phases
use constant ELE_THRESHOLD    => 25; # threshold for detecting/accepting ascent or descent (in m);
                                     # should be larger than short-term altitude measurement error
my $PHASES_REPORT_THRESHOLD = ELE_THRESHOLD; # default

# merge control
use constant MERGE_MIN_GAP    =>  20; # minimum gap duration in seconds for starting to switch to alternative input
use constant MERGE_MIN_MORE   =>   7; # minimum number of more points on alternative input to allow switching to it
use constant MERGE_WARN_DIST  => 100; # warn if first point of merged portion is more than given number of meters away from regular track

# smoothing control
use constant ENABLE_SMOOTHING  =>  0; # by default disable smoothing
use constant SMOOTHING_MAX_GAP => 60; # maximum number of seconds between trackpoint to be smoothened and its neighbors

# control insertion of points within long time gaps
my $insert_points                =     1; # whether to enable inserting points with interpolated data in long time gaps
use constant INSERT_SEGMENT      =>    0; # enable inserting segment for long time gaps just before segment start
use constant INSERTION_MAX_GAP   => 1800; # length of gap in seconds between consecutive trackpoints before
# insertion should be used (1800 for exiftool default GeoMaxIntSecs; for TrailGuru 3600 would be sufficient)

# info warning output control
my $output_info                  =     1; # whether outputting information is enabled
my $output_warnings              =     1; # whether outputting warnings is enabled
my $info_wpts                    =     0; # whether to output information points in addition as waypoints
my $warn_wpts                    =     0; # whether to output warning     points in addition as waypoints
my $cmt_wpts = ""; # pattern to match for adding info waypoints for trackpoints with matching 'cmt' element
my $src_wpts = ""; # pattern to match for adding wanrning waypoints for trackpoints with matching 'src' element
use constant WARNING_TPT_DIST    => 2000; # threshold for trackpoint distance warning (in meters)
use constant WARNING_WPT_DIST    =>  100; # threshold for waypoint distance warning (in meters)
use constant MAX_SPEED_DEVIATION => 0.25; # threshold for speed measuring deviation warning (fraction of MAX_PLAUSIBLE_SPEED)
# using absolute value for MAX_SPEED_DEVIATION because we compute speed from absolute values (coordinates and time)

# sanitization control
use constant   SANITIZE_TRKPTS   =>    1; # whether trackpoint sanitization is enabled by default
my $sanitize = SANITIZE_TRKPTS;
my $MIN_TIMEDIFF          ; # in seconds
my $MAX_PLAUSIBLE_SPEED   ; # maximal speed in km/h
my $MAX_PLAUSIBLE_ACCEL   ; # maximal acceleration m/s/s
my $MAX_PLAUSIBLE_DECEL   ; # maximal deceleration m/s/s
my $MAX_PLAUSIBLE_ELE_GAIN; # maximal ascent rate in m/h
my $MAX_PLAUSIBLE_ELE_LOSS; # maximal descent rate in m/h
use constant MIN_TUNNEL_GAP           =>   20; # threshold for assuming exit of tunnel or building before elevation gain/loss correction (in seconds)
use constant MAX_PLAUSIBLE_ANGLE_DIFF =>  170; # maximal plausible turning angle above speed threshold
use constant MAX_PLAUSIBLE_ANGLE_SPD_THRESHOLD => 15; # speed threshold for maximal plausible turning angle in km/h
use constant MIN_PLAUSIBLE_ELE        => -450; # minimal plausible altitude; actual elevation values on Earth may be as low as -450 m at Dead Sea

# time control
use constant TIME_CORRECTION => 0; # number of seconds to add to time stamps in trkpt and wpt
my $MAX_UNDEF_TIME = str_to_epoch("9999-12-31T23:59:59Z"); # as used in RFC 5820 section 4.1.2.5

# elevation correction control
use constant GEOID_ELE_CORRECTION   => 0; # whether to correct elevation wrt. geoid height retrieved online
use constant DEFAULT_ELE_CORRECTION => 0; # 47; # if used should be -(average geoid height)

### other constants

use constant H2S            =>    3600 ; # factor for converting hours to seconds
use constant M2KM           =>    0.001; # factor for converting m to km
use constant MPS2KMPH       =>H2S * M2KM;# factor for converting m/s to km/h
use constant METERS_PER_DEGREE_LAT => 111133; # at 45 degrees, cf. https://en.wikipedia.org/wiki/Latitude#Length_of_a_degree_of_latitude
use constant LAT_PRECISION  =>   "%.5f"; #  latitude resolution = 0.00001 degrees (<= 1.11 meters)
use constant LON_PRECISION  =>   "%.5f"; # longitude resolution = 0.00001 degrees (<= 1.11 meters)
use constant ELE_PRECISION  =>  "%5.0f"; #  altitude resolution = 1 meters
use constant ELE_PRECISION0 =>   "%.0f";
use constant ELE_PRECISION3 =>   "%.3f";
use constant DIF_PRECISION  => "%+5.0f"; #  altitude difference          resolution = 1 meters
use constant SLOPE_PRECISION=>  "%2.0f"; #     slope resolution = 1 %
use constant RAT_PRECISION  => "%+5.0f"; #  altitude ascent/descent rate resolution = 1 meters
use constant RAT_PRECISION2 => "%+6.0f"; #  altitude ascent/descent rate resolution = 1 meters
use constant DIS_PRECISION  =>   "%.3f"; #  distance resolution = 0.001 km
use constant DIS_PRECISION2 =>  "%6.3f"; #  distance resolution = 0.001 km
use constant SEC_PRECISION  =>   "%.0f"; #   seconds resolution = 1 second
use constant SPD_PRECISION  =>  "%4.0f"; #     speed resolution = 1 km/h
use constant AVG_PRECISION  =>   "%.1f"; # avg speed resolution = 0.1 km/h
use constant SPD_PRECISION0 =>   "%.0f";
use constant ACC_PRECISION  =>   "%.1f"; # acceleration resolution = 0.1 m/s/s

use constant EMOJI_PREFIX   =>       ""; # or "emoji-" for GAIA GPS

# trkpt flags
use constant DUPLICATED_POINT   => 1 << 0; # no more actually uused
use constant INSERTED_POINT     => 1 << 1; # its use could be removed
use constant INTERPOLATED_ELE   => 1 << 2;
use constant INTERPOLATED_TIM   => 1 << 3;
use constant SUBSTITUTED_ELE    => 1 << 4;
use constant SUBSTITUTED_TIM    => 1 << 5;
use constant COMPUTED_SPEED     => 1 << 6;
use constant CORRECTED_ELE      => 1 << 7;
use constant START_MERGE        => 1 << 8;
use constant PAUSE_MERGE        => 1 << 9;
use constant STATISTICS_WPT     => 1 <<10; # used to prevent removal on pruning trkpts

# global option variables with related num_* state variables

my $debug = 0;
my $merge;
my $num_ignored_diff = 0;
my $num_ignored_hdrs = 0;
my $begin_sec;
my $num_ignored_before = 0;
my $end_sec;
my $min_input_timediff = MIN_INPUT_TIMEDIFF;
my $num_ignored_after  = 0;
my $max_elem;
my $max_elem_limit;
my $num_elem_above  = 0;
my $min_elem;
my $min_elem_limit;
my $num_elem_below  = 0;
my $max_ext;
my $max_ext_limit;
my $num_ext_above  = 0;
my $min_ext;
my $min_ext_limit;
my $num_ext_below  = 0;
my $min_dist;
my $num_removed_min_dist = 0;
my $min_angle;
my $num_removed_min_angle = 0;
my $num_points_output = 0;
my $ele_corr         = DEFAULT_ELE_CORRECTION;
my $smoothing        = ENABLE_SMOOTHING;
my $num_smoothened   = 0;
my $phases           = ANALYZE_PHASES;
my $weight;
my $prune_wpts       = PRUNE_WPTS;
my $num_pruned_wpts  = 0;
my $num_auto_wpts    = 0;
my $prune_cmts       = PRUNE_CMTS;
my $num_pruned_cmts  = 0;
my $num_pruned_descs = 0;
my $num_pruned_links = 0;
my $prune_exts       = PRUNE_EXTENSIONS;
my $num_pruned_exts  = 0;

# trackpoint data from regular infile(s), followed by data from alternative input
# array with doubly linked list structure for efficient element insertion and deletion
my @PREV; # index of point before (initially current index - 1)
my @NEXT; # index of point after  (initially current index + 1)
sub PREV { return defined $_[0] ? $PREV[$_[0]] : undef; }
sub NEXT { return defined $_[0] ? $NEXT[$_[0]] : undef; }
my @TXT; # text (containing trkpt output) to prepend when printing
my @IGN; # number of points ignored immediately before this one
my @SEG;
my @FLG;
my @LAT;
my @LON;
my @ELE; # or "" if not available
my @TIM; # or "" if not available
my @SEC; # computed from @TIM, possibly adapted, or undefined if not available
my @DIS2;# square sum of the LAT and COS(LAT)*LON differences from point before, 0 for first point
my @DIS; # derived from LAT and LON, and ELE if available, also from point before, 0 for first point
my @ANG; # horizontal direction change (i.e., turning angle) between last, current, and next point, undef if not available
my @DIR; # the direction arriving at each point from the previous one, undef if not (yet) computed or invalidated
my @SPD; # as recorded or computed (average from previous point if FLG contains COMPUTED_SPEED), or "" if not available
my @CMT; # as given in input but without <cmt> and </cmt> tags, or "" if not available
my @SRC; # as given in input but without <src> and </src> tags, or "" if not available
my @FIX; # as given in input but without <fix> and </fix> tags, or "" if not available
my @EXT; # or "" if not available

# waypoint data

my @WLAT;
my @WLON;
my @WELE; # or "" if not available
my @WTIM; # or "" if not available
my @WSEC; # computed from @WTIM, possibly adapted, or undefined if not available
my @WTXT; # inner waypoint text after </time>, or "" if not available
my @WSTR; # for diagnostics only: full original but indented waypoint text surrounded by '...', or undefined for statistics waypoint
my @WIDX; # index of trkpt related to this wpt belongs, or -1

# ascent/descent phase data

my @PHASE_DURATION;
my @PHASE_DIFF;
my @PHASE_DIST;
my @PHASE_SPD;
my @PHASE_RATE;
my @PHASE_MAX_SPD_INDEX;
my @PHASE_MAX_RATE;
my @PHASE_MAX_RATE_INDEX;
my @PHASE_END_INDEX;
my $ascent__phases_suppressed = 0;
my $descent_phases_suppressed = 0;

# global state variables

my $HEAD;
my @NAME;
my $found_corr = 0;

my $num_merges = 0;
my $num_merged_points = 0;
my $num_skip_points = 0;
my $num_merged_points_last = 0;
my $num_skip_points_last = 0;

my $min_lat, my $max_lat;
my $min_lon, my $max_lon;

my ($min_sec, $max_sec);
my ($min_tim, $max_tim) = ("", ""); # $max_tim is not really used

my ($tz, $sec_offset);
my ($segs, $days);
my ($part_start, $part_end, $part_offset);
my $split;
my $lean_stat_wpts = 0;
# the following arrays hold statistic data per part, with index ($part_start - 1) holding the overall values
my (@MAX_SPD      , @MIN_ELE      , @MAX_ELE      , @MAX__ASC,        @MAX_DESC      , @MAX_GAIN      , @MAX_LOSS      );
my (@MAX_SPD_INDEX, @MIN_ELE_INDEX, @MAX_ELE_INDEX, @MAX__ASC_INDEX,  @MAX_DESC_INDEX, @MAX_GAIN_INDEX, @MAX_LOSS_INDEX);
my (@LAST_GAIN_INDEX, @LAST_LOSS_INDEX, @PART_START_INDEX, @PART_END_INDEX);
my (@SUM__ASCENT, @TIME__ASCENT, @MISSING__ASCENT_TIME,
    @SUM_DESCENT, @TIME_DESCENT, @MISSING_DESCENT_TIME);
my (@SUM_DIS, @SUM_DIS_MOV, @SUM_TIMEDIFF_MOV, @DURATION);

 # for calculating relative deviation between recorded and computed speed
my $sum_speed = 0;
my $sum_speed_deviation = 0;

my $sum_energy = 0;

my $out = *STDOUT;  # https://stackoverflow.com/questions/10478884/are-there-rules-which-tell-me-what-form-of-stdout-stderr-sdtin-i-have-to-choose

### various subprocedures

#http://www.perlmonks.org/?node_id=406883
#sub max { return $_[0] > $_[1] ? $_[0] : $_[1]; }
use List::Util qw[min max];

# http://stackoverflow.com/questions/178539/how-do-you-round-a-floating-point-number-in-perl
use Math::Round qw(nearest);
use Math::Trig  qw(deg2rad pi2);   # use Math::Trig 'great_circle_distance';
use constant    RAD2DEG => 360 / pi2;
use File::Temp  qw(tempfile);

# use DateTime::Format::ISO8601;
use Time::ParseDate ();

# use Time::PrintDate;
use Time::gmtime qw(gmtime);

use Scalar::Util qw(looks_like_number);

sub print_line { print STDERR "$_[0]\n"; }

sub warning {
    print_line("WARNING: $_[0]") if $output_warnings;
}

sub info {
    print_line("INFO   : $_[0]") if $output_info;
}

sub abort {
    print_line("$_[0]\naborting");
    exit 1;
}

sub debug {
    print_line("### DEBUG: $_[0]");
}

sub debug_log {
    print_line("### LOG: $_[0]") if $debug;
}

# https://stackoverflow.com/questions/229009/how-can-i-get-a-call-stack-listing-in-perl
# alternatively to making explicit individual calls to stacktrace, may use: perl -d:Confess
sub stacktrace {
    use Devel::StackTrace;
    print_line("");
    debug(Devel::StackTrace->new->as_string);
    # { # further alternative:
    #     use Carp qw<longmess>;
    #   # local $Carp::CarpLevel = -1;
    #     debug longmess();
    # }
    # { # further alternative:
    #     for (print STDERR "Stack Trace:\n", my $i = 0; my @call_details = caller($i); $i++) {
    #         print STDERR $call_details[1].":".$call_details[2]." in function ".$call_details[3]."\n";
    #     }
    # }
}

# str_to_epoch("1970-01-01T00:00:00Z") = 0; returns undef on error
# due to parsedate(), allows for weird values, e.g., 2000-01-01T00:00:63Z is taken as 2000-01-01T00:01:03Z
sub str_to_epoch {
    stacktrace() unless $_[0] eq '' || $_[0];
    my $s = $_[0];
    return undef unless $s =~ m/^\d{2,4}-\d{1,2}-\d{1,2}(T(\d{1,2}(:\d{1,2}(:\d{1,2}(\.\d{1,})?)?)?)?)?\s*([A-Z]+|[+-]\d+)$/; # pre-parse because parsedate() is not very strict

# within a day, may use: seconds + 60 * (minute + 60 * hour)
# return DateTime::Format::ISO8601->parse_datetime($s)->epoch();
# # not used due to Perl library bug on Mac OS X: "dyld: lazy symbol binding failed: Symbol not found: _Perl_newSVpvn_flags"
# # http://www.en8848.com.cn/Reilly%20Books/perl3/cookbook/ch03_08.htm
# use Time::Local;
# $date is "1998-06-03" (YYYY-MM-DD form).
# ($yyyy, $mm, $dd) = ($date =~ /(\d+)-(\d+)-(\d+)/;
# # calculate epoch seconds at midnight on that day in this timezone
# $epoch_seconds = timegm(0, 0, 0, $dd, $mm, $yyyy);
    $s =~ s/-/\//g;
    $s =~ s/(T\d+)Z/$1:00Z/; # workaround for the case that only hour is given
    $s =~ s/T/ /;
    $s =~ s/Z/+0000/;
    my $sec = Time::ParseDate::parsedate($s);
    $sec += $1 if defined $sec && $s =~ m/:\d\d(\.\d+)/; # fractional seconds
    return $sec;
}

sub epoch_to_str {
    stacktrace() unless defined $_[0];

# use DateTime; # not used due to Perl library bug on Mac OS X:
                # "dyld: lazy symbol binding failed: Symbol not found: _Perl_newSVpvn_flags"
#  my $dt = DateTime->from_epoch(epoch => $_[0]);
#  return $dt->ymd."T".$dt->hms."Z";
#use Date::Manip qw(ParseDate UnixDate);
#$date = ParseDate("18 Jan 1973, 3:45:50");
#  return UnixDate($_[0], "%Y-%m-%dT%H:%M:%SZ");
    my $sec = shift;
    my $millisec = nearest(.001, $sec - int($sec)); # fractional seconds
    my $millisecs = substr(substr($millisec, 1)."00", 0, 4);
    my $tm = gmtime($sec);
    return sprintf(
        "%04d-%02d-%02dT%02d:%02d:%02d%sZ",
        $tm->year + 1900,
        $tm->mon + 1,
        $tm->mday, $tm->hour, $tm->min, $tm->sec,
        $millisecs eq "00" ? "" : $millisecs);
}

sub round_tim { # round to nearest second
    stacktrace() unless $_[0];
    my $tim = shift;
    # $tim =~ s/(:\d\d)(\.\d+)/$1/;  # strip fractional seconds
    return $tim =~ m/\.(\d+)Z$/ && $1 != 0 ? epoch_to_str(nearest(1, str_to_epoch($tim))) : $tim;
}

sub timediff_str {
    stacktrace() if !defined $_[1] || $_[1] eq "";
    my $t = $_[1] + 0.5; # for rounding to nearest second
    my $s = $t % 60;
    $t = ($t - $s) / 60;
    my $m = $t % 60;
    $t = ($t - $m) / 60;
    return sprintf($_[0].":%02d:%02d h", $t, $m, $s);
}
sub timediff_string { return timediff_str("%d", $_[0]); }

sub ele_string {
    my ($len, $ele) = @_;
    stacktrace() if !defined $ele || $ele eq "";
    $ele = (sprintf ELE_PRECISION0, $ele) =~ s/^\s*-(0(\.0*)?)$/$1/r if $len >= 0;
    return (" " x max(0, abs($len) - length($ele))) . $ele;
}
sub dis_str {
    stacktrace() if !defined $_[1] || $_[1] eq "";
    return sprintf $_[0]." km", M2KM * $_[1];
}
sub dis_string { return dis_str(DIS_PRECISION, $_[0]); }

sub point_str {
    my $i = $_[0];
    return "point $i: ". ($TIM[$i] ne ""
        ? "tim=$TIM[$i]"
        : "($LAT[$i], $LON[$i], ".($ELE[$i] eq "" ? "none" : $ELE[$i]).")");
}

sub spd_or_none { return $_[0] eq "" ? "none" : (sprintf "%.3f", MPS2KMPH * $_[0]) =~ s/\.?0+$//r; }

my $spd_prec_len = 1; # used only for spd_string
sub spd_full {
    stacktrace() if !defined $_[0] || $_[0] eq "";
    return sprintf SPD_PRECISION, MPS2KMPH * $_[0];
}
sub spd_str {
    stacktrace() if !defined $_[0] || $_[0] eq "";
    return sprintf SPD_PRECISION0, MPS2KMPH * $_[0];
}
sub spd_string {
    my $s = spd_str($_[0]);
    return " " x ($spd_prec_len - length($s)) . $s;
}

# used only for formatting:
my $lat_prec_len = my $lat_full_len = lat_str(0);
my $lon_prec_len = my $lon_full_len = lon_str(0);
my $ele_prec_len = my $ele_full_len = length((sprintf ELE_PRECISION0, 0));
my $tim_len = length(epoch_to_str(0));
sub lat_str { return (sprintf LAT_PRECISION, $_[0]) =~ s/^\s*-(0(\.0*)?)$/$1/r; }
sub lon_str { return (sprintf LON_PRECISION, $_[0]) =~ s/^\s*-(0(\.0*)?)$/$1/r; }
sub lat_string {
    my $s = lat_str($_[0]);
    return " " x max(0, $lat_prec_len - length($s)) . $s;
}
sub lon_string {
    my $s = lon_str($_[0]);
    return " " x max(0, $lon_prec_len - length($s)) . $s;
}


my $bride_gap_insert_txt    = "bridging time gap: inserted point";
my $bride_gap_duplicate_txt = "bridging time gap: duplicated point ";

sub parse_point {
    my $str = $_[0];
    my $is_trkpt = $str =~ m/^<trkpt /;

    my $lat, my $lon;
    abort("FATAL: cannot find lat/lon part of point: $_[0]")
        unless $str =~ s#^(<\w+)(\s+(\w+\s*=\s*"[^\"]*"\s*)+)\s*>#$1>#s;
    my $lat_lon_str = $2;
    if ($lat_lon_str =~ m#\slat\s*=\s*"(-?\d*(\.\d+)?)"\s*lon\s*=\s*"(-?\d*(\.\d+)?)"#s) {
        ($lat, $lon) = ($1, $3);
    } elsif ($lat_lon_str =~ m#\slon\s*=\s*"(-?\d*(\.\d+)?)"\s*lat\s*=\s*"(-?\d*(\.\d+)?)"#s) {
        ($lat, $lon) = ($3, $1);
    } else {
        abort("FATAL: error parsing lat/lon part '$lat_lon_str' of point: $_[0]");
    }
    $lat_full_len = max($lat_full_len, length($lat));
    $lon_full_len = max($lon_full_len, length($lon));

    my $ele_str = my $ele = "";
    if ($str =~ s#[ \t]*(<ele>\s*([^\s<]*)\s*</ele>)\n?##s) {
        ($ele_str, my $ele_inner) = ($1, $2);
        abort("FATAL: error parsing elevation value '$ele_inner' in point: $_[0]")
            unless $ele_inner =~ m/^\s*(-?\d*(\.\d+)?)(\.\d*)?\s*$/; # ignores any decimals past any second decimal dot, e.g, takes ".12" on ".12.34" input
        $ele = $1; # may be empty or 0
        $ele_prec_len = max($ele_prec_len, length((sprintf ELE_PRECISION0, $ele)));
        $ele_full_len = max($ele_full_len, length($ele));
    }

    my $tim = "", my $sec;
    if ($str =~ s#[ \t]*(<time>\s*([^\s<]*)\s*</time>)\n?##s) {
        $tim = my $orig_tim = $2;
        $tim =~ s/\.0+Z/Z/;
        $tim_len = max($tim_len, length($tim));
        $sec = str_to_epoch($tim);
        $sec = 0 if $orig_tim eq "NO" || $orig_tim eq "NO_TIME" || $orig_tim =~ m/^\s*$/;
        abort("FATAL: error parsing time value '$orig_tim' in point: $_[0]")
            unless defined $sec;
        $sec += 0.001 if $tim =~ m/\.999Z/; # for some reason, at least OruxMaps on Android often reports 0.001 seconds less than a full second
        if (TIME_CORRECTION) {
            $sec += TIME_CORRECTION;
            $tim = epoch_to_str($sec);
        } else {
            $tim =~ s/(T\d\d)Z/$1:00Z/;
            $tim =~ s/(T\d\d:\d\d)Z/$1:00Z/;
        }
    }
    $ele = "" if !$ele_str && $tim eq ""; # for routes generated, e.g., using Google My Maps

    my $cmt = "";
    if ($is_trkpt && !$prune_cmts && $str =~ s#[ \t]*<cmt>\s*(.*?)\s*?</cmt>\r?\n?##s) {
        $cmt = $1;
        $cmt =~ s#\n##sg;
    }

    my $src = "";
    if ($is_trkpt && $str =~ s#[ \t]*<src>\s*(.*?)\s*?</src>\r?\n?##s) {
        $src = $1;
        $src =~ s#\r?\n##sg; # TODO match \r? also in other similar situations
    }

    my $fix = "";
    if ($is_trkpt && $str =~ s#[ \t]*<fix>\s*(.*?)\s*?</fix>\r?\n?##s) {
        $fix = $1;
        $fix =~ s#\r?\n##sg;
    }

    my $ext = "";
    while ($str =~ s#[ \t]*<extensions>(.*?)</extensions>\n?#($ext .= $1, "")#es) {
        $ext =~ s#\n##sg;
    }

    my $spd = "";
    if ($is_trkpt && $ext =~ m#<(\w+:)?speed>([^<]*)</(\w+:)?speed>#s) {
        abort("FATAL: inconsistent tags '<$1speed>' and '</$3speed>' in speed extension of in point: $_[0]") if ($1 // "") ne ($3 // "");
        my ($prefix, $spd_str) = ($1, $2);
        abort("FATAL: error parsing speed extension value '$spd_str' in point: $_[0]")
            unless $spd_str =~ m/^\s*(-?\d*(\.\d+)?)(\.\d*)?\s*$/; # ignores any decimals past any second decimal dot, e.g, takes ".12" on ".12.34" input
        # we assume that recorded speed takes altitude into account (i.e., is not just horizontal/lateral speed),
        # see also https://forums.geocaching.com/GC/index.php?/topic/209525-gps-speed-3d-or-2d/
        $spd = $1;
        $spd *= 1 / MPS2KMPH unless $prefix; # usually, speed is given m/s, while our abbreviated extension uses km/h
        $spd_prec_len = max($spd_prec_len, length(spd_str($spd))); # used only for spd_string
    }

    my $rest = $str;
    $rest =~ s#\n##sg;
    $rest =~ s#>[ \t]+<#><#g;
    # clean up potential errors in input:
    $rest =~ s#<ele></ele>##;
    $rest =~ s#<time></time>##;

    my $ignore = 0;
    $ignore = 1 if
        (diff_defined($sec, $begin_sec) // 0) < 0 && ++$num_ignored_before ||
        (diff_defined($sec,   $end_sec) // 0) > 0 && ++$num_ignored_after  ||
         $max_elem  && $_[0] =~ m#<$max_elem>\s*(-?[\.\d]+)\s*</$max_elem>#s && $1 > $max_elem_limit && ++$num_elem_above ||
         $min_elem  && $_[0] =~ m#<$min_elem>\s*(-?[\.\d]+)\s*</$min_elem>#s && $1 < $min_elem_limit && ++$num_elem_below ||
        # silently ignore any point inserted or duplicated by earlier run of this tool:
        ($src =~ m#^^($bride_gap_insert_txt|$bride_gap_duplicate_txt(after|before))$# && $fix eq "none");

    return ($lat, $lon, $ele, $tim, $sec, $spd, $rest, $cmt, $src, $fix, $ext, $ignore);
}

sub plural { my ($n, $str) = @_; return "$n $str".($n == 1 ? "" : "s"); }
sub points { return plural ($_[0], "point"); }

sub info_trkpt {
    stacktrace() unless defined $_[1];
    print_trkpt(shift @_, "INFO   : " . shift @_, @_) if $output_info;
}

sub warn_trkpt {
    stacktrace() unless defined $_[1];
    print_trkpt(shift @_, "WARNING: " . shift @_, @_) if $output_warnings;
}

sub abort_trkpt {
    stacktrace() unless defined $_[1];
    print_trkpt(shift @_, "FATAL: " . shift @_, @_);
}

sub print_trkpt {
    stacktrace() unless defined $_[0];
    my ($i, $head, $tail, $lat, $lon, $ele, $tim, $sec) = @_;
    # $head may be empty, $tail may be empty or undef
    $tail = " ".$tail if $tail && !($tail =~ m/^;/);
    ($lat    , $lon    , $ele    , $tim    , $sec    ) =
    ($LAT[$i], $LON[$i], $ELE[$i], $TIM[$i], $SEC[$i]) if $i >= 0;
    stacktrace() if (!defined $lat || !defined $lon || !defined $ele || !defined $tim);

    # using original precision here for better traceability
    # $lat = lat_str($lat);
    # $lon = lon_str($lon);
    my $elev = $ele;
    if ($elev eq "") {
        $elev         = "NONE";
        $ele_prec_len = max($ele_prec_len, length($elev));
        $ele_full_len = max($ele_full_len, length($elev));
    } elsif ($ele_full_len > $ele_prec_len) {
        my $decimals = -1;
        $decimals = length($1) if $elev =~ m/\.(\d+)/;
        my $n_spaces = $ele_full_len - $ele_prec_len - $decimals - 1;
        $elev .= ' ' x $n_spaces if $n_spaces > 0;
    }
    my $time = $tim eq "NO" ? "" : " ".($tim eq "" ? "NO_TIME" : $tim);
    print_line("$head "
               . "<trkpt"
               . ' lat="' . ($lat . '"' . (' ' x max(0, $lat_full_len - length($lat))))
               . ' lon="' . ($lon . '"' . (' ' x max(0, $lon_full_len - length($lon))))
               . " <ele>" . ele_string(-$ele_full_len, $elev)
               . $time
               .($tail ? (" " x max(0, $tim_len + length(" ") - length($time))).$tail : ""));
    abort("\n\n\n") if $head =~ "^FATAL:";

    $head =~ s|^(INFO)\s*:|$1:|;
    $head =~ s/(.*):\s*$/$1/ ||
    $head =~ s/the( most recently merged)$/this$1 point/ ||
    $head =~ s/(skipping) (regular)$/$1 this $2 point/ ||
    $head =~ m/(start|end)\s*$/ ||
    $head =~ s/ at$// ||
    $head =~ s/$/ this point/;
    $tail = "" unless $tail; # enforce output also in case $lean_stat_wpts
    add_stat_wpt($i, $head, $tail, undef, $lat, $lon, $ele, $tim, $sec) if $head =~ m|^INFO:|    && $info_wpts;
    add_stat_wpt($i, $head, $tail, undef, $lat, $lon, $ele, $tim, $sec) if $head =~ m|^WARNING:| && $warn_wpts;
}

sub diff_defined  { return defined $_[0] && defined $_[1] ? $_[0] - $_[1] : undef; }
sub diff_nonempty { stacktrace() unless defined $_[0] && defined $_[1];
                    return $_[0] ne ""   && $_[1] ne ""   ? $_[0] - $_[1] : undef; }
sub eq_nonempty   { return $_[0] eq "" ? $_[1] eq ""
                        :  $_[1] eq "" ? $_[0] eq "" : $_[0] == $_[1] }

sub section { return ($FLG[$_[0]] & (START_MERGE | PAUSE_MERGE)) != 0 } # switched between regular and merged track

sub rate { return defined $_[0] && $_[1] ? H2S * $_[0] / $_[1] : ""; } # ascent/descent rate in m/h

sub recalculate_DIS {
    my ($i, $const_2d) = @_;
    my $prev_i = PREV($i);
    return if !defined $prev_i || $SEG[$i];
    $DIS2[$i] = sq_distance_2d($LAT[     $i], $LON[     $i],
                               $LAT[$prev_i], $LON[$prev_i]) unless $const_2d;
    $DIS[$i] = distance_3d($DIS2[$i], $ELE[$i], $ELE[$prev_i]);
}
sub recalculate_SPD {
    my $i = $_[0];
    return if !defined PREV($i) || $SEG[$i] || section($i) || !($FLG[$i] & COMPUTED_SPEED);
    my $diff_time = diff_defined($SEC[$i], $SEC[PREV($i)]);
    $SPD[$i] = $DIS[$i] / $diff_time if $diff_time;
}
sub recalculate_all {
    my ($i, $old_i) = @_;
    $DIR[$i] = undef; # invalidate cached direction
    my $const_2d = defined $old_i &&
        ($old_i == $i || $LAT[$old_i] == $LAT[$i] && $LON[$old_i] == $LON[$i]);
    recalculate_DIS($i, $const_2d);
    recalculate_SPD($i);
}

#my @COS2; # approximation of cos(2phi) for efficiency
#for (my $phi2 = 0; $phi2 <= 180; $phi2++) {
#    $COS2[-$phi2] = $COS2[$phi2] = abs(cos(deg2rad(($phi2 + 1) / 2)));
#    # error in cos(deg2rad(($phi2 / 2))) - $COS2[$phi2]
#    # max absolute error is < 0.0088
#    # max relative error is < 0.027  for latitudes below North Cape
#}
# https://stackoverflow.com/questions/639695/how-to-convert-latitude-or-longitude-to-meters
my @COS;
sub sq_distance_2d {
    my ($lat, $lon, $lat2, $lon2) = @_;
    my $diff_lat = ($lat - $lat2) * METERS_PER_DEGREE_LAT;
    return $diff_lat * $diff_lat if $lon == $lon2; # optimization in special case, no need to use cos()
    # with the following two approximations, cache hit rate is around 90% and resulting error is rarely more than a meter:
    my $mean_lat = ($lat + $lat2) / 2;
    # cosine cache input resolution 0.0001 deg corresponding to < 10 meters:
    my $rounded_scaled_abs_mean_lat = ($mean_lat < 0 ? -$mean_lat : $mean_lat) * 10000 + 0.5;
    my $diff_lon = ($lon - $lon2) * METERS_PER_DEGREE_LAT *
        # $COS2[$lat + $lat2]; # optimization would provide almost no efficiency advantage
        ($COS[$rounded_scaled_abs_mean_lat] // ($COS[$rounded_scaled_abs_mean_lat] = cos(deg2rad($mean_lat))));
    # my $diff_lon = ($lon * cos(deg2rad($lat)) - $prev_lon * cos(deg2rad($prev_lat))) * METERS_PER_DEGREE_LAT;
    return $diff_lat * $diff_lat + $diff_lon * $diff_lon;
}
my @SQR;
sub sqrt_with_cache {
    my $sq_dis = $_[0];
    # with the following, cache hit rate is around 98% and resulting error is rarely more than a meter:
    return sqrt($sq_dis) if $sq_dis >= 10000; #  do not cache for distances >= 100 m
    my $rounded_sq_dis = $sq_dis + 0.5;
    return $SQR[$rounded_sq_dis] // ($SQR[$rounded_sq_dis] = sqrt($sq_dis));
}
sub distance_3d {
    my ($sq_dis2d, $ele, $ele2) = @_;
    my $diff_ele = diff_nonempty($ele, $ele2) // 0; # assuming no altitude change if no altitude available
    return sqrt_with_cache($sq_dis2d + $diff_ele * $diff_ele);
}
sub distance {
    my ($lat, $lon, $ele, $lat2, $lon2, $ele2) = @_;
    stacktrace()
        if !defined $lat  || !defined $lon  || !defined $ele
        || !defined $lat2 || !defined $lon2 || !defined $ele2;
    return distance_3d(sq_distance_2d($lat, $lon, $lat2, $lon2), $ele, $ele2);
#
#$distance = Math::Trig::great_circle_distance( #does not account for $diff_ele!
#						  deg2rad($lon) , deg2rad(90 - $lat ),
#						  deg2rad($lon2), deg2rad(90 - $lat2),
#						  40*1000*1000/pi/2); #http://perldoc.perl.org/Math/Trig.html
## debug "diff_lat=$diff_lat, diff_lon=$diff_lon, dis=$dis, distance=$distance";
}

sub comp_diffs {
    my $i = shift; # calculate differences relative to previous point PREV($i)
    # $prev_comp_spd is assumed to be average speed computed between previous points
    my $prev_i = PREV($i);
    my ($prev_ele, $prev_sec, $prev_dis, $prev_comp_spd) = ($ELE[$prev_i], $SEC[$prev_i], $DIS[$prev_i], $SPD[$prev_i]);
    my (     $ele,      $sec,      $dis,           $spd) = ($ELE[     $i], $SEC[     $i], $DIS[     $i], $SPD[     $i]);
    my $diff_ele = diff_nonempty($ele, $prev_ele);

    #http://forums.howwhatwhy.com/showflat.php?Cat=&Board=scigen&Number=-208125
    my $timediff = diff_defined($sec, $prev_sec);
    my $rate= rate($diff_ele, $timediff);
    my $comp_spd = defined $timediff ? ($timediff > 0 ? $dis / $timediff : $prev_comp_spd) : "";
    $spd_prec_len = max($spd_prec_len, length(spd_str($comp_spd))) if $comp_spd ne ""; # used only for spd_string
    my $acc = $comp_spd ne "" && $prev_comp_spd ne "" && $timediff ? ($comp_spd - $prev_comp_spd) / $timediff : "";
    return ($diff_ele // "", $timediff // "", $rate, $dis, $comp_spd, $acc);
}

my @SLO; # cached slope; no need to use quick approximation (such as quotient of elevation difference and 3d distance)
sub slope {
    my ($prev_i, $i) = @_; # both are assumed to be defined
    my $value = $SLO[$i];
    debug_log "$i $TIM[$i]     slope=$value (cached)" if defined $value;
    return $value if defined $value;
    return undef if $ELE[$i] eq "" || $ELE[$prev_i] eq "";
    my $ele_diff = $ELE[$i] - $ELE[$prev_i];
    my $sq_dis_2d = $DIS2[$i]; # not taking into account: $ele_avg = ($ELE[$i] + $ELE[$prev_i]) / 2
    $value = $ele_diff == 0 ? 0
        : $sq_dis_2d == 0 ? ($ele_diff < 0 ? -90 : 90)
        : atan2($ele_diff, sqrt_with_cache($sq_dis_2d)) * RAD2DEG;
    debug_log "$i $TIM[$i]     slope=$value";
    return ($SLO[$i] = $value);
}

use constant N_AVG_SLOPE => 4;
sub avg_slope { # between the last N_AVG_SLOPE points, assuming 0 where not available
    my $sum = 0;
    my $i = $_[0];
    my $prev_i = PREV($i);
    my $n;
    for ($n = 0; $n < N_AVG_SLOPE; $n++) {
        last if $SEG[$i] || section($i);
        my $slope = slope($prev_i, $i) // 0;
        $sum += $slope;
        ($i, $prev_i) = ($prev_i, PREV($prev_i));
    }
    # if ($n == 1) { # use as fallback the average with slope of next point, if available
    #     my $i2 = NEXT($_[0]);
    #     my $slope2 = defined $i2 ? slope($i2) : undef;
    #     return ($sum + $slope2) / 2 if defined $slope2;
    # }
    debug_log "$_[0] $TIM[$_[0]] avg_slope=".($n == 0 ? "undef" : $sum / N_AVG_SLOPE);
    return $n == 0 ? undef : $sum / N_AVG_SLOPE
}

sub get_direction { # compute and store in @DIR the direction from PREV($i) to $i, and return it
    my $i = $_[0];
    return $DIR[$i] if defined $DIR[$i]; # return cached value
    my $prev_i = PREV($i);
    return undef unless defined $prev_i;
    my $curr_lat_diff = $LAT[$i] - $LAT[$prev_i];
    my $curr_lon_diff = $LON[$i] - $LON[$prev_i];
    return undef unless $curr_lat_diff || $curr_lon_diff;
    return $DIR[$i] = atan2($curr_lon_diff, $curr_lat_diff);
}

sub comp_angle { # calculate horizontal direction change (i.e., turning angle)
    # between last, current, and next point, or undef if not available
    # provides almost no efficiency penalty
    my ($i, $timediff, $dis) = @_;
    my ($theta_diff, $spd2) = (undef, "");

    my $next_i = NEXT($i);
    if (defined $next_i && !$SEG[$next_i] && !section($i) && !section($next_i)) { # next point exists and is not at segment start and not at merge start/pause
        my $next_dis      =              $DIS[$next_i];
        my $next_timediff = diff_defined($SEC[$next_i], $SEC[$i]);
        my $timediff2 = $timediff eq "" || !defined $next_timediff  ? 0 : $timediff + $next_timediff;
        my      $dis2 = $dis      eq "" ||          $next_dis eq "" ? 0 : $dis      + $next_dis;
        $spd2 = $timediff2 <= 0 || $dis2 <= 0 ? "" : $dis2 / $timediff2; # average speed
        if ($spd2 ne "" && MPS2KMPH * $spd2 > MAX_PLAUSIBLE_ANGLE_SPD_THRESHOLD) {
            my ($theta1, $theta2) = (get_direction($i), get_direction($next_i));
            if (defined $theta1 && defined $theta2) {
                # debug (($theta1 * RAD2DEG)." <$i> ".($theta2* RAD2DEG));
                $theta_diff = ($theta2 - $theta1) * RAD2DEG % 360;
                $theta_diff -= 360 if $theta_diff >= 180;
            }
        }
    }
    return ($theta_diff, $spd2);
}


### parse command line
# TODO clean up use of ARGV

for (my $i = 0; $#ARGV - $i >= 0; ) {
    my $opt  = $ARGV[$i];
    my $arg1 = $ARGV[$i + 1];
    my $arg2 = $ARGV[$i + 2];

    # not checking for duplicate options

    if ($opt eq "-h" || $opt =~ /^-*help$/) {
        print_line($usage);
        exit 0;
    } elsif ($opt eq "-version") {
        print_line($TOOL_ID);
        exit 0;
    } elsif ($opt eq "-swim"  ||
        $opt eq "-walk"  ||
        $opt eq "-cycle" ||
        $opt eq "-drive" ||
        $opt eq "-fly") {
        $activity = $opt;
    } elsif ($opt eq "-no_sanitize") {
        $sanitize = !SANITIZE_TRKPTS;
    } elsif ($opt eq "-no_insert") {
        $insert_points = 0;
    } elsif ($opt eq "-smooth") {
        $smoothing = 1;
    } elsif ($opt eq "-phases") {
        $phases = 1;
    } elsif ($opt eq "-merge") {
        abort("missing file argument for -$opt option") if $#ARGV - $i < 1;
        $merge = $arg1;
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-segs" || $opt eq "-days") {
        if ($opt eq "-segs") {
            abort("cannot use both -days and -$opt option") if $days;
            $segs = 1;
        } else {
            abort("cannot use both -segs and -$opt option") if $segs;
            $days = 1;
        }
        if ($#ARGV - $i >= 1 && $arg1 =~ m/^\d/) {
            abort("argument '$arg1' of -$opt option is not of the form <n>..[<m>][(+|-)<d>] with natural numbers n, m, and integer d")
            unless $arg1 =~ m/^(\d+)\.\.(\d+)?([-+]\d+)?+$/ && $1 > 0;
            ($part_start, $part_end, $part_offset) = ($1, $2, $3); # $part_end < $part_start will lead to empty selection, just like without $opt
            splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
        }
    } elsif ($opt eq "-tz") {
        abort("missing value argument for -$opt option") if $#ARGV - $i < 1;
        abort("cannot parse value argument in '$opt $arg1'") unless looks_like_number($arg1);
        $tz = $arg1 + 0;
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-split") {
        abort("missing base file name argument for -$opt option") if $#ARGV - $i < 1;
        $split = $arg1;
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-lean_stat_wpts") {
        $lean_stat_wpts = 1;
    } elsif ($opt eq "-info_wpts") {
        $info_wpts = 1;
    } elsif ($opt eq "-warn_wpts") {
        $warn_wpts = 1;
    } elsif ($opt eq "-cmt_wpts") {
        abort("missing pattern argument for -$opt option") if $#ARGV - $i < 1;
        $cmt_wpts = $arg1;
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-src_wpts") {
        abort("missing pattern argument for -$opt option") if $#ARGV - $i < 1;
        $src_wpts = $arg1;
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-weight") {
        abort("missing value argument for -$opt option") if $#ARGV - $i < 1;
        abort("cannot parse value argument in '$opt $arg1'") unless looks_like_number($arg1);
        $weight = $arg1 + 0;
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-begin" || $opt eq "-end") {
        abort("missing time argument for $opt option") if $#ARGV - $i < 1;
        my $sec = str_to_epoch($arg1);
        abort("cannot parse time argument in '$opt $arg1'") unless defined $sec;
        $begin_sec = $sec if $opt eq "-begin";
        $end_sec   = $sec if $opt eq "-end";
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-resolution") {
        abort("missing number of seconds argument for $opt option") if $#ARGV - $i < 1;
        abort("illegal natural number argument in '$opt $arg1'") unless $arg1 =~ m/^\d+$/;
        $min_input_timediff = $arg1 + 0;
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-max" || $opt eq "-min") {
        abort("missing element name argument for $opt option") if $#ARGV - $i < 1;
        abort("missing limit argument for -$opt option") if $#ARGV - $i < 2;
        abort("cannot parse limit argument in '$opt $arg1 $arg2'") unless looks_like_number($arg2);
        ($max_elem, $max_elem_limit) = ($arg1, $arg2 + 0) if $opt eq "-max";;
        ($min_elem, $min_elem_limit) = ($arg1, $arg2 + 0) if $opt eq "-min";;
        splice @ARGV, $i + 1, 2; # remove from ARGV the option arguments
    } elsif ($opt eq "-max_ext" || $opt eq "-min_ext") {
        abort("missing extension name argument for $opt option") if $#ARGV - $i < 1;
        abort("missing limit argument for -$opt option") if $#ARGV - $i < 2;
        abort("cannt parse limit argument in '$opt $arg1 $arg2'") unless looks_like_number($arg2);
        ($max_ext, $max_ext_limit) = ($arg1, $arg2 + 0) if $opt eq "-max_ext";
        ($min_ext, $min_ext_limit) = ($arg1, $arg2 + 0) if $opt eq "-min_ext";
        splice @ARGV, $i + 1, 2; # remove from ARGV the option arguments
    } elsif ($opt eq "-min_dist") {
        abort("missing distance argument for -$opt option") if $#ARGV - $i < 1;
        $min_dist = $arg1;
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-min_angle") {
        abort("missing angle argument for -$opt option") if $#ARGV - $i < 1;
        $min_angle = $arg1;
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-prune_wpts") {
        $prune_wpts = 1;
    } elsif ($opt eq "-prune_cmts") {
        $prune_cmts = 1;
    } elsif ($opt eq "-prune_exts") {
        $prune_exts = 1;
    } elsif ($opt eq "-ni") {
        $output_info = 0;
    } elsif ($opt eq "-nw") {
        $output_warnings = 0;
    } elsif ($opt eq "-o") {

     # not checking for operlap with shell-level redirection of STDOUT using '>'
        abort("missing outfile argument for -o option") if $#ARGV - $i < 1;
        open($out, "> $arg1") || abort("FATAL: cannot open output file $arg1: $.");
        splice @ARGV, $i + 1, 1; # remove from ARGV the option argument
    } elsif ($opt eq "-debug") {
        $debug = 1;
    } elsif ($opt =~ m/^-/) {
        abort("unknown option: $opt");
    } else {
        $i++;
        next; # infile
    }
    splice @ARGV, $i, 1; # remove from ARGV the option just handled
}
if ($prune_exts) {
    warning("-min_ext option has little effect since -prune_exts is given") if $min_ext;
    warning("-max_ext option has little effect since -prune_exts is given") if $max_ext;
}
abort "cannot smoothen tracks without sanitizing" if $smoothing && !$sanitize; # TODO maybe fix this limitation

$part_start  = 1 unless defined $part_start;
$part_offset = 0 unless defined $part_offset;
my $all = $part_start - 1;


if ($activity eq "-swim") {
    $MIN_TIMEDIFF           = 4;
    $MAX_PLAUSIBLE_SPEED    = 7;
    $MAX_PLAUSIBLE_ACCEL    = 1;
    $MAX_PLAUSIBLE_DECEL    = 2;
    $MAX_PLAUSIBLE_ELE_GAIN = 1;
    $MAX_PLAUSIBLE_ELE_LOSS = 1;
} elsif ($activity eq "-walk") {
    $MIN_TIMEDIFF           = 4;
    $MAX_PLAUSIBLE_SPEED    = 20;
    $MAX_PLAUSIBLE_ACCEL    = 0.5;
    $MAX_PLAUSIBLE_DECEL    = 1;
    $MAX_PLAUSIBLE_ELE_GAIN = 1800;
    $MAX_PLAUSIBLE_ELE_LOSS = 3600;
} elsif ($activity eq "-cycle") {
    $MIN_TIMEDIFF           = 2;
    $MAX_PLAUSIBLE_SPEED    = 80;
    $MAX_PLAUSIBLE_ACCEL    = 4;
    $MAX_PLAUSIBLE_DECEL    = 8;
    $MAX_PLAUSIBLE_ELE_GAIN = 3000;
    $MAX_PLAUSIBLE_ELE_LOSS = 9000;
} elsif ($activity eq "-drive") { # default
    $MIN_TIMEDIFF           = 2;     # in seconds
    $MAX_PLAUSIBLE_SPEED    = 200;   # maximal speed in km/h
    $MAX_PLAUSIBLE_ACCEL    = 5;     # maximal acceleration m/s/s
    $MAX_PLAUSIBLE_DECEL    = 10;    # maximal deceleration m/s/s
    $MAX_PLAUSIBLE_ELE_GAIN = 20000; # maximal ascent rate in m/h
    $MAX_PLAUSIBLE_ELE_LOSS = 20000; # maximal descent rate in m/h
    $PHASES_REPORT_THRESHOLD = 100;
} elsif ($activity eq "-fly") {
    $MIN_TIMEDIFF           = 1;
    $MAX_PLAUSIBLE_SPEED    = 1200;
    $MAX_PLAUSIBLE_ACCEL    = 10;
    $MAX_PLAUSIBLE_DECEL    = 10;
    $MAX_PLAUSIBLE_ELE_GAIN = 30000;
    $MAX_PLAUSIBLE_ELE_LOSS = 50000;
    $PHASES_REPORT_THRESHOLD = 100;
} else {
    abort("unknown activity: $activity\naborting");
}
my $INSERTION_MAX_SPEED = $MAX_PLAUSIBLE_SPEED / 2;
my $MAX_SPEED_DEVIATION = $MAX_PLAUSIBLE_SPEED * MAX_SPEED_DEVIATION;


### main loops

## read all trackpoints from all track segments, optionally completing them

sub read_tracks
{
    our ($lat, $lon, $ele, $tim, $sec, $spd);
    sub info_trkpt_curr  { info_trkpt (-1, $_[0], $_[1], $lat, $lon, $ele, $tim, $sec); }
    sub warn_trkpt_curr  { warn_trkpt (-1, $_[0], $_[1], $lat, $lon, $ele, $tim, $sec); }
    sub print_trkpt_curr { print_trkpt(-1, $_[0], $_[1], $lat, $lon, $ele, $tim, $sec); }
    sub distance_curr {
        my $sq_dis2d = sq_distance_2d($lat, $lon, $_[0], $_[1]);
        return ($sq_dis2d, distance_3d($sq_dis2d, $ele, $_[2]));
    }

    my $merge = $_[0];
    my $i_start = $#TIM + 1;
    my $ele_miss_start  = $i_start; # first index of recent stretch with missing ele ; TODO restrict scope
    my $tim_miss_start  = $i_start; # first index of recent stretch with missing time; TODO restrict scope
    my $num_no_ele      = 0;
    my $ele_gap_length  = 0;
    my $num_no_time     = 0;
    my $time_gap_length = 0;

    $/ = "<trkpt ";
    open M, $merge || abort("FATAL: cannot open alternative GPX source file $merge: $.") if $merge;

    my $state = 0;
    my ($prev_lat, $prev_lon, $prev_ele, $prev_sec, $prev_dis);
    my $prev_tim = "";
    my ($prev_avail_tim, $prev_avail_sec);
    my $prev_avail_ele = "";
    (my $ele_miss_length = 0, my $tim_miss_length = 0);
    my $i = $i_start;

    { local $SIG{"__WARN__"} = sub { abort $_[0] }; # https://www.perlmonks.org/?node_id=11149744 Diamond errors
    while ($_ = ($merge ? <M> : <>)) {
        # anywhere (in headers, tracks, waypoints, trackpoints, ...):
        push @NAME, $2 if m|<metadata>\s*<name>(<!\[CDATA\[)?(.*?)(\]\]>)?</name>|s;
        $num_pruned_cmts  += s#[ \t]*<cmts>.*?</cmt>\n?##sg              if $prune_cmts;
        $num_pruned_descs += s#[ \t]*<desc>.+?</desc>\n?##sg             if PRUNE_DESCS;
        $num_pruned_links += s#[ \t]*<link.*?>.+?</link>\n?##sg          if PRUNE_LINKS;
        $num_pruned_exts  += s#[ \t]*<extensions>.*?</extensions>\n?##sg if $prune_exts;
        if (!$prune_exts && ($min_ext || $max_ext)) {
            s@(<extensions>)(.*?)(</extensions>\n?)@ my ($start, $ext, $end) = ($1, $2, $3);
                $ext =~ s#(<$min_ext>\s*(-?[\.\d]+)\s*</$min_ext>)#$2 < $min_ext_limit ? ($num_ext_below++, "") : $1#esg if $min_ext;
                $ext =~ s#(<$max_ext>\s*(-?[\.\d]+)\s*</$max_ext>)#$2 > $max_ext_limit ? ($num_ext_above++, "") : $1#esg if $max_ext;
                $ext eq "" ? "" : "$start$ext$end"
            @esg;
        }

        #collect waypoints, possibly in later regular input headers
        my $prev_wpt_time_stripped = "";
        while (s#(<wpt .+?)(</wpt>)\n?##s)
        {
            my $wpt = $1 . $2;

            # workaround for OruxMaps partly duplicating waypoints of <type>Finishing Point</type> when manually starting new segment
            my $wpt_time_stripped = $wpt;
            $wpt_time_stripped =~ s#<time>.*</time>##;
            next if $wpt =~ m#<type>Finishing Point</type># && $wpt_time_stripped eq $prev_wpt_time_stripped;
            $prev_wpt_time_stripped = $wpt_time_stripped;

            ($lat, $lon, $ele, $tim, my $sec, my $spd_unused, my $rest, my $cmt_unused,
             my $src_unused, my $fix_unused, my $exts, my $ignore) =  parse_point($wpt); # collect some basic data even if pruned
            $rest =~ s#^<wpt[^>]*>##; # ignoring any further wpt attributes, yet there should not be any (besides lat and lon)
            $rest =~ s#</wpt>$##;
            ++$num_pruned_wpts if $prune_wpts && !$ignore;
            if (!$ignore && !$prune_wpts
                && !(PRUNE_AUTOMATIC_WPTS &&
                     $rest =~ m#<name><!\[CDATA\[((day )?\d+)?\]\]></name># && ++$num_auto_wpts)
                # better not leave out waypoints added previously by this tool because this would disturb the ordering and duplicates are now taken care of by add_stat_point():
                # && !($rest =~ m#\[(start|end.*?)\]#) # start/end waypoint added previously by this tool
                # && # waypoint with statistics added previously by this tool:
                # !($rest =~ m#\[((max|min)\s+(altitude|elevation|height)|max\s+speed|max\s+(climb|ascent|descent)\s+rate|total\s+(
                #  ascent|descent|gain|loss))\s+=\s+[+-]?[\.\d]+ k?m(\/h)?( at -?[\.\d]+ km/h|;\s+avg\s+rate\s+=\s+[+-]?[\d.]+ m/h)?\]#)
                )
            {
                push @WLAT, $lat;
                push @WLON, $lon;
                push @WELE, $ele;
                push @WTIM, $tim;
                push @WSEC, $sec; # may be 0
                push @WTXT, $rest . ($exts ? "<extensions>$exts</extensions>" : "");
                $wpt =~ s/(^|\n)/\n  /sg; # indent waypoint lines, just for better readability of diagnostics
                push @WSTR, "'$wpt\n'";
                push @WIDX, -1;
            }
        }

        sub warn_smoothing {
            warning("smoothing has already been done on input")
                if $smoothing && $_[0] =~ m#(\d+) points .*?smoothened($|\W)# && $1 > 0;
        }
        if ($state == 0) {
            # processing first header (of first input file)

            warn_smoothing($1) if s#\r?\nGPXConv statistics (.*?)GPXConv statistics end\r?\n##s; # remove any earlier statistics section
            if (m#$/#s) {
                s#$/##s;    #remove trailing "<trkpt "
            } else {
                s#(^.*</type>).*#$1\n<trkseg>\n#s; #ignore rest of (header-only) file having no trkpt
            }
            s#[ \t]+(\r?\n)$#$1#s; # remove any trailing spaces on line before
            s#[ \t]*(<trkseg>)#$1#s;
            $HEAD  = $_ unless $merge; # ignore header of alternative input
            $state = 1; # expecting the very first trkpt
        } elsif (m#(<gpx[ >])#s) { # processing further GPX header, typically from next file
            ++$num_ignored_hdrs;
            warn_smoothing($1) if m#\r?\nGPXConv statistics (.*?smoothened)(.*?)GPXConv statistics end\r?\n#;
            if (!$merge && $#TIM >= 0) { # typically, new file
                info_trkpt($i - 1, "end of  segment at") if $i != $i_start;
                info("further track, from file $ARGV");
            }
            $state = 2 if $state > 2; # expecting new trkseg
        } else {
            # when $state == 1 || $state == 2: processing first trkpt in new trkseg
            # when $state == 3: processing further trkpt in trkseg
            my ($sec, $spd, $cmt, $src, $fix, $exts);

            m#(.*?</trkpt>)#s;
            my $trkpt = $1;
            abort("FATAL: cannot find end of trackpoint: <trkpt $_") unless $trkpt;
            ($lat, $lon, $ele, $tim, $sec, $spd, my $rest_unused, $cmt, $src, $fix, $exts, my $ignore) = parse_point("<trkpt $trkpt");
            my $seg = $state < 3; # first trkpt of a segment
            goto continuation if $ignore
                || !($ARGV =~ m/test(\.out\d?|_sanitization\.gpx)$/) && !$seg
                && ((diff_defined($sec, $prev_sec) // $min_input_timediff) < $min_input_timediff
                    || $lat == $prev_lat && $lon == $prev_lon && $ele eq $prev_ele && $tim eq $prev_tim)
                && ++$num_ignored_diff; # silently ignore
            ## debug_log "read non-ignored point $i from with lat=$lat lon=$lon ele=$ele tim=$tim sec=".($sec // "none")." spd=".spd_or_none($spd)."; state = $state";

            my $diff_time = diff_defined($sec, $prev_avail_sec) // 0;

            if ($state > 1) {
                my $severity = $seg ? "WARNING" : "FATAL"; # at segment start, just warn
                print_trkpt_curr("$severity: before", "negative time difference (-".timediff_string(-$diff_time).")")
                    if $diff_time < 0; # this implies defined $sec && defined $prev_avail_sec
                if (defined $sec && defined $prev_avail_sec && $tim eq $prev_avail_tim) { # really equal (before rounding seconds)
                    print_trkpt_curr("$severity: point before      ",
                                     "has same time but different position: $prev_lat, $prev_lon")
                        if $lat != $prev_lat || $lon != $prev_lon;
                    print_trkpt_curr("WARNING: point before      ",
                                     "has same time but different altitude: $prev_ele m")
                        if $ele != $prev_ele;
                }
            }

            my $flags = 0;
            $flags |= START_MERGE if $cmt =~ m/start alternative/;
            $flags |= PAUSE_MERGE if $cmt =~ m/pause alternative/;
            my ($sq_dis2d, $dis) = $i == $i_start || $state == 1 ? (0, 0) : distance_curr($prev_lat, $prev_lon, $prev_ele);
            if ($seg) { # at the first point in trkseg, first handle last point of segment before
                $prev_avail_ele = "";
                if (CARRY_OVER_ELE && $ele ne "") {
                    if ($i != $i_start && $ELE[$i - 1] eq "") {
                        # fill in altitude in last point of segment before
                        $FLG[$i - 1] |= SUBSTITUTED_ELE;
                        $ELE[$i - 1]  = $prev_ele = $ele;
                        recalculate_all($i - 1, $i - 1); # 2d part did not change
                        $dis = distance_3d($sq_dis2d, $ele, $prev_ele); # recalculation of 3d part
                        $ele_miss_length += $DIS[$i - 1] if $ele_miss_start != $i_start;
                        warn_trkpt($i - 1, "at segment end    ",
                                   "no altitude was given; using value from begin of next segment");
                    }
                } elsif ($ele_miss_start != $i_start) {
                    warn_trkpt($i - 1, "             up to",
                               "was unable to interpolate missing altitude within segment");
                    $ele_miss_start = $i_start;
                }
                if (CARRY_OVER_TIME && $tim ne "") {
                    if ($i != $i_start && $TIM[$i - 1] eq "") {
                        # fill in time in last point of segment before
                        $FLG[$i - 1] |= SUBSTITUTED_TIM;
                        $TIM[$i - 1]  = $tim;
                        $SEC[$i - 1]  = $sec;
                        recalculate_SPD($i - 1);
                        $tim_miss_length += $DIS[$i - 1] - $prev_dis if $tim_miss_start != $i_start;
                        warn_trkpt($i - 1, "at segment end    ",
                                   "no time     was given; using value from begin of next segment");
                    }
                } elsif ($tim_miss_start != $i_start) {
                    warn_trkpt($i - 1, "             up to",
                               "was unable to interpolate missing time within segment");
                    $tim_miss_start = $i_start;
                }
            }

            if ($state == 1) {
                my $name = $ARGV =~ s/^(test.*?\.out)\d$/$1/r; # avoid variation in test log output due to iteration count in file name
                info("initial track".($merge ? " from alternative input $merge" : ", from file $name"));
                info_trkpt_curr("first   segment at");
            } elsif ($state == 2) {
                info_trkpt_curr("further segment at", "after a time gap of ".timediff_string($diff_time)." and ".dis_string($dis));
            }

            if ($seg) { # now handle first point in trkseg
                my $resolution;
                if ($ele eq "") {
                    if ($i == $i_start || $prev_ele eq "") { # $prev_avail_ele
                        $resolution = " and no previous value available";
                    } elsif (CARRY_OVER_ELE) {
                        $ele = $prev_ele; # $prev_avail_ele; # will become $ELE[$i]
                        $resolution = "; using previous available value";
                        $flags |= SUBSTITUTED_ELE;
                    } else {
                        $resolution = " and carrying over previous available value is disabled";
                    }
                    warn_trkpt_curr("at segment start  ",
                                    "no altitude ".($resolution =~ / using / ? "was " : "")."given$resolution");
                }
                if ($tim eq "") {
                    if ($prev_tim eq "") {
                        $resolution = " and no previous value available";
                    } elsif (CARRY_OVER_TIME) {
                        $tim = $prev_tim; # will become $TIM[$i]
                        $sec = $prev_sec; # will become $SEC[$i]
                        $flags |= SUBSTITUTED_TIM;
                        $resolution = "; using previous value";
                    } else {
                        $resolution = " and carrying over previous value is disabled";
                    }
                    warn_trkpt_curr("at segment start  ",
                                    "no time     ".($resolution =~ / using / ? "was " : "")."given$resolution");
                }
            } else { # state == 3, further (i.e., not the first) point in trkseg
                if ($ele eq "" && FILL_MISSING_ELE && ($ele_miss_start != $i_start || $prev_ele ne "")) { # $prev_avail_ele
                    ($ele_miss_start, $ele_miss_length) = ($i, 0)
                      if $ele_miss_start == $i_start; # $i is != $i_start
                    $ele_miss_length += $dis;
                }

                # substitute missing time if posssible using previous time, recorded speed, and distance from previous point
                if (!defined $sec && defined $prev_sec && $spd ne "" && $dis > 0 && $spd > 0) {
                    my $comp_sec = $prev_sec + $dis / $spd;
                    $tim = epoch_to_str($comp_sec); # will become $TIM[$i]
                    $sec = $comp_sec;               # will become $SEC[$i]
                    warn_trkpt_curr("at                ", " no time recognized, " .
                                    "filled in using time difference calculated by distance / speed");
                    $FLG[$i] |= SUBSTITUTED_TIM;
                } # any inconsistency where either $dis is (nearly) 0 or $spd is (nearly) 0 are covered by warning w.r.t. MAX_SPEED_DEVIATION

                if ($tim eq "" && FILL_MISSING_TIME && ($tim_miss_start != $i_start ||
                                                        $prev_tim ne "")) # not $prev_avail_tim, which may be an earlier time
                {
                    ($tim_miss_start, $tim_miss_length) = ($i, 0) if $tim_miss_start == $i_start; # $i is != $i_start
                    $tim_miss_length += $dis;
                }
            }
            if ($ele ne "" && $ele_miss_start != $i_start) { # fill in missing altitude
                my $i_end = $i;
                if ($seg)
                { # have been carrying over ele to end of last section
                    $i_end--;
                } else {   # state == 3, further (i.e., not the first) point in trkseg
                    $ele_miss_length += $dis;
                }
                my $ele_curr = $ELE[$ele_miss_start - 1]; # $prev_avail_ele
                my $ele_diff = $ele - $ele_curr;
                # $ele_miss_length == 0 is special case of missing altitude over a stretch with all same lat/lon
                my $ele_gradient = $ele_miss_length > 0 ? $ele_diff / $ele_miss_length : $ele_diff / ($i - $ele_miss_start + 1);
                if ($ele_gradient != 0) {
                    for (my $j = $ele_miss_start; $j < $i_end; $j++) {
                        $ele_curr += ($ele_miss_length > 0 ? $DIS[$j] : 1) * $ele_gradient;
                        $ELE[$j]  = sprintf ELE_PRECISION3, $ele_curr;
                        $FLG[$j] |= INTERPOLATED_ELE;
                        $ele_full_len = max($ele_full_len, length($ELE[$j]));
                        my $prev_dis = $DIS[$j];
                        recalculate_all($j, $j); # 2d part did not change
                        $tim_miss_length += $DIS[$j] - $prev_dis if $tim_miss_start != $i_start && $TIM[$j] eq "";
                        # not changing $ele nor $prev_ele such that below time interpolation is not affected
                        my $Tim = $TIM[$j] ne "" ? $TIM[$j] : "NO_TIME";
                        debug_log "interpolating missing ele at point $j $Tim: new ele=$ELE[$j] dis=$DIS[$j]";
                    }
                    if ($i_end < $i) { # have been carrying over ele to end of last section
                        recalculate_all($i_end, $i_end); # 2d part did not change
                    } else {
                        $dis = distance_3d($sq_dis2d, $ele, $ELE[$i_end - 1]); # recalculation of 3d part
                    }
                }
                $ele_miss_start = $i_start;
            }
            if ($tim ne "" && $tim_miss_start != $i_start) { # fill in missing time by interpolation, assuming constant speed
                my $i_end = $i;
                if ($seg)
                { # have been carrying over time to end of last section
                    $i_end--;
                } else {   # state == 3, further (i.e., not the first) point in trkseg
                    $tim_miss_length += $dis;
                }
                my $sec_curr = my $sec_prev = $SEC[$tim_miss_start - 1]; # corresponds to $prev_tim
                my $sec_diff = $sec - $sec_curr;
                my $i_diff = $i_end - $tim_miss_start + 1;
                # $tim_miss_length == 0 is special case of missing time over a stretch with all same position;
                # since altitude is interpolated in same iteration, a stretch with all same lat/lon looks like all same position
                my $sec_gradient = $tim_miss_length > 0 ? $sec_diff / $tim_miss_length : $sec_diff / $i_diff;
                my $const_spd = $sec_diff == 0 ? "" : $tim_miss_length > 0
                    ? $tim_miss_length / $sec_diff
                    : ($ele ne "" && $ELE[$tim_miss_start] ne "") ? abs($ele - $ELE[$tim_miss_start]) / $sec_diff : 0;
                if ($sec_gradient != 0) {
                    debug_log "interpolating missing time over ".dis_string($tim_miss_length).", sec_diff = $sec_diff";
                    for (my $j = $tim_miss_start; $j < $i_end; $j++) {
                        $sec_curr += ($tim_miss_length > 0 ? $DIS[$j] : 1) * $sec_gradient;
                        if (!$seg && $sec_curr - $sec_prev < $min_input_timediff) {
                            ++$num_ignored_diff; # silently ignore
                            cut_trkpt($j);
                            $DIS[$j + 1] += $DIS[$j];
                            next;
                        }
                        $SEC[$j]  = int($sec_curr);
                        $TIM[$j]  = epoch_to_str(int($sec_curr));
                        $FLG[$j] |= INTERPOLATED_TIM | COMPUTED_SPEED;
                        $SPD[$j] = $const_spd;
                        debug_log "interpolating missing time at point $j with ele=$ELE[$j] new time=$TIM[$j] dis=$DIS[$j]";
                        $sec_prev = $sec_curr;
                    }
                    recalculate_SPD($i_end) if $i_end < $i;
                }
                $tim_miss_start = $i_start;
            }
            push @PREV, $#TIM;
            push @NEXT, $#TIM + 2;
            push @TXT, undef;
            push @IGN, 0;
            push @SEG, $seg;
            push @FLG, $flags;
            push @LAT, $lat;
            push @LON, $lon;
            push @ELE, $ele;
            push @TIM, $tim;
            push @SEC, $sec;
            push @DIS2,$sq_dis2d;
            push @DIS, $dis;
            push @SPD, $spd; # initially as recorded, "" if not available
            push @CMT, $cmt;
            push @SRC, $src;
            push @FIX, $fix;
            push @EXT, $exts;

            # detect and report missing altitude
            my $ele_solution  = "cannot calculate ascent/descent";
            $ele_solution  = "trying to interpolate altitude" if $ele_miss_start != $i_start;
            if ($ele eq "") {
                $num_no_ele++;
                $ele_gap_length += $dis; # using $dis == 0 for first point
                warn_trkpt($i, "no altitude     at", "; $ele_solution") if !$SEG[$i] && $num_no_ele == 1;
            } elsif ($num_no_ele) {
                # unfortunately, this message may show after initial ones on next segment:
                warn_trkpt($i - 1, "no altitude  up to",
                           "for ".points($num_no_ele)." over a distance of ".dis_string($ele_gap_length));
                $num_no_ele     = 0;
                $ele_gap_length = 0;
            }

            # detect and report missing time
            my $time_solution = "cannot calculate speed and acceleration";
            $time_solution = "trying to interpolate time" if $tim_miss_start != $i_start;
            unless (defined $sec) {
                $num_no_time++;
                $time_gap_length += $dis; # using $dis == 0 for first point
                warn_trkpt($i, "no time         at", "; $time_solution") if !$SEG[$i] && $num_no_time == 1;
            } elsif ($num_no_time) {
                # unfortunately, this message may show after initial ones on next segment:
                warn_trkpt($i - 1, "no time      up to",
                           "for ".points($num_no_time)." over a distance of ".dis_string($time_gap_length));
                $num_no_time     = 0;
                $time_gap_length = 0;
            }

            my $sec_str = defined $sec ? epoch_to_str($sec) : "";
            abort_trkpt($i, "internally inconsistent time data: tim=$tim != $sec_str (=$sec) at")
                if $tim =~ s/(\.\d*?)0+Z/$1Z/r ne $sec_str =~ s/(\.\d*?)0+Z/$1Z/r &&
                !($tim =~ m/\.999Z/ && defined $sec && $tim eq epoch_to_str($sec - 0.001));
            $i++;

            # will go wrong if previous trkpt gets ignored: $ELE[$#$ELE] = $ELE[$#$ELE - 1] if $ele eq "" && !$SEG[$#$ELE]; # implies $#$ELE > 0
            ($prev_lat, $prev_lon, $prev_ele, $prev_tim, $prev_sec, $prev_dis) =
            (     $lat,      $lon,      $ele,      $tim,      $sec,      $dis);
            $prev_avail_ele = $ele if $ele ne "";

            # enables checking negative time difference also across missing time entries:
            ($prev_avail_tim, $prev_avail_sec) = ($tim, $sec) if defined $sec;

            $state = 3; # ready to process further trkpt(s) in trkseg (or new trkseg or new track)
        }
      continuation:
        if (m#</trkseg>.*<trkseg>#s) {
            info_trkpt($i - 1, "end of  segment at") if $i != $i_start;
            $state = 2; # starting new trkseg
            # reset intra-segment values
        }
        if ($i != $i_start && m#<trk>#) { # further track
            info("further track, again from file $ARGV") unless m#<gpx[ >]#s;
        }
    }
    } # end of local
    close M if $merge;

    my $Len = dis_string($ele_gap_length); # TODO restrict scope
    warn_trkpt($#TIM, "no altitude  up to",
               "for " . $num_no_ele . " points over a distance of $Len") if $num_no_ele > 1;
    $Len = dis_string($time_gap_length);
    warn_trkpt($#TIM, "no time      up to",
               "for " . $num_no_time . " points over a distance of $Len") if $num_no_time > 1;
    info_trkpt($#TIM, "last   segment end") if $#TIM >= $i_start;
    warn_trkpt($#TIM, "                at",
               "was unable to interpolate missing altitude by end of last segment") if $ele_miss_start != $i_start;
    warn_trkpt($#TIM, "                at",
               "was unable to interpolate missing time by end of last segment") if $tim_miss_start != $i_start;
}

sub check_PREV_NEXT {
    my $i = $_[0];
    my $prev_i = PREV($i);
    my $npi = NEXT($prev_i);
    if (defined $i && defined $prev_i) {
        stacktrace() unless defined $npi && $i == $npi;
        abort_trkpt($i, "internal inconsistency: NEXT(PREV($i)) not defined")   if !defined $npi;
        abort_trkpt($i, "internal inconsistency: NEXT(PREV($i)) = ".$npi." != $i") if $i != $npi;
    }
    my $next_i = NEXT($i);
    my $pni = PREV($next_i);
    if (defined $i && defined $next_i) {
        stacktrace() unless defined $pni &&  $i == $pni;
        abort_trkpt($i, "internal inconsistency: PREV(NEXT($i)) not defined")   if !defined $pni;
        abort_trkpt($i, "internal inconsistency: PREV(NEXT($i)) = ".$pni." != $i") if $i != $pni;
    }
}

sub cut_trkpt { return remove_trkpt(shift, -1); }
sub remove_trkpt {
    my ($i, $count_ignored) = @_;
    $count_ignored //= 0;
    my ($prev_i, $next_i) = (PREV($i), NEXT($i));
    $PREV[$next_i] = $PREV[$i] if defined $next_i;
    $NEXT[$prev_i] = $NEXT[$i] if defined $prev_i;
    $PART_START_INDEX[$all] = $next_i unless defined $prev_i;
    $PART_END_INDEX  [$all] = $prev_i unless defined $next_i;
    # check_PREV_NEXT($prev_i);
    # check_PREV_NEXT($next_i);
    return if $count_ignored == -1;

    # optionally indicate that previous point has been ignored, in addition to any points ignored before the removed point:
    # debug "IGN[".($next_i // 0)."]==".($IGN[$next_i // 0])." += IGN[".$i."]==$IGN[$i] + $count_ignored";
    $IGN[$next_i // 0] += $IGN[$i] + $count_ignored; # do this also past end (wrapping around in array)
    return unless defined $next_i;

    $TXT[$next_i] = defined $TXT[$i] ? $TXT[$i] . ($TXT[$next_i] // "") : $TXT[$next_i]; # not actually needed before points get inserted on time gaps
    $SEG[$next_i] = 1 if !defined $prev_i || $SEG[$i]; # copy over any segment start
    recalculate_all($next_i, $i);
}

sub merge_tracks
{
    my $in_merge = 0; # 1 means using trkpt from $merge, < 0: EOF in $merge, -1: should report its EOF
    my $infile_EOF = 0; # how often infile EOF has been hit
    my $skip_point = 0;
    my $prev_input_sec;
    our $i = 0; # current regular point
    our $j = $PART_END_INDEX[$all] + 1; # current alternative point
    abort("FATAL: no trkpt found in alternative input $merge") if $j > $#TIM;

    sub warn_skipped_merged_ratio
    {
        my $i = $_[0];
        my $merged  = $num_merged_points - $num_merged_points_last;
        my $skipped = $num_skip_points   - $num_skip_points_last;
        warn_trkpt($i, "skipped $skipped regular points while merging $merged ones before")
            if $skipped > $merged;
    }

    sub insert_trkpt { # insert $j virtually before $i
        # my ($i, $j) = @_;
        my $before_i = defined $i ? PREV($i) : $PART_END_INDEX[$all];
        $PREV[$j]= $before_i;
        $NEXT[$j]= $i;
        $TXT[$j] = undef;
        $SEG[$j] = defined $i ? $SEG[$i] : 0; # ignore segments of alternative input
        $SEG[$i] = 0 if defined $i; # has been moved over to inserted point just before
        $NEXT[$before_i] = $j if defined $before_i;
        if (defined $i) {
            $PREV[$i] = $j;
            if ($PART_START_INDEX[$all] == $i) {
                $PART_START_INDEX[$all]  = $j;
            }
        } else {
            $PART_END_INDEX[$all] = $j;
        }
        # check_PREV_NEXT($i);
        # check_PREV_NEXT($j);
        recalculate_all($j, $i);
    }

    while (1) {
        my ($lat , $lon , $ele ) = ($LAT[$i], $LON[$i], $ELE[$i]) if defined $i;      # of current point from regular input
        my ($lat2, $lon2, $ele2) = ($LAT[$j], $LON[$j], $ELE[$j]) if $in_merge != -2; # of current point from alternative input
        my ($tim , $sec ) = ($TIM[$i], $SEC[$i])                  if defined $i;
        my ($tim2, $sec2) = ($TIM[$j], $SEC[$j])                  if $in_merge != -2;
        debug_log "at ".(defined $i ? point_str($i) : "infile EOF")." vs. ".($in_merge < 0 ? "alternative EOF" : point_str($j));

        our ($prev_i, $next_i) = (PREV($i), NEXT($i));
        if (!$skip_point && defined $prev_i && defined $SEC[$prev_i] && defined $sec) {
            my $diff_time = $sec - $SEC[$prev_i]; # difference from previous, possibly merged point
            abort_trkpt($i, "before", "negative time difference (-".timediff_string(-$diff_time).")")
                if $diff_time < 0;
            if ($tim eq $TIM[$prev_i]) { # really equal (before rounding seconds)
                abort_trkpt($i, "point before",
                            "has same time but different position: $LAT[$prev_i], $LON[$prev_i] at")
                    if $lat != $LAT[$prev_i] || $lon != $LON[$prev_i];
                warn_trkpt($i, "point before      ",
                           "has same time but different altitude: $ELE[$prev_i] m")
                    if $ele != $ELE[$prev_i];
            }
        }

        # main strategy : use alternative input when it has more points during same period
        sub num_points_more { # number of extra alternative points before next regular point if exists,
            # otherwise all remaining points or 0 if time information (seconds) not available for current points
            my $k = $#TIM + 1;
            if (defined $next_i) { # not already at last regular input point
                $k = $j; # the current alternative point
                my $next_sec = $SEC[$next_i]; # time of next regular point
                if (defined $next_sec) {
                    $next_sec -= $MIN_TIMEDIFF; # do not count points too close to next regular point
                    while ($k <= $#TIM && defined $SEC[$k] && $SEC[$k] < $next_sec) {
                        $k++;
                    }
                }
            }
            return $k - $j;
        }
        my $num_points_more;
        my $next_sec = defined $next_i && defined $SEC[$next_i] ? $SEC[$next_i] : $sec;

        my $start_merge = "start using alternative input";
        if ($in_merge == 0 || $in_merge == -2) {
            # skipping regular input, or currently not merging from alternative input, which not just reached its EOF
            my $timediff_merge = diff_defined($sec2, $sec) // 0;
            debug_log "timediff_merge $timediff_merge";

            # skipping continues until point from regular input is past next alternative input
            if ($in_merge == 0 && $skip_point && !($timediff_merge < 0)) {
                warn_trkpt($i, "due to time diff ".timediff_string($timediff_merge).
                           ", skipping further regular");
                remove_trkpt($i);
                $num_skip_points++;
                goto continuation;
            }

            my $gap_time = diff_defined($sec, $prev_input_sec) // 0; # difference from previous point from infile(s);
            if ($skip_point ||
                ($timediff_merge <= -$MIN_TIMEDIFF # regular input went sufficiently past current alternative input
                 && (!defined $prev_i # alternative input started before regular one
                     || ($gap_time >= MERGE_MIN_GAP &&
                         ($num_points_more = num_points_more()) >= MERGE_MIN_MORE)))) {
                # (re-)start merge
                $in_merge = 1;
                $num_merges++ unless $skip_point;
                $num_merged_points++;
                unless (defined $prev_i) {
                    info("$start_merge before regular input");
                    debug_log point_str($j).": switched to merge before regular input";
                    insert_trkpt();
                    $FLG[$j] |= START_MERGE;
                } elsif ($skip_point) {
                    debug_log point_str($j).": continue merge";
                    insert_trkpt();
                    $skip_point = 0;
                } else { # certainly defined $prev_i
                    # difference of next point from previous point (both from regular input):
                    my $gap_time2 = timediff_string(diff_defined($next_sec, $prev_input_sec) // 0);
                    info_trkpt($i, "$start_merge due to $num_points_more more points during $gap_time2 around");
                    my ($prev_lat , $prev_lon , $prev_ele ) = ($LAT[$prev_i], $LON[$prev_i], $ELE[$prev_i]);
                    my $dis1_from_merge = distance($lat2, $lon2, $ele2, $prev_lat, $prev_lon, $prev_ele);
                    my $dis2_from_merge = distance($lat2, $lon2, $ele2,      $lat,      $lon,      $ele);
                    my $dis_from_merge = ($dis1_from_merge + $dis2_from_merge) / 2;
                    insert_trkpt();
                    $FLG[$j] |= START_MERGE;
                    warn_trkpt($j, "points on regular track are ".dis_string($dis_from_merge)." away from first point being merged:  ")

                        if $dis_from_merge > MERGE_WARN_DIST;
                    if ($skip_point) {
                        debug_log point_str($i).": skipping regular point while switching to merging $tim2";
                        remove_trkpt($i);
                        $num_skip_points++;
                        $i = $next_i;
                        $infile_EOF++ unless defined $i;
                        goto continuation;
                    }
                }
            } else {
                debug_log point_str($i).": using regular input";
                if ($in_merge == 0) {
                    # advance alternative input until next point no more too close to regular point
                    while ($in_merge != -2 &&
                           (((!defined $sec || !defined $sec2) &&
                             $lat == $lat2 && $lon == $lon2 && eq_nonempty($ele, $ele2))
                            || $sec2 < $sec + $MIN_TIMEDIFF)) {
                        # advance $merge as far as possible to get past regular input
                        debug_log point_str($j).": skipping merge";
                        $in_merge = -2 if ++$j > $#TIM;
                        ($tim2, $sec2) = ($TIM[$j], $SEC[$j]);
                    }
                }
            }
        } else { # currently merging from alternative input $merge or just reached its EOF
            if ($in_merge == -1 # EOF in merge, implies defined $i (not at infile EOF)
                || (defined $i && defined $sec && defined $sec2 &&
                    $sec2 + $MIN_TIMEDIFF > $sec) # alternative input from $merge is no more at sufficiently earlier time than regular input
               ) {
                # return to next point from regular input unless current regular point to be skipped
                if ($in_merge == 1) {
                    $in_merge = 0;

                    my ($prev_lat2, $prev_lon2, $prev_ele2) = ($LAT[$j - 1], $LON[$j - 1], $ELE[$j - 1]);
                    my (     $lat2,      $lon2,      $ele2) = ($LAT[$j    ], $LON[$j    ], $ELE[$j    ]);
                    my $dis1_from_infile = distance($lat, $lon, $ele, $prev_lat2, $prev_lon2, $prev_ele2);
                    my $dis2_from_infile = distance($lat, $lon, $ele,      $lat2,      $lon2,      $ele2);
                    my $Dis = dis_string(($dis1_from_infile + $dis2_from_infile) / 2);
                    $num_points_more = num_points_more();
                    if ($num_points_more > 0) { # alternative input still has more points
                        info_trkpt($i, "alternative input is $Dis away and has $num_points_more more points, skipping regular");
                        debug_log point_str($i).": skipping regular input";
                        $skip_point = 1;
                        remove_trkpt($i);
                        $num_skip_points++;
                        goto continuation; # skip current point from regular input
                    }
                    warn_trkpt($prev_i, "next point on regular track is " .
                               dis_string($dis1_from_infile) . " away from the most recently merged")
                        if $dis1_from_infile > MERGE_WARN_DIST;
                }
                my $end_diff = timediff_string(diff_defined($sec, $SEC[$prev_i]) // 0);
                info_trkpt($prev_i, ($in_merge == 0 ? "pause using" : "     end of").
                           " alternative input,          next regular point $end_diff after");
                $in_merge = -2 if $in_merge == -1;
                recalculate_all($i, $prev_i);
                $FLG[$i] |= PAUSE_MERGE;
                debug_log point_str($i).": switched back to base input";
                warn_skipped_merged_ratio($i);
                ($num_merged_points_last, $num_skip_points_last) =
                ($num_merged_points     , $num_skip_points     );
            } else {
                insert_trkpt();
                $num_merged_points++;
                if ($infile_EOF == 1) {
                    $infile_EOF = 2;
                    $num_merges++;
                    my $inserted = $PART_END_INDEX[$all];
                    info_trkpt($inserted, "$start_merge after end of regular input;   first point at", "") unless $skip_point;
                    $FLG[$inserted] |= START_MERGE unless $skip_point;
                    debug_log point_str($j).": switched to merge after regular input";
                } else {
                    debug_log point_str($j).": continued from merge";
                }
            }
        }

        $prev_input_sec = $SEC[$i] if defined $i && defined $SEC[$i];
      continuation:
        if ($in_merge == 1) {
            if (++$j > $#TIM) { # EOF in $merge
                $in_merge = -1;
                last unless defined $i;
            }
        } else {
            unless (defined (($i = $next_i))) {
                $infile_EOF++;
                last if $in_merge < 0; # EOF in $merge
                $in_merge = 1;
            }
      }
    }
    warn_skipped_merged_ratio($PART_END_INDEX[$all]);
    info_trkpt($PART_END_INDEX[$all], "     end of alternative input                                        after")
        if $in_merge == -1;
}

sub read_all
{
    read_tracks(); # regular input file(s)
    $PART_START_INDEX[$all] = 0;
    $PART_END_INDEX  [$all] = $#TIM;
    $PREV[$PART_START_INDEX[$all]] = undef;
    $NEXT[$PART_END_INDEX[  $all]] = undef;
    if ($merge) {
        read_tracks($merge); # alternative input
        merge_tracks();
    }
}
read_all();

## sanitize trackpoints with implausible data

my $avg_timediff if $smoothing;

my ($num_removals_ele, $num_removals_dir, $num_removals_spd, $num_removals_acc, $num_removals_dec,
    $num_removals_gain, $num_removals_loss, $num_removals_diff) = (0, 0, 0, 0, 0, 0, 0, 0);
sub sanitize
{
    my $num_removals_in_row = 0;
    my $sum_timediff_mov0 = 0;
    my $num_pts_mov       = 0;
    for (my $i = $PART_START_INDEX[$all]; defined $i; $i = NEXT($i)) {
        # check_PREV_NEXT($i);
        sub ignore_trkpt {
            my ($i, $Spd, $reason) = @_;
            my $Spd_txt = $Spd eq "" ? "" : "at speed = $Spd km/h ";
            warn_trkpt(shift @_, "ignoring", $Spd_txt."as $reason");
        }

        my ( $lat    , $lon    , $ele    , $tim    , $sec    , $spd    ) =
            ($LAT[$i], $LON[$i], $ELE[$i], $TIM[$i], $SEC[$i], $SPD[$i]);
        my ($prev_i, $next_i) = (PREV($i), NEXT($i));

        debug_log "consider sanitizing ".point_str($i)." spd=".spd_or_none($spd // "");
        my $ignore = 0;

        # handle cases that do not involve earlier points and lead to removing current point
        if ($ele ne "" && $ele < MIN_PLAUSIBLE_ELE) {
            $num_removals_ele++;
            $ignore = 1;
            ignore_trkpt($i, "", "altitude is less than ".MIN_PLAUSIBLE_ELE);
        } elsif ($spd ne "" && ($spd < 0 || MPS2KMPH * $spd > $MAX_PLAUSIBLE_SPEED)) {
            $num_removals_spd++;
            $ignore = 1;
            ignore_trkpt($i, spd_string($spd), "recorded speed is ". ($spd < 0 ? "negative" : "higher than ".$MAX_PLAUSIBLE_SPEED));
        }

        if (!$ignore && !$SEG[$i]) { # not already ignoring current point, not at segment start

            # correct implausible elevation gain/loss after somewhat long gap, typically on exit of tunnel or building
            if (!section($i) && # not at merge start/pause
                !$SEG[$prev_i] && !section($prev_i)) {
                my $prev_timediff = diff_defined($SEC[$prev_i], $SEC[PREV($prev_i)]);
                my $rate = rate(diff_nonempty($ele, $ELE[$prev_i]), diff_defined($sec, $SEC[$prev_i]));
                if (defined $prev_timediff && $prev_timediff >= MIN_TUNNEL_GAP &&
                    $rate ne "" && ($rate > $MAX_PLAUSIBLE_ELE_GAIN / 4 || $rate < -$MAX_PLAUSIBLE_ELE_LOSS / 4)) {
                    warn_trkpt($prev_i, "      at",
                               "correcting elevation from ".ele_string(0, $ELE[$prev_i])." to ".ele_string(0, $ELE[$i])
                               ." m after gap of ".timediff_string($prev_timediff)." and ".dis_string($DIS[$prev_i]));
                    $ELE[$prev_i]  = $ELE[$i];
                    $FLG[$prev_i] |= CORRECTED_ELE;
                    recalculate_all($prev_i, $prev_i); # 2d part did not change
                    recalculate_all($i, $i); # 2d part did not change
                }
            }

            # handle cases that involve point(s) before and may lead to removing current point
            my ($diff_ele, $timediff, $rate, $dis, $comp_spd, $acc) = comp_diffs($i);
            my ($theta_diff, $spd2) = comp_angle($i, $timediff, $dis); # for spike detection
            $ANG[$i] = $theta_diff;

            if (!section($i) && # not at merge start/pause
                !($FLG[$i] & COMPUTED_SPEED)) { # not already computed speed (due to interpolation of missing time)
                if ($comp_spd ne "" && $spd ne "" && $timediff ne "" && $timediff <= MAX_TIMEDIFF_RECORDED_SPEED) {
                    my $spd_diff = abs($comp_spd - $spd);
                    $spd_prec_len = max($spd_prec_len, length(spd_str($spd_diff))); # used only for spd_string
                    $sum_speed_deviation += $spd_diff;
                    $sum_speed += $comp_spd + $spd; # will later divide by 2 for average;
                    warn_trkpt($i, "before  ",
                               "speed recorded: ".spd_string($spd)." computed: ".spd_string($comp_spd)." ".
                               "deviation: "    . spd_string($spd_diff) . " km/h" .
                               " (more than " . $MAX_SPEED_DEVIATION . ")") if $spd_diff > ($MAX_SPEED_DEVIATION + 0.5) / MPS2KMPH;
                }
                # from now on, use for $spd and $SPD[$i] the computed speed (average since previous point) if recorded speed is not available or previous point is not very recent
                if (($spd eq "" || $timediff ne "" && $timediff > MAX_TIMEDIFF_RECORDED_SPEED)) {
                    $SPD[$i]  = $spd = $comp_spd;
                    $FLG[$i] |= COMPUTED_SPEED if $comp_spd ne "";
                }
            }

            my $Spd      = $spd ne "" ? spd_string($spd) : "";
            my $Spd_txt  = $spd ne "" ? "at speed = $Spd km/h " : "";
            my $Dis      = dis_string($dis);
            my $Rate     = $rate ne "" ? sprintf RAT_PRECISION, $rate : "";
            my $Acc      = $acc ne "" ? sprintf ACC_PRECISION, $acc : "";
            my $Dec      = $acc ne "" ? sprintf ACC_PRECISION, -$acc : "";
            warn_trkpt($i, "before  ", $Spd_txt."distance between points = $Dis (more than ".(WARNING_TPT_DIST / 1000).")") if $dis > WARNING_TPT_DIST;

            $ignore = 1; # tentatively indicate removal of current point
            if (MAX_PLAUSIBLE_ANGLE_DIFF < abs($theta_diff // 0)) { # implies !section($i) etc.
                $num_removals_dir++;
                $spd_prec_len = max($spd_prec_len, length(spd_str($spd2))); # used only for spd_string
                my $Theta_diff = sprintf "%+.0f", $theta_diff;
                ignore_trkpt($i, spd_string($spd2), "turned by $Theta_diff degrees (".($theta_diff >= 0 ? "more than  " : "less than -").MAX_PLAUSIBLE_ANGLE_DIFF.")");
            } elsif (!section($i) && $spd ne "" && $Spd > $MAX_PLAUSIBLE_SPEED) {
                $num_removals_spd++;
                ignore_trkpt($i, $Spd, "speed is higher than ".$MAX_PLAUSIBLE_SPEED);
            } elsif (!section($i) && $acc ne "" &&  $acc > $MAX_PLAUSIBLE_ACCEL &&
                     ($FLG[$prev_i] & INTERPOLATED_TIM) == ($FLG[$i] & INTERPOLATED_TIM)) {
                $num_removals_acc++;
                ignore_trkpt($i, $Spd, "acceleration = $Acc m/s/s (more than $MAX_PLAUSIBLE_ACCEL)");
            } elsif (!section($i) && $acc ne "" && -$acc > $MAX_PLAUSIBLE_DECEL &&
                     ($FLG[$prev_i] & INTERPOLATED_TIM) == ($FLG[$i] & INTERPOLATED_TIM)) {
                $num_removals_dec++;
                ignore_trkpt($i, $Spd, "deceleration = $Dec m/s/s (more than $MAX_PLAUSIBLE_DECEL)");
            } elsif (!section($i) && $rate ne "" && $rate > $MAX_PLAUSIBLE_ELE_GAIN) {
                $num_removals_gain++;
                ignore_trkpt($i, $Spd, " ascent rate = $Rate m/h (more than  $MAX_PLAUSIBLE_ELE_GAIN)");
            } elsif (!section($i) && $rate ne "" && $rate < -$MAX_PLAUSIBLE_ELE_LOSS) {
                $num_removals_loss++;
                ignore_trkpt($i, $Spd, "descent rate = $Rate m/h (less than -$MAX_PLAUSIBLE_ELE_LOSS)");
            } elsif ($timediff ne "" && $timediff < $MIN_TIMEDIFF &&
                     ($FLG[$prev_i] & INTERPOLATED_TIM) == ($FLG[$i] & INTERPOLATED_TIM))
            {
                $num_removals_diff++;
                my $inter = $FLG[$i] & INTERPOLATED_TIM ? "interpolated " : "";
                ignore_trkpt($i, $Spd, $inter."time difference = $timediff sec (less than $MIN_TIMEDIFF)"); # " over a distance of $Dis"
            } else {
                $ignore = 0;
                if ($timediff ne "" && $spd ne "" && $spd >= MIN_SPEED_MOVING / MPS2KMPH)
                {
                    $sum_timediff_mov0 += $timediff; # used for smoothing
                    $num_pts_mov++;                  # used for smoothing
                }
            }
        }

        if ($ignore) {
            $num_removals_in_row++;
            remove_trkpt($i, 1); # also recalculates speed etc. of former point $next_i relative to its new previous point
        }

        # debug "$i $tim $num_removals_in_row vs. $IGN[$i] $IGN[$next_i]";
        my $num_removals_in_row2 = $IGN[$i];
        my $curr_not_ignored = ($IGN[$next_i // 0] // 0) == 0; # also past end (wrapping around in array)
        # warning("ignored $num_removals_in_row2 points in a row") if $num_removals_in_row2 > 1 && $curr_not_ignored;
        if (!$ignore) {
            abort_trkpt($i, "$tim internal inconsistency on number of most recent removals in a row: ".
                        "$num_removals_in_row vs. $num_removals_in_row2, condition 1 vs. $curr_not_ignored at")
                if $num_removals_in_row != $num_removals_in_row2 || !$curr_not_ignored;
            warning("ignored $num_removals_in_row points in a row") if $num_removals_in_row > 1;
            $num_removals_in_row = 0;
        } elsif ($curr_not_ignored) {
            abort_trkpt($i, "internal inconsistency on condition for reporting number of most recent removals in a row at");
        }
    }
    my $num_removals_in_row_final = $IGN[$PART_START_INDEX[$all]]; # past end, wrapping around in array
    # warning("ignored $num_removals_in_row_final points in a row") if $num_removals_in_row_final > 1;
    abort_trkpt($PART_END_INDEX[$all], "internal inconsistency on reporting number of final removals in a row: ".
                "$num_removals_in_row vs. $num_removals_in_row_final at")
        if       $num_removals_in_row !=  $num_removals_in_row_final;
    warning("ignored $num_removals_in_row points in a row") if $num_removals_in_row > 1;
    # note that in the GPX output information on points sanitzed at the end of the track will get lost

    $avg_timediff = $num_pts_mov ? $sum_timediff_mov0 / $num_pts_mov : 0 if $smoothing;
    # $SUM_TIMEDIFF_MOV[$all] will be calculated independently, as points may have been ignored or smoothened
}
sanitize() if $sanitize;

abort("FATAL: no trackpoint (possibly after ignoring some)") if $#TIM < 0;


sub check_waypoints {
    for (my $j = 0; $j <= $#WSEC; $j++) {
        my $id = $WSTR[$j];
        next unless defined $id; # just info/warning wpt
        my $sec = $WSEC[$j];
        unless (defined $sec && $sec) {
            warning("waypoint $id does not have time information; skipping proximity check")
                unless $id =~ m/\[(INFO|WARNING): /;
            next;
        }
        if ($j > 0 && !($WSTR[$j] =~ m|<src>GPXConv</src>|) && (diff_defined($sec, $WSEC[$j - 1]) // 0) < 0) {
            warning("waypoint $id has earlier time than the waypoint found before: $WSTR[$j - 1]");
        }
        if ($#TIM < 0) { # can not happen due to above check
            warning("no trackpoint found for checking waypoints");
            last;
        }
        # find trackpoint $i2 at or just after wpt
        my $i2 = 0; # would be more efficient to do this just before the wpt loop, but then cannot deal with wpts not in time order
        my $i1; # most recent trackpoint having time info, otherwise undefined; ulimately becomes same as $i2 if no such trackpoint exists
        while (defined NEXT($i2) && (!defined $SEC[$i2] || $SEC[$i2] < $sec))
        {
            $i1 = $i2 if defined $SEC[$i2];
            $i2 = NEXT($i2);
        }
        $i1 = $i2 unless defined $i1;

        my ($wlat, $wlon, $wele) = ($WLAT[$j], $WLON[$j], $WELE[$j]);
        next if $LAT[$i1] == $wlat && $LON[$i1] == $wlon && $ELE[$i1] eq $wele; # late finishing point
        next if $LAT[$i2] == $wlat && $LON[$i2] == $wlon && $ELE[$i2] eq $wele; # early starting point
        # due to time possibly carried over across segement end and start of next segment, wpt may match next trkpt:
        my $next_i2 = NEXT($i2);
        next if defined $next_i2 && $SEC[$i2] == $sec && defined $SEC[$next_i2] && $SEC[$next_i2] == $sec
            && $LAT[$next_i2] == $wlat && $LON[$next_i2] == $wlon && $ELE[$next_i2] eq $wele;
        warning("no trackpoint found at same time or later than waypoint $id; will ".
                ($i1 < $i2 ? "extrapolate distance from the one before" : "take the difference from the last one"))
            unless defined $next_i2;

        # calculate and report distance to nearest trackpoint
        my $diff_sec_i2_i1 = diff_defined($SEC[$i2], $SEC[$i1]);
        my $linear_factor = $diff_sec_i2_i1 ? ($sec - $SEC[$i1]) / $diff_sec_i2_i1 : 0;
        my $lat = $LAT[$i1] + ($LAT[$i2] - $LAT[$i1]) * $linear_factor;
        my $lon = $LON[$i1] + ($LON[$i2] - $LON[$i1]) * $linear_factor;
        my $ele2 = $ELE[$i2] eq "" ? ($WELE[$j] eq "" ? $ELE[$i1] : $WELE[$j]) : $ELE[$i2];
        my $ele1 = $ELE[$i1] eq "" ? $ele2 : $ELE[$i1];
        my $ele = $ele1 eq "" || $ele2 eq "" ? ""
            : $ele1 + ($ele2 - $ele1) * $linear_factor;
        my $dis = distance($lat, $lon, $ele, $wlat, $wlon, $wele);
        my $tim = "NO"; # do not print time for interpolated position, since done already for wpt
        warn_trkpt(-1, "according to its time information, waypoint $id appears not close to track "
                   . "(interpolated distance = ".(sprintf "%.0f", $dis)." m > " . WARNING_WPT_DIST . "); "
                   . "it should be approximately at", "", $lat, $lon, $ele, $tim)
          if $dis > WARNING_WPT_DIST;
    }
}
check_waypoints();

## smoothen trackpoints

# smoothen each point, by weightened average with (smoothened) previous point and (original) next point
sub average {

    sub weight {
        #   my $d = 1.5*$avg_timediff; # max dist of influence
        #   return max(0, ($d-$_[0])/$d);
        return 1 / (1 + $_[0] / $avg_timediff);
    }
    my $x1 = $_[0];
    my $x  = $_[1];
    my $x2 = $_[2];
    my $diff_time1 = $_[3];
    my $diff_time2 = $_[4];

    #return $x;
    return ($x1 * weight($diff_time1) + $x + $x2 * weight($diff_time2)) /
      (weight($diff_time1) + 1 + weight($diff_time2));
}

if ($smoothing) {
    my $pending_recalculation = 0;
    my ($prev_i, $i, $next_i);
    for ($i = $PART_START_INDEX[$all]; defined $i && defined $next_i; $prev_i = $i, $i = $next_i) {
        # check_PREV_NEXT($i);
        $next_i = NEXT($i);
        next if !defined $prev_i;

        my $sec = $SEC[$i];
        my $diff_time1 = diff_defined($sec, $SEC[$prev_i]);
        my $diff_time2 = diff_defined($SEC[$next_i], $sec);
        if (defined $diff_time1 && defined $diff_time2 && !$SEG[$i] && !$SEG[$next_i] # && $section($i) && !section($next_i)
            && max($diff_time1, $diff_time2) <= SMOOTHING_MAX_GAP) {
            my ($prev_ele, $ele, $next_ele) = ($ELE[$prev_i], $ELE[$i], $ELE[$next_i]);
            my ($old_lat, $old_lon) = ($LAT[$i], $LON[$i]);
            $LAT[$i] = average($LAT[$prev_i], $LAT[$i], $LAT[$next_i], $diff_time1, $diff_time2);
            $LON[$i] = average($LON[$prev_i], $LON[$i], $LON[$next_i], $diff_time1, $diff_time2);
            $ELE[$i] = average($prev_ele, $ele, $next_ele , $diff_time1, $diff_time2)
                if $prev_ele ne "" && $ele ne "" && $next_ele ne "";
            recalculate_all($i, $old_lat == $LAT[$i] && $old_lon == $LON[$i] ? $i : undef);
            $pending_recalculation = 1;
            # individiual points are not marked as smoothened; this is done globally
            $num_smoothened++;
        } elsif ($pending_recalculation) {
            recalculate_all($i);
            $pending_recalculation = 0;
        }
    }
    recalculate_all($i) if $pending_recalculation;
}

## perform general statistics, optionally also for parts (segements or days)
## with optional ascent/descent phase analysis and optional energy calculation

# TODO support partial energy computation also on missing elevatin or speed
my $num_energy_missing = 0;
my $energy_missing_ele = 0;
my $energy_missing_spd = 0;
my $cannot_energy = "; cannot calculate energy here";
sub energy { # calculate energy difference between trackpoints  # TODO simplify provisioning of parameters using index
    my $i  = $_[0];                # index of end of the current span between points with assumed linear movement
    my $t0 = $_[1];                # start time in s (relative to begin of track)
    my $t1 = $_[2];                # end   time in s (relative to begin of track)
    my $dt = $t1 - $t0;            # duration in s
    my $d  = $_[3];                # 3-dimensional distance in m
    my $h0 = $_[4];                # altitude before in m
    my $h1 = $_[5];                # altitude after  in m
    if ($h0 eq "" || $h1 eq "") {
        # $h0 = $h1 = 500;         # default altitude in m assumed for calculating density
        warn_trkpt($h1 eq "" ? $i : PREV($i), "no altitude    at", $cannot_energy) if $energy_missing_ele < 2;
        $energy_missing_ele++;
        $num_energy_missing++;
        return 0;
    }
    my $dh =  $h1 - $h0;           # altitude difference (or 0 if default elevation used)
    my $h  = ($h0 + $h1) / 2;      # average altitude, used for air density
    my $v0 = $_[6];                # speed before in m/s
    my $v1 = $_[7];                # speed after  in m/s
    if ($v0 eq "" || $v1 eq "") {  # $v1 eq "" cannot happen because times are present, lat/lon are guaranteed to be there, just ele may be missing
        # $v0 = $v1 = 0;           # default speed in m/s assumed for calculating acceleration and air drag
        warn_trkpt($v1 eq "" ? $i : PREV($i), "no speed info  at", $cannot_energy) if $energy_missing_spd < 2;
        $energy_missing_spd++;
        $num_energy_missing++;
        return 0;
    }
    my $v  = ($v0 + $v1) / 2;      # average speed, used for air drag, assumed with no wind (i.e., calm air)
    my $g = 9.81;                  # gravity constant in m/(s^2)

# values and formulas adapted from https://www.sheldonbrown.com/rinard/aero/formulas.html
    my $m = $weight;               # total mass including vehicle and rider in kg
    my $cr = 0.003;                # coefficient of rolling friction for racing bike      - adapt
    my $cdA = 0.39;                # coefficient of wind resistance * frontal area in m^2 - adapt
    my $r = 1.2;                   # air density on average weather (20 degree Celsius) at sea level in kg/(m^3)
    my $density = $r - (0.55 * $r) * $h / 5000; # roughly corrected air density depending on altitude

    # energy components and their sum along the distance $d in Joule
    my $roll  = $cr * $m * $g * $d;
    my $slope = $m * $g * $dh;
    my $accel = $m / 2 * ($v1 * $v1 - $v0 * $v0);
    my $drag  = $density * $cdA / 2 * $v * $v * $d;
    my $sum   = $roll + $slope + $accel + $drag;

    ## my $debug_inputs = "dt=" . int($dt) . " d=" . int($d) . " dh=" . int($dh) . " v0=" . int($v0) . " v1=" . int($v1);
    ## my $debug_outputs = "roll:" . int($roll) . " + slope:" . int($slope) . " + accel:" . int($accel) . " + drag:" . int($drag) . " = sum:" . int($sum);
    ## debug "energy @ time=$t  $debug_inputs   =>   $debug_outputs";

    return $sum;
}

sub push_phase {
    my $start_index    = $_[0];
    my $end_index      = $_[1];
    my $sec1           = $SEC[$start_index];
    my $sec2           = $SEC[$end_index];
    my $duration        = diff_defined($sec2, $sec1);
    my $diff_ele       = $ELE[$end_index] - $ELE[$start_index];
    my $dist           = 0;
    my $rate           = 0;
    my $spd_           = 0;
    my $max_rate       = 0;
    my $max_rate_index;
    my $max__spd       = 0;
    my $max__spd_index;

    my $i = NEXT($start_index);
    for (my $prev_i = PREV($i);
         $i <= $end_index;
         $prev_i = $i, $i = NEXT($i)) {
        # check_PREV_NEXT($i);
        $dist += $DIS[$i];
        next if section($i);

        $rate = rate(diff_nonempty($ELE[$i], $ELE[$prev_i]), diff_defined($SEC[$i], $SEC[$prev_i]));
        $spd_ = $SPD[$i];
        ($max_rate, $max_rate_index) = ($rate, $i) if $rate ne "" && ($diff_ele >= 0 ? $rate > $max_rate : $rate < $max_rate);
        ($max__spd, $max__spd_index) = ($spd_, $i) if $spd_ ne "" && $spd_ > $max__spd;
    }
    my $avg_spd  = $duration ? $dist / $duration : "";
    my $avg_rate = rate($diff_ele, $duration);
    push @PHASE_DURATION,       $duration;
    push @PHASE_DIFF,           $diff_ele;
    push @PHASE_DIST,           $dist;
    push @PHASE_SPD,            $avg_spd;
    push @PHASE_RATE,           $avg_rate;
    push @PHASE_MAX_SPD_INDEX,  $max__spd_index;
    push @PHASE_MAX_RATE,       $max_rate;
    push @PHASE_MAX_RATE_INDEX, $max_rate_index;
    push @PHASE_END_INDEX,      $end_index;
}

my $part;
sub next_part {
    my $day_or_seg = $_[0];
    $part = $day_or_seg if $day_or_seg >= $part_start;
    $part = undef       if defined $part_end && $day_or_seg > $part_end;
}
my $seg = 0 if $segs;
my $day = 0, my $prev_date, my $date if $days;
my ($part_start_sec, $part_end_sec);

# phase analysis ignores merge start/pause
for (my $direction, my $prev_sec, my $timediff,
     my $prev_turn_index, my $prev_turn_ele, my $prev_turn_sec,
     my $local_max_index = 0, my $local_min_index = 0,
     my $local_max_sec, my $local_min_sec,
     my $local_max_ele = "", my $local_min_ele = "",
     my $prev_i, my $next_i,
     my $i = $PART_START_INDEX[$all]; defined $i; $prev_i = $i, $i = $next_i) {
    # check_PREV_NEXT($i);
     $next_i = NEXT($i);

    my ( $lat    , $lon    , $ele    , $tim    , $sec    , $dis    , $spd    ) =
        ($LAT[$i], $LON[$i], $ELE[$i], $TIM[$i], $SEC[$i], $DIS[$i], $SPD[$i]);
    if (!$SEG[$i]) {
        $prev_sec = $SEC[$prev_i];
        $timediff = diff_defined($sec, $prev_sec) // "";
    }

    # local inititialization of optional statistics for parts (segements or days)
    my $new_part = 0;
    if ($days) {
        if (defined $tz && !defined $min_sec && defined $sec) {
            # https://stackoverflow.com/questions/2143528/whats-the-best-way-to-get-the-utc-offset-in-perl
            use Time::Local;
            my @t = localtime($sec);
            my $system_offset_from_UTC_in_seconds = timegm(@t) - timelocal(@t);
            $sec_offset = $tz * H2S - $system_offset_from_UTC_in_seconds;
            # more accurate but more involed alternative: https://metacpan.org/pod/Geo::Location::TimeZone
        }
        $date = defined $sec ? join("-", (localtime(defined $sec_offset ? $sec + $sec_offset : $sec))[3,4,5]) : undef;
        $new_part = 1
            if !defined $prev_i
            # assuming that dates, as far as available, are ascending:
            || (defined $date && defined $prev_date && $date ne $prev_date);
        $prev_date = $date if $days;
    }
    $new_part = 1 if $segs && $SEG[$i];
    $DURATION[$part] = diff_defined($part_end_sec, $part_start_sec)
        if defined $part && ($new_part || !defined NEXT($i)); # implies $i != 0
    if ($new_part) {
        next_part(++$day) if $days;
        next_part(++$seg) if $segs;
    }
    if (defined $part && $new_part) {
        $PART_START_INDEX[$part] = $i;
        $part_start_sec = $part_end_sec = $SEC[$i];
    } else {
        $part_start_sec = $part_end_sec = $SEC[$i] unless defined $part_start_sec;
    }
    # preliminary values are updated on each iteration:
    $part_end_sec = $SEC[$i] if defined $SEC[$i];
    $PART_END_INDEX[$part] = $i if defined $part;

    if ($SEG[$i]) { # track segment start (which includes the start: !defined $prev_i)
        # (re-)start climb phase recognition
        $direction       = 0;
        $prev_turn_ele   = $local_max_ele   = $local_min_ele = "";
    }
    if ($ele ne "") {
        ($prev_turn_ele, $prev_turn_index) = ($ele, $i) if $prev_turn_ele eq "";
        ($local_max_ele, $local_max_index) = ($ele, $i) if $direction >= 0
            && ($local_max_ele eq "" || $ele > $local_max_ele);
        ($local_min_ele, $local_min_index) = ($ele, $i) if $direction <= 0
            && ($local_min_ele eq "" || $ele < $local_min_ele);
    }
    if (!$SEG[$i] && defined $prev_turn_index) { # not at segment start
      check_next:
        $prev_turn_sec = $SEC[$prev_turn_index];
        $local_max_sec = $SEC[$local_max_index];
        $local_min_sec = $SEC[$local_min_index];
        my $at_end = !defined $next_i || $SEG[$next_i] || $ELE[$next_i] eq ""; # at track or segment end or just before missing altitude data
        if ($direction >= 0 && $local_max_index > $prev_turn_index
             && ($at_end || $ele ne "" && $ele <= $local_max_ele - ELE_THRESHOLD))
        {    # at end or when actually starting going down
            my $diff_ele = $local_max_ele - $prev_turn_ele;
            if ($diff_ele >= ELE_THRESHOLD) {
                if ($phases && $diff_ele >= $PHASES_REPORT_THRESHOLD) {
                    push_phase($prev_turn_index, $local_max_index);
                } else {
                    $ascent__phases_suppressed++;
                }
                $SUM__ASCENT[$all ] += $diff_ele;
                $SUM__ASCENT[$part] += $diff_ele if defined $part;
                if (defined $local_max_sec && defined $prev_turn_sec) {
                    my $diff_sec = $local_max_sec - $prev_turn_sec;
                    $TIME__ASCENT[$all ] += $diff_sec;
                    $TIME__ASCENT[$part] += $diff_sec if defined $part;
                } else {
                    $MISSING__ASCENT_TIME[$all ]++;
                    $MISSING__ASCENT_TIME[$part]++ if defined $part;
                }
                $LAST_GAIN_INDEX[$all ] = $local_max_index;
                $LAST_GAIN_INDEX[$part] = $local_max_index if defined $part;
            }
            ($prev_turn_ele, $prev_turn_index) = ($local_max_ele, $local_max_index);
            ($local_min_ele, $local_min_index) = ($ele, $i);
            $direction = -1;
            goto check_next;
        } elsif ($direction <= 0 && $local_min_index > $prev_turn_index
            && ($at_end || $ele ne "" && $ele >= $local_min_ele + ELE_THRESHOLD))
        {    #  at end or when actually starting going up
            my $diff_ele = $local_min_ele - $prev_turn_ele;
            if (-$diff_ele >= ELE_THRESHOLD) {
                if ($phases && -$diff_ele >= $PHASES_REPORT_THRESHOLD) {
                    push_phase($prev_turn_index, $local_min_index);
                } else {
                    $descent_phases_suppressed++;
                }
                $SUM_DESCENT[$all ] += $diff_ele;
                $SUM_DESCENT[$part] += $diff_ele if defined $part;
                my $diff_sec = diff_defined($local_min_sec, $prev_turn_sec);
                if (defined $diff_sec) {
                    $TIME_DESCENT[$all ] += $diff_sec;
                    $TIME_DESCENT[$part] += $diff_sec if defined $part;
                } else {
                    $MISSING_DESCENT_TIME[$all ]++;
                    $MISSING_DESCENT_TIME[$part]++ if defined $part;
                }
                if ($ele ne "") {
                    $LAST_LOSS_INDEX[$all ] = $local_min_index;
                    $LAST_LOSS_INDEX[$part] = $local_min_index if defined $part;
                }
            }
            ($prev_turn_ele, $prev_turn_index) = ($local_min_ele, $local_min_index);
            ($local_max_ele, $local_max_index) = ($ele, $i);
            $direction = 1;
            goto check_next;
        }
    }
    if (!$SEG[$i]) { # not at segment start, implies defined $prev_i, but ignoring merge start/pause
        if ($debug) {
            (my $diff_ele_unused, my $td, my $rate, my $dis, my $comp_spd_unused, my $acc_unused) = comp_diffs($i);
            my $rate_str = $rate eq "" ? "" : sprintf "%.3f", $rate;
            debug_log "phase analysis   at ".point_str($i)." dis=".(sprintf "%.3f", $dis)." rate=$rate_str spd=".spd_or_none($spd);
            my $msg = "internal inconsistency on ";
            abort_trkpt($i, "$msg time difference: $timediff vs. $td at") if $timediff ne $td;
            abort_trkpt($i, "$msg distance value: $DIS[$i] vs. $dis at" ) if abs($DIS[$i] -  $dis) > 0.001;
        }

        if ($weight) {
            unless (defined $sec && defined $prev_sec) {
                warn_trkpt(defined $sec ? PREV($i) : $i, "missing time   at", $cannot_energy);
                $num_energy_missing++;
            } elsif (section($i)) {
                warn_trkpt($i, ($FLG[$i] & START_MERGE ? "start of merge" : "end of merge  ")." at", $cannot_energy);
                $num_energy_missing++;
            } else {
                my $start_sec = $SEC[0] // 0; # used just for smaller time arguments to sub energy()
                $sum_energy += energy($i, $prev_sec - $start_sec, $sec - $start_sec,
                                      $dis, $ELE[$prev_i], $ele, $SPD[$prev_i], $spd);
            }
        }
    }

    # perform general statistics, optionally also for parts (segments or days)
    $min_lat = $lat if !(defined $prev_i) || $min_lat > $lat;
    $max_lat = $lat if !(defined $prev_i) || $max_lat < $lat;
    $min_lon = $lon if !(defined $prev_i) || $min_lon > $lon;
    $max_lon = $lon if !(defined $prev_i) || $max_lon < $lon;
    # calculate min and max tim/sec values rather efficiently asumming ascending order:
    $min_tim = $tim if $min_tim eq ""   ; # || defined $sec && $sec < $min_sec;
    $min_sec = $sec if !defined $min_sec; # || defined $sec && $sec < $min_sec;
    $max_tim = $tim if $tim ne ""  ; # if $max_tim eq ""    || defined $sec && $sec > $max_sec;
    $max_sec = $sec if defined $sec; # if !defined $max_sec || defined $sec && $sec > $max_sec;
    sub set_min_or_max {
        my ($M, $M_IDX, $val, $i, $all, $part, $min) = @_;
        ($$M[$all ], $$M_IDX[$all ]) = ($val , $i) if !defined $$M[$all ] || ($min ? $val < $$M[$all ] : $val > $$M[$all ]);
        return unless defined $part;
        ($$M[$part], $$M_IDX[$part]) = ($val , $i) if !defined $$M[$part] || ($min ? $val < $$M[$part] : $val > $$M[$part]);
    }
    sub set_min { push @_, 1; set_min_or_max(@_); }
    sub set_max { push @_, 0; set_min_or_max(@_); }
    if ($ele ne "") {
        set_min(\@MIN_ELE, \@MIN_ELE_INDEX, $ele, $i, $all, $part);
        set_max(\@MAX_ELE, \@MAX_ELE_INDEX, $ele, $i, $all, $part);
    }
    if ($spd ne "") { # may be meaningful also at segment start: recorded speed
        set_min_or_max(\@MAX_SPD, \@MAX_SPD_INDEX, $spd, $i, $all, $part, 0);
    }
    if (!$SEG[$i]) {
        $SUM_DIS[$all ] += $dis;
        $SUM_DIS[$part] += $dis if defined $part;
        if ($spd ne "" && $timediff ne "" && $spd >= MIN_SPEED_MOVING / MPS2KMPH) {
            $SUM_DIS_MOV     [$all ] += $dis;
            $SUM_DIS_MOV     [$part] += $dis if defined $part;
            $SUM_TIMEDIFF_MOV[$all ] += $timediff;
            $SUM_TIMEDIFF_MOV[$part] += $timediff if defined $part;
        }
        my $slope = avg_slope($i);
        if (defined $slope) {
            set_max(\@MAX__ASC, \@MAX__ASC_INDEX, $slope, $i, $all, $part) if $slope > 0;
            set_min(\@MAX_DESC, \@MAX_DESC_INDEX, $slope, $i, $all, $part) if $slope < 0;
        }
        my $rate = rate(diff_nonempty($ele, $ELE[$prev_i]), $timediff);
        if ($rate ne "" && !section($i)) {
            set_max(\@MAX_GAIN, \@MAX_GAIN_INDEX, $rate, $i, $all, $part) if $rate >= 0;
            set_min(\@MAX_LOSS, \@MAX_LOSS_INDEX, $rate, $i, $all, $part) if $rate <= 0;
        }
    }
}
$lat_prec_len = max(length(lat_str($max_lat)), length(lat_str($min_lat)));
$lon_prec_len = max(length(lon_str($max_lon)), length(lon_str($min_lon)));

#$sec_diff = $SEC[$max_ele_index] - $SEC[$min_ele_index];
#$avg_gain = rate($max_ele - $min_ele, $sec_diff);


## optional elevation correction by geoid height

sub neg_geoid_height { # potentially $updates $found_corr and $ele_corr
    my ($lat, $lon, $ele, $tim) = @_;
    return $ele_corr if !GEOID_ELE_CORRECTION;
    (my $fh, my $tmp_GeoidEval) = tempfile();
    system "wget --quiet -O $tmp_GeoidEval "
      . "http://geographiclib.sourceforge.net/cgi-bin/GeoidEval?input=$lat+$lon";
    open G, $tmp_GeoidEval || abort("FATAL: cannot open geoid data file $tmp_GeoidEval: $.");

    while (<G>) {
        if (m#EGM84</a>\s*=\s*<font.*?>\s*(-?[\.\d]+)\s*</font>#i) {
            $found_corr = 1;
            $ele_corr   = -$1;
        }
    }
    close G;
    if ($found_corr) {
        info_trkpt(-1, "   at    ", "geoid height is $ele_corr m", @_);
    } else {
        warn_trkpt(-1, "at    ", "geoid height is assumed $ele_corr m", @_);
    }
    return $ele_corr;
}

sub linear_ele_correction {
    return 0;                                   # linear correction not wanted
    my $sec  = $_[0];
    my $sec1 = 1296543617;                      #start time
    my $sec2 = 1296569335;                      #end   time
    return 0 unless ($sec1 <= $sec && $sec <= $sec2);   # linear correction not in scope
    my $ele1 = 4720 - 4677;                     #actual and wanted start value
    my $ele2 = 1925 - 1635;                     #actual and wanted end   value
    return
      int(($sec - $sec1) / ($sec2 - $sec1) * ($ele2 - $ele1) + $ele1);
}

sub corrected_ele {
    my ($ele, $sec) = @_;
    return $ele + $ele_corr + linear_ele_correction($sec) if $ele_corr && $ele ne "" && defined $sec;
    return $ele;
}

sub correct_ele { # may update $ele_corr; uses $ele_corr
    my $j = $_[0]; # index of waypoint if >= 0
    my $i = $_[1]; # otherwise, index of trackpoint or < 0 for none
    return "" if $j < 0 && $i < 0;
    my ($lat, $lon, $ele, $sec, $tim) =
      $j >= 0
      ? ($WLAT[$j], $WLON[$j], $WELE[$j], $WTIM[$j], $WSEC[$j])
      : ( $LAT[$i],  $LON[$i],  $ELE[$i],  $TIM[$j], $SEC[$i]);
    return "" if $ele eq "";
    my $prev_i = PREV($i);
    $ele_corr = neg_geoid_height($lat, $lon, $ele, $tim) if ($j >= 0 || !defined $prev_i || defined $SEC[$prev_i]) ;
    # do expensive operation only at specific points or beginning of track segment

    return corrected_ele($ele, $sec);
}


## prepare output header

# compensate for any missing essential start tags
$HEAD =~ s/(<\??xml.*?>\r?\n?)?/($1 ? $1 : "")."<gpx>\n"/se unless $HEAD =~ m/<gpx[ >]/;
$HEAD = $HEAD."<trkseg>\n" unless $HEAD =~ m/<trkseg>/;
$HEAD =~ s/(<trkseg>)/<trk>\n$1/s unless $HEAD =~ m/<trk>/;

$HEAD =~ s|(<gpx)([ >])|$1 version="1.1"$2|      unless $HEAD =~ m|<gpx.*? version\s*=.*?>|s;
$HEAD =~ s|(<gpx)([ >])|$1 creator="$TOOL_ID"$2| unless $HEAD =~ m|<gpx.*? creator\s*=.*?>|s;
# $HEAD =~ s|(<gpx[ >])|$1xmlns:x="http://www.garmin.com/xmlschemas/TrackPointExtension/v2" |; # would not fit for our simple speed extension, and if so, should not be duplicated
#$speed_extension = '<xsd:element name="speed" type="xsd:decimal" minOccurs="0"><xsd:annotation><xsd:documentation>The current speed relative to the previous point, in meters per second</xsd:documentation></xsd:annotation></xsd:element>';
#$HEAD =~ s|(( )*</metadata>)|$2$2<extensions>$speed_extension</extensions>\n$1|s if $speed_extension ne ""; # add declaration for speed extension

$HEAD =~ s|(<gpx.*?>)(\s*)|$1$2<metadata></metadata>$2|s unless $HEAD =~ m|<gpx.*?>\s*<metadata>.*?</metadata>|s;
if ($#NAME >= 0) {
    $HEAD =~ s|(<metadata>)(\s*)|$1$2<name></name>$2|s unless $HEAD =~ m|<metadata>\s*<name>.*?</name>.*?</metadata>|s;
    my $names = "<![CDATA[".join("\n", @NAME)."]]>";
    $HEAD =~ s|(<metadata>\s*<name>).*?(</name>)|$1$names$2|s;
}
if ($HEAD =~ m|(\r?\n?[ \t]*)</metadata>|s) { # match multiple lines with 's' option
    my $indent = $1;

    my $bounds = '<bounds'
      .  ' minlat="' . lat_str($min_lat)
      . '" minlon="' . lon_str($min_lon)
      . '" maxlat="' . lat_str($max_lat)
      . '" maxlon="' . lon_str($max_lon) . '"/>';
    $HEAD =~ s|\r?\n?[ \t]*<bounds .*?/>||sg; # remove any pre-existing bounds entry
    abort("FATAL: cannot add <bounds> to header")
        unless $HEAD =~ s|(\r?\n?[ \t]*<extensions>)|$indent$bounds$1|s
        ||     $HEAD =~ s|(\r?\n?[ \t]*</metadata>)|$indent$bounds$1|s;

    if ($min_tim ne "") {
        $min_tim = round_tim($min_tim);
        $HEAD =~ s|\r?\n?[ \t]*<time>.*?</time>||sg; # remove any pre-existing time entry
        abort("FATAL: cannot add <time> to header")
            unless $HEAD =~ s|(\r?\n?[ \t]*<keywords>)|$indent<time>$min_tim</time>$1|s
            ||     $HEAD =~ s|(\r?\n?[ \t]*<bounds[ >])|$indent<time>$min_tim</time>$1|s;
    }
} else {
    abort("FATAL: cannot find </metadata>, which is usually part of header");
}
$HEAD =~ s/(\r?\n)*[ \t]*(<trk>)/$1$2/sg; # remove empty lines (after metadata) and indentation before <trk>
$HEAD =~ s/ *$//;

# statistics must be generated before interpolating, which may shift indexes

my @STAT; # collecting texts for all statistics points
add_stat_point($PART_START_INDEX[$all], "start", "", EMOJI_PREFIX."▶️");

sub wpt {
    my $j = $_[0]; # index of waypoint
    my ($ele, $tim) = (correct_ele($j, -1), $WTIM[$j]);
    return '<wpt lat="' . lat_str($WLAT[$j]) . '" ' . 'lon="' . lon_str($WLON[$j]) . '">'
        .($ele eq "" ? "     " .(" " x $ele_prec_len) ."      "  : "<ele>" .ele_string($ele_prec_len, $ele) . "</ele>")
        .($tim eq "" ? "      "."                    "."       " : "<time>".round_tim($tim)                 ."</time>")
      . $WTXT[$j] . "</wpt>\n";
}

sub point {
    my ($i, $text) = @_;
    my $ele      = correct_ele(-1, $i);
    my $ele_val  = $ele     eq "" ? "NONE".(" " x max(0, $ele_prec_len - length("NONE"))) : ele_string($ele_prec_len, $ele);
    my $time_val = $TIM[$i] eq "" ? "NONE".(" " x max(0,      $tim_len - length("NONE"))) : round_tim($TIM[$i]);
    return "At lat,lon,ele = " . lat_string($LAT[$i]) . "," . lon_string($LON[$i]).",$ele_val time = $time_val: $text\n";
}

sub add_stat_point { # generate trackpoint-related statistics output and waypoint with same data
    my ($i, $txt, $extra, $sym, $part) = @_; # $i is index of trackpoint containing the location and time
    return unless defined $i;

    my $point_txt = point($i, $txt);
    $STAT[$all ] .= $point_txt;
    $STAT[$part] .= $point_txt if defined $part && $part != $all;
    add_stat_wpt(@_) if !$lean_stat_wpts || defined $extra;
}

sub add_stat_wpt {
    my ($i, $txt, $extra, $sym, $lat, $lon, $ele, $tim, $sec) = @_;
    # unless $i < 0, it is the index of trackpoint containing the location and time:
    ($lat    , $lon    , $ele    , $tim    , $sec    ) =
    ($LAT[$i], $LON[$i], $ELE[$i], $TIM[$i], $SEC[$i]) if $i >= 0;
    stacktrace() unless defined $tim;

    # adapt text to be used as name in statistics wpt
    $txt .= $extra if defined $extra;
    $txt =~ s/\n/; /g;
    $txt =~ s/; $//;
    $txt =~ s/\s+/ /g;
    $txt = "$tim $txt" if !$lean_stat_wpts && $tim ne "";
    $txt = "<name>[$txt]</name>" if $txt ne "";
    $txt .= "<src>".TOOL_NAME."</src>";
    $txt .= "<sym>$sym</sym>" if defined $sym;

    # adapt to avoid rouding issues on comparison
    $lat = lat_str($lat);
    $lon = lon_str($lon);
    $ele = ele_string($ele_prec_len, $ele) if $ele ne "";
    $tim = round_tim($tim) if $tim ne "";

    # insert wpt in existing array, avoiding duplicates,
    # in ascending time order after any entries with same time, where missing time is considered infinite,
    # assuming that the array is already ordered

    my $sec_or_inf = $sec // $MAX_UNDEF_TIME;
    my $j = 0;
    while ($j <= $#WSEC && ($WSEC[$j] // $MAX_UNDEF_TIME) < $sec_or_inf) {
        $j++;
    }
    while ($j <= $#WSEC && ($WSEC[$j] // $MAX_UNDEF_TIME) == $sec_or_inf
           && !($txt =~ m/(Z |\[)start\]/) # make sure that "start" wpt is placed first
           && ($txt =~ m/(Z |\[)end; / || !($WTXT[$j] =~ m/(Z |\[)end; /))) { # make sure that "end" wpt is placed last
        return if # do not add duplicates before point of insertion
            $WLAT[$j] == $lat && $WLON[$j] == $lon && eq_nonempty ($WELE[$j], $ele) &&
            $WTIM[$j] eq $tim && $WTXT[$j] eq $txt;
        $j++;
    }
    for (my $k = $j; $k <= $#WSEC; $k++) {
        return if # do not add duplicates at point of insertion or after it
            $WLAT[$k] == $lat && $WLON[$k] == $lon && eq_nonempty ($WELE[$k], $ele) &&
            $WTIM[$k] eq $tim && $WTXT[$k] eq $txt;
    }
    splice @WLAT, $j, 0, $lat;
    splice @WLON, $j, 0, $lon;
    splice @WELE, $j, 0, $ele;
    splice @WTIM, $j, 0, $tim;
    splice @WSEC, $j, 0, $sec;
    splice @WTXT, $j, 0, $txt;
    splice @WSTR, $j, 0, undef;
    splice @WIDX, $j, 0, $i;
    $FLG[$i] |= STATISTICS_WPT if $i >= 0;

}

# add overall statistics and for parts

my $part_last = $part_end // $part // $all;
sub part_number {
    my ($char, $i) = @_;
    return ($char x (length($part_last + $part_offset) - length($i + $part_offset))).($i + $part_offset);
}

for (my $overall_stat, my $p = $all; $p <= $part_last; $p++) {
    $STAT[$p] = $overall_stat if $p != $all; # prepend overall statistics to all parts

    my $Part = $p == $all ? "" : ($segs ? "segment" : "day")." ".part_number(' ', $p)." ";
    my $Min_ele      = defined $MIN_ELE_INDEX[$p] ? sprintf ELE_PRECISION, correct_ele(-1, $MIN_ELE_INDEX[$p]) : "";
    my $Max_ele      = defined $MAX_ELE_INDEX[$p] ? sprintf ELE_PRECISION, correct_ele(-1, $MAX_ELE_INDEX[$p]) : "";
    my $Max_spd      = defined $MAX_SPD_INDEX [$p] ? " ".spd_full($SPD[$MAX_SPD_INDEX [$p]]) : "";
    my $Max_gain_spd = defined $MAX_GAIN_INDEX[$p] ?     spd_str ($SPD[$MAX_GAIN_INDEX[$p]]) : "";
    my $Max_loss_spd = defined $MAX_LOSS_INDEX[$p] ?     spd_str ($SPD[$MAX_LOSS_INDEX[$p]]) : "";
    while (substr($Max_gain_spd, 0, 1) eq " " &&
           substr($Max_loss_spd, 0, 1) eq " ")
    {
        $Max_gain_spd = substr($Max_gain_spd, 1);
        $Max_loss_spd = substr($Max_loss_spd, 1);
    }

    # make sure to report also non-presence of max ascent & descent rates
    $MAX__ASC_INDEX[$p] //= $PART_START_INDEX[$p];
    # print STDERR "#### slope max_desc ($MAX_DESC[$p], $MAX_DESC_INDEX[$p]) = (0, $PART_END_INDEX[$p])\n" if !defined $MAX_DESC_INDEX[$p] && defined $MAX__ASC_INDEX[$p]; 
    $MAX_DESC_INDEX[$p] //= $PART_START_INDEX[$p];

    sub Avg_slope { return defined $_[0] ? (sprintf SLOPE_PRECISION, $_[0])." degrees" : "none"; }
    my $Max__asc =  Avg_slope($MAX__ASC[$p], $MAX__ASC_INDEX[$p]);
    my $Max_desc = (Avg_slope($MAX_DESC[$p], $MAX_DESC_INDEX[$p])) =~ s/\+0/0/r;

    # make sure to report max ascent & descent rates always pairwise
    ($MAX_GAIN[$p], $MAX_GAIN_INDEX[$p]) = (0, $PART_END_INDEX[$p]) if !defined $MAX_GAIN_INDEX[$p] && defined $MAX_LOSS_INDEX[$p];
    ($MAX_LOSS[$p], $MAX_LOSS_INDEX[$p]) = (0, $PART_END_INDEX[$p]) if !defined $MAX_LOSS_INDEX[$p] && defined $MAX_GAIN_INDEX[$p];
    my $Max_gain = sprintf RAT_PRECISION, $MAX_GAIN[$p] // 0; # relevant only if defined $MAX_GAIN_INDEX[$p]
    my $Max_loss = sprintf RAT_PRECISION, $MAX_LOSS[$p] // 0; # relevant only if defined $MAX_LOSS_INDEX[$p]
    $Max_loss =~ s/\+0/ 0/;

    #my $Avg_gain=sprintf RAT_PRECISION, $avg_gain;
    sub sum_avg_rate {
        my ($sum, $time, $miss) = @_;
        my $avg_rate = rate($sum, $time);
        my $Avg_rate = $avg_rate ne "" ? " = ".(sprintf RAT_PRECISION, $avg_rate)." m/h"
            . ($miss ? " with ".points($miss)." not having time information" : "")
            : "unknown due to missing time";
        return (sprintf DIF_PRECISION, $sum // 0)." m; avg rate $Avg_rate";
    }
    my $text__ascent = sum_avg_rate($SUM__ASCENT[$p], $TIME__ASCENT[$p], $MISSING__ASCENT_TIME[$p]);
    my $text_descent = sum_avg_rate($SUM_DESCENT[$p], $TIME_DESCENT[$p], $MISSING_DESCENT_TIME[$p]);

    my $max_spd_text  = "max speed        = $Max_spd km/h"; # may use "" if not defined
    my $min_ele_text  = "min altitude     = $Min_ele m";    # may use "" if not defined
    my $max_ele_text  = "max altitude     = $Max_ele m";    # may use "" if not defined
    my $ascent_text   = "total  ascent    = $text__ascent";
    my $descent_text  = "total descent    = $text_descent";
    my $max__asc_text = "max rolling average  ascent (slope) = $Max__asc";
    my $max_desc_text = "max rolling average descent (slope) = $Max_desc";
    my $max_gain_text = "max  ascent rate = $Max_gain m/h".
        ($Max_gain_spd eq "" ? "" : " at $Max_gain_spd km/h"); # relevant only if defined $MAX_GAIN_INDEX[$p]
    my $max_loss_text = "max descent rate = $Max_loss m/h".
        ($Max_loss_spd eq "" ? "" : " at $Max_loss_spd km/h"); # relevant only if defined $MAX_LOSS_INDEX[$p]
    my $part_end_text = $Part."end";
    my $sep = $p == $all ? "\ntotal" : ";";
    $part_end_text .= "$sep length ".($p == $all ? "              ": "     ")."= ".dis_string($SUM_DIS[$p]) if defined $SUM_DIS[$p];
    if ($p == $all) {
        $part_end_text .= "\ngross duration             = ".timediff_string(diff_defined($max_sec, $min_sec) // 0) if defined $SUM_DIS[$all];
    } else {
        $part_end_text .= "$sep duration = ".timediff_string($DURATION[$p]) if $DURATION[$p];
    }
    if ($SUM_TIMEDIFF_MOV[$p]) {
        $part_end_text .= "$sep moving time ".($p == $all ? "         ": "")."= ".
            timediff_string($SUM_TIMEDIFF_MOV[$p]); # ", with speed being at least ${\MIN_SPEED_MOVING} km/h\n";
        $part_end_text .= "$sep average moving speed = ".
            (sprintf AVG_PRECISION, MPS2KMPH * $SUM_DIS_MOV[$p] / $SUM_TIMEDIFF_MOV[$p])." km/h";
    }
    if ($p == $all && $weight) {
        $part_end_text .= "\ntotal biking energy        = ".(sprintf "%.0f", $sum_energy / 1000)." kJ";
        $part_end_text .= ", leaving out ".points($num_energy_missing)." where time or altitude information is missing or at start or end of merge sections" if $num_energy_missing;
    }
    add_stat_point($MAX_SPD_INDEX  [$p], $Part.$max_spd_text , undef, EMOJI_PREFIX."⏩", $p);
    add_stat_point($MIN_ELE_INDEX  [$p], $Part.$min_ele_text , undef, EMOJI_PREFIX."🔽", $p);
    add_stat_point($MAX_ELE_INDEX  [$p], $Part.$max_ele_text , undef, EMOJI_PREFIX."🔼", $p);
    add_stat_point($LAST_GAIN_INDEX[$p], $Part.$ascent_text  , undef, EMOJI_PREFIX."⬆️", $p);
    add_stat_point($LAST_LOSS_INDEX[$p], $Part.$descent_text , undef, EMOJI_PREFIX."⬇️", $p);
    add_stat_point($MAX__ASC_INDEX [$p], $Part.$max__asc_text, undef, EMOJI_PREFIX."↗️", $p);
    add_stat_point($MAX_DESC_INDEX [$p], $Part.$max_desc_text, undef, EMOJI_PREFIX."↘️", $p);
    add_stat_point($MAX_GAIN_INDEX [$p], $Part.$max_gain_text, undef, EMOJI_PREFIX."⏫", $p);
    add_stat_point($MAX_LOSS_INDEX [$p], $Part.$max_loss_text, undef, EMOJI_PREFIX."⏬", $p);
    $LAT[$PART_END_INDEX[$all]] += 0.00001 if $p != $all && $p == $part_last; # make sure that last "part end" wpt is placed separately
    my $extra = $lean_stat_wpts ? ("; $max_spd_text; $min_ele_text; $max_ele_text") =~ s/\s+/ /gr : "";
    add_stat_point($PART_END_INDEX [$p],       $part_end_text, $extra, EMOJI_PREFIX.($p == $all ? "⏹️" : $days ? "😴" : "🈁"), $p);
    # missing a "pause" symbol or the like for segment ends etc., was using "🎬" before instead of "🈁"
    $LAT[$PART_END_INDEX[$all]] -= 0.00001 if $p != $all && $p == $part_last; # make sure that last "part end" wpt is placed separately
    $overall_stat = $STAT[$all] if $p == $all;
}
#average ".($avg_gain >= 0 ?  "ascent rate between lowest and highest"
#	                  : "descent rate between highest and lowest")." point = $Avg_gain m/h

# add statistics for phases

sub order_postfix {
    my $n = $_[0];
    my $d = $n % 10;
    return 11 <= $n && $n <= 13 ? "th" :
        $d == 1 ? "st" :
        $d == 2 ? "nd" :
        $d == 3 ? "rd" :
        "th";
}
sub phase_text {
    my $n = $_[0] + 1;
    my $d = $n % 10;
    return (" " x max(0, length($#PHASE_END_INDEX + 1) - length($n)))
        .$n.order_postfix($n)." ";
}

for (my $p = 0; $p <= $#PHASE_END_INDEX; $p++) {
    my $type = $PHASE_DIFF[$p] >= 0 ? "ascent " : "descent";
    my $Len = phase_text($p) . "$type   = " . (sprintf DIF_PRECISION, $PHASE_DIFF[$p]) . "  m;   "
      . "distance   =".($PHASE_DIST[$p] / 1000 >= 10 ? "" : " ").dis_string($PHASE_DIST[$p]);
    my $duration = $PHASE_DURATION[$p];
    $Len .= "; duration = " . timediff_string($duration) if $duration;
    my $max_spd_idx = $PHASE_MAX_SPD_INDEX[$p];
    my $Spd = phase_text($p)
        . "max speed =  ".spd_full($SPD[$max_spd_idx])." km/h; "
        . "avg speed  = ".spd_full($PHASE_SPD[    $p])." km/h"
        if defined $max_spd_idx;
    my $max_rate_idx = $PHASE_MAX_RATE_INDEX[$p];
    my $Rate = phase_text($p)
        . "max rate  =" . (sprintf RAT_PRECISION2, $PHASE_MAX_RATE[$p]) . "  m/h; "
        . "avg rate   =" . (sprintf RAT_PRECISION, $PHASE_RATE    [$p]) . "  m/h"
        if defined $max_rate_idx;
    add_stat_point($PHASE_END_INDEX     [$p], $Len );
    add_stat_point($PHASE_MAX_SPD_INDEX [$p], $Spd ) if defined $max_spd_idx;
    add_stat_point($PHASE_MAX_RATE_INDEX[$p], $Rate) if defined $max_rate_idx;
}


## fill long time gaps as far as needed for exiftool etc. by inserting intermediate points with linear interpolation of data

my ($num_inserted_trkpts, $num_duplicated_trkpts) = (0, 0);
sub insert_points_in_gaps { # only done virtually, as string to be printed before the respective points
    my ($prev_i, $i);
    for ($i = $PART_START_INDEX[$all];
         defined $i;
         $prev_i = $i, $i = NEXT($i)) {
        # check_PREV_NEXT($i);
        next if !defined $prev_i;

        my ( $lat    , $lon    , $ele    , $tim    ) =
            ($LAT[$i], $LON[$i], $ELE[$i], $TIM[$i]);
        # do also across segments:
        # my $sq_dis2d = $DIS2[$i];
        my $dis = $DIS[$i];
        my $diff_time = diff_defined($SEC[$i], $SEC[$prev_i]) // 0;
        my $dis_text = defined $dis ? dis_str(DIS_PRECISION2, $dis) : "with no distance available";
        debug_log "considering  gap at ".point_str($i)." of $diff_time seconds and $dis_text";
        if (int($diff_time) > INSERTION_MAX_GAP) {
            # note that this is done only between adjacent track points that both have time information
            my $lin_spd = $dis / $diff_time;
            my $gap_text = "gap = ".timediff_str("%2d", $diff_time)." over $dis_text";
            if (MPS2KMPH * $lin_spd > $INSERTION_MAX_SPEED) {
                warn_trkpt($i, "before", "$gap_text; not inserting points because linear speed would be ".
                           spd_str($lin_spd)." km/h, higher than $INSERTION_MAX_SPEED");
                next;
            }
            # insert equidistant points with interpolated data, plus
            # in case of new segment: an extra trkpt that equals the end of the previous segment and one that equals the start of the next segment
            my $n_spans  = int($diff_time / INSERTION_MAX_GAP) + ($diff_time % INSERTION_MAX_GAP == 0 ? 0 : 1);
            my $seg = INSERT_SEGMENT ? $SEG[$i] : 0; # if gap was at segment start, optionally have segment starts just before and after inserted points
            my $start = $seg ? 0 : 1;
            my $n_end = $n_spans - ($seg ? 0 : 1);

            my $p = $prev_i;
            my $flags_before = $FLG[$p];
            my $flags_after  = $FLG[$i];
            # my $sq_dis2d_delta  = $sq_dis2d / $n_spans / $n_spans;
            # my      $dis_delta  =    $dis   / $n_spans;
            my $lat_delta = ($lat - $LAT[$p]) / $n_spans;
            my $lon_delta = ($lon - $LON[$p]) / $n_spans;
            my $sec_delta = $diff_time / $n_spans; # 0 < $sec_delta <= INSERTION_MAX_GAP
            my $ele_delta = diff_nonempty($ele, $ELE[$p]);
            $ele_delta /= $n_spans if defined $ele_delta;
            info_trkpt($i, "before", "$gap_text; inserting ".points($n_spans - 1)." each ".timediff_string($sec_delta));

            # point $i will become the point past the insertion (possibly with duplication before and after it)
            # need to move info about previous ignored points from point past the insertion to first new point
            my $ign = $IGN[$i];
            # for the point that will be the one past the insertion:
            $IGN[$i] = 0;
            # ($DIS2[$i], $DIS[$i]) = $seg ? (0, 0) : ($sq_dis2d_delta, $dis_delta); # recalculate distance
            # ascent/decent rate as well as speed stays the same, due to linear interpolation

            # https://www.informit.com/articles/article.aspx?p=1582932&seqNum=4
            # You can build up a string with a filehandle
            open my ($str_fh), '>>', \my $inserted_points;
            my ( $lat    , $lon    , $ele    , $sec     ) =
                ($LAT[$p], $LON[$p], $ELE[$p], $SEC[$p]);
            for (my $n = 0; $n <= $n_end; $n++) {
                if ($n >= $start) {
                    $tim = epoch_to_str($sec);
                    my $flags;
                    if ($n == 0) {
                        # place a new segment start, with contents to be duplicated from previous point (at segment end just before)
                        $flags = $flags_before | DUPLICATED_POINT;
                        $num_duplicated_trkpts++;
                    } elsif ($n == $n_spans) {
                        $flags = $flags_after | DUPLICATED_POINT;
                        # place a new segment end, with contents duplicated from next point at segment start just after insertion
                        $num_duplicated_trkpts++;
                    } else {
                        $flags = INSERTED_POINT;
                        $num_inserted_trkpts++;
                    }
                    my $cmt = "";
                    $cmt = plural($ign, "previous point")." ignored" if !DISABLE_TRKPT_CMT && $n == $start && $ign != 0;
                    my $src = "";
                    $src = $flags & INSERTED_POINT
                        ? $bride_gap_insert_txt
                        : $bride_gap_duplicate_txt .(defined $prev_i && $FLG[$prev_i] & INSERTED_POINT ? "after" : "before")
                        if !DISABLE_TRKPT_SRC;
                    print_point($str_fh, $n == 0, $lat, $lon, corrected_ele($ele, $sec),
                                $tim, $cmt, $src, DISABLE_TRKPT_FIX ? "" : "none", "");
                    # TODO check; was before: $fix = $flags == 0 ? "" : ($FLG[$i] & (INSERTED_POINT | INTERPOLATED_TIM | SUBSTITUTED_TIM) == 0 ? "2d" : "none");
                    add_stat_wpt(-1, "INFO: "    . $cmt, undef, undef, $lat, $lon, $ele, $tim, $sec)
                        if $cmt_wpts && $cmt && $cmt =~ m|$cmt_wpts|;
                    add_stat_wpt(-1, "WARNING: " . $src, undef, undef, $lat, $lon, $ele, $tim, $sec)
                        if $src_wpts && $src && $src =~ m|$src_wpts|;
                }
                $lat += $lat_delta;
                $lon += $lon_delta;
                $ele += $ele_delta if defined $ele_delta;
                $sec += $sec_delta;
            }
            $TXT[$i] = $inserted_points;
        }
    }
}
insert_points_in_gaps() if $insert_points;


### generate final output

my $final_text = "";
$final_text .=
    "$ascent__phases_suppressed ascent phases and $descent_phases_suppressed descent phases " .
    "not reported because they are below the threshold of +/-$PHASES_REPORT_THRESHOLD m\n"
    if $phases && $PHASES_REPORT_THRESHOLD > ELE_THRESHOLD;
$final_text .= "average deviation between recorded and computed speed (for points at most ".MAX_TIMEDIFF_RECORDED_SPEED." seconds apart) = ".
    int(0.5 + 100 * $sum_speed_deviation / ($sum_speed / 2))." %\n" if $sum_speed;
$final_text .= plural($num_pruned_wpts,         "waypoint")    . " pruned\n" if $num_pruned_wpts;
$final_text .= plural($num_auto_wpts, "automatic waypoint")." pruned\n" if $num_auto_wpts;
$final_text .= plural($num_ignored_hdrs,      "GPX header")." ignored (apart from any included names and waypoints)\n" if $num_ignored_hdrs;
$final_text .= plural($num_pruned_cmts,              "cmt")." pruned\n" if $num_pruned_cmts;
$final_text .= plural($num_pruned_descs,     "description")." pruned\n" if $num_pruned_descs;
$final_text .= plural($num_pruned_links,            "link")." pruned\n" if $num_pruned_links;
$final_text .= plural($num_pruned_exts,        "extension")." pruned\n" if $num_pruned_exts;
$final_text .= plural($num_merges,               "section")." merged from alternative input\n" if $num_merges;
$final_text .= plural($num_merged_points,   "merged point")."\n" if $merge;
$final_text .= plural($num_skip_points,    "regular point")." skipped during merge\n" if $num_skip_points;
$final_text .= points($num_ignored_diff  )." less than $min_input_timediff sec apart filtered out\n" if $num_ignored_diff;
$final_text .= points($num_ignored_before)." with time before ".epoch_to_str($begin_sec)." filtered out\n" if defined $begin_sec;
$final_text .= points($num_ignored_after )." with time after  ".epoch_to_str($end_sec  )." filtered out\n" if defined $end_sec;
$final_text .= points($num_elem_above    )." with $max_elem above $max_elem_limit"      ." filtered out\n" if $max_elem;
$final_text .= points($num_elem_below    )." with $min_elem below $min_elem_limit"      ." filtered out\n" if $min_elem;
$final_text .= points($num_removals_ele  )." ignored due to altitude less than "            .MIN_PLAUSIBLE_ELE                                      ." m\n" if $num_removals_ele;
$final_text .= points($num_removals_dir  )." ignored due to direction change out of range -".MAX_PLAUSIBLE_ANGLE_DIFF."..+".MAX_PLAUSIBLE_ANGLE_DIFF." degrees\n" if $num_removals_dir;
$final_text .= points($num_removals_spd  )." ignored due to speed out of range 0..$MAX_PLAUSIBLE_SPEED km/h\n"    if $num_removals_spd;
$final_text .= points($num_removals_acc  )." ignored due to acceleration more than $MAX_PLAUSIBLE_ACCEL m/s/s\n"  if $num_removals_acc;
$final_text .= points($num_removals_dec  )." ignored due to deceleration more than $MAX_PLAUSIBLE_DECEL m/s/s\n"  if $num_removals_dec;
$final_text .= points($num_removals_gain )." ignored due to ascent rate more than $MAX_PLAUSIBLE_ELE_GAIN m/s\n"  if $num_removals_gain;
$final_text .= points($num_removals_loss )." ignored due to descent rate more than $MAX_PLAUSIBLE_ELE_LOSS m/s\n" if $num_removals_loss;
$final_text .= points($num_removals_diff )." ignored due to time difference less than $MIN_TIMEDIFF sec\n"        if $num_removals_diff;
$final_text .= plural($num_ext_below, "extension element")." with name '$min_ext' and value below $min_ext_limit pruned\n" if $min_ext;
$final_text .= plural($num_ext_above, "extension element")." with name '$max_ext' and value above $max_ext_limit pruned\n" if $max_ext;
$final_text .= points($num_smoothened)." being at most ${ \SMOOTHING_MAX_GAP } seconds apart smoothened\n" if $smoothing;
$final_text .= points($num_inserted_trkpts)." inserted ".(INSERT_SEGMENT ? "and ".points($num_duplicated_trkpts)." duplicated " : "").
    "due to time gaps > ".INSERTION_MAX_GAP." seconds\n" if $insert_points;

my $num_sanitized = $num_removals_ele + $num_removals_dir + $num_removals_spd + $num_removals_acc
    + $num_removals_dec + $num_removals_gain + $num_removals_loss + $num_removals_diff;
my $sanitized_text = $num_sanitized ? plural($num_sanitized, "sanitized trackpoint")." nor " : "";
my $stats_start = TOOL_NAME." statistics (not including $sanitized_text"."periods between segments)\n";
my $stats_end = "elevation ". (!GEOID_ELE_CORRECTION ? "not corrected" :
                               "corrected by " . ($found_corr ? "last retrieved" : "assumed") . " geoid height = $ele_corr m") . "\n"
    . TOOL_NAME." statistics end";

## expand metadata in trackpoints and optionally add waypoints for matching 'cmt' and 'src' elements
for (my $i = $PART_START_INDEX[$all]; defined $i; $i = NEXT($i)) {
    # check_PREV_NEXT($i);

    sub append_unless_present { my ($str, $msg) = @_; return $str =~ m/$msg/ ? "" : ($str eq "" ? "" : "; ") . $msg; }

    if (!DISABLE_TRKPT_CMT) {
        my $cmt = $CMT[$i];
        if ($IGN[$i]) {
            my $old_count = ($cmt =~ s#(\d+) (previous|final) points? ignored(; )?##s ? $1 : 0); # remove and previous instance of our comment
            my $ign_str = plural(($ARGV =~ m/^(test.*?\.out)\d$/ ? 0 : $old_count) + $IGN[$i], ($i ? "previous" : "final")." point")." ignored";
            $cmt = $cmt eq "" ? $ign_str : "$ign_str; $cmt";
        }
        $cmt .= append_unless_present($cmt, "start alternative"    ) if $FLG[$i] & START_MERGE;
        $cmt .= append_unless_present($cmt, "pause alternative"    ) if $FLG[$i] & PAUSE_MERGE;
        $CMT[$i] = $cmt;
        add_stat_wpt($i, "INFO: " . $cmt) if $cmt_wpts && $cmt && $cmt =~ m|$cmt_wpts|;
    }

    if (!DISABLE_TRKPT_SRC) {
        my $src = $SRC[$i];
        # cannot place info on interpolation etc. in <trkseg> element because it cannot hold a comment
        $src .= append_unless_present($src, "interpolated altitude") if $FLG[$i] & INTERPOLATED_ELE;
        $src .= append_unless_present($src, "interpolated time"    ) if $FLG[$i] & INTERPOLATED_TIM;
        $src .= append_unless_present($src, "substituted altitude" ) if $FLG[$i] & SUBSTITUTED_ELE;
        $src .= append_unless_present($src, "substituted time"     ) if $FLG[$i] & SUBSTITUTED_TIM;
        $src .= append_unless_present($src, "computed speed"       ) if $FLG[$i] & COMPUTED_SPEED && INCLUDE_SPEED;
        $src .= append_unless_present($src, "corrected altitude"   ) if $FLG[$i] & CORRECTED_ELE;
        $SRC[$i] = $src;
        add_stat_wpt($i, "WARNING: " . $src) if $src_wpts && $src && $src =~ m|$src_wpts|;
    }
    # as far as needed, metadata on points inserted/duplicated due to time gaps already handled by insert_points_in_gaps()

    my $exts = $EXT[$i];
    $exts =~ s#<gpxtpx:TrackPointExtension>\s*<gpxtpx:speed>.*?</gpxtpx:speed>\s*</gpxtpx:TrackPointExtension>##sg if PRUNE_RECORDED_SPEED;
    if (INCLUDE_SPEED && $SPD[$i] ne "") { # may be recorded or computed, depending on $FLG[$i] & COMPUTED_SPEED
        my $spdstr = "<speed>".spd_string($SPD[$i])."</speed>";
        $exts .= $spdstr # add new speed extension at the end of existing extensions
            unless $exts =~ s|<speed>.*?</speed>|$spdstr|sg; # substitute previous speed extension at same place
    }
    $EXT[$i] = $exts;
}

abort("FATAL: end of track header not found") unless $HEAD =~ m/<trk>(.*?)<trkseg>/s;
# add preliminarily empty comment in track; # if not already existing:
abort("FATAL: end of track header not found") unless $1 =~ m/<cmt>/ || $HEAD =~ s|(<trk>(\s*<name>.*?</name>)?)|$1\n<cmt></cmt>|s;
# HEAD includes <gpx> <trk> <trkseg>


### write final output

sub print_point {
    my ($out, $seg, $lat, $lon, $ele, $tim, $cmt, $src, $fix, $ext) = @_;
        print $out "</trkseg>\n<trkseg>\n" if $seg;
        print $out '<trkpt lat="' . lat_str($lat) . '" lon="' . lon_str($lon) . '">';
        print $out "<ele>"  . ele_string($ele_prec_len, $ele) . "</ele>" if $ele ne "";
        print $out "<time>" . round_tim($tim) . "</time>"                           if $tim;
        print $out "<cmt>$cmt</cmt>"                          if !DISABLE_TRKPT_CMT && $cmt;
        print $out "<src>$src</src>"                          if !DISABLE_TRKPT_SRC && $src;
        print $out "<fix>$fix</fix>"                          if !DISABLE_TRKPT_FIX && $fix;
        print $out "<extensions>$ext</extensions>"                                  if $ext;
        print $out "</trkpt>\n";
}

sub print_part {
    my ($out, $p, $head, $wpts) = @_;
    my $stats = $stats_start.$STAT[$p].$final_text.$stats_end;
    info($stats) if $p == $all;
    $head =~ s/(<trk>.*?<cmt>)/$1\n$stats\n/s; # prepend stats to track comment
    print $out ($head =~ s/(<trk>)/$wpts$1/sr); # prepend collected and computed waypoints to track

    my $recent_dist = 0 if defined $min_dist;
    my $recent_angle = 0 if defined $min_angle;
    my $safetynet_count = $#TIM;
    for (my $i = $PART_START_INDEX[$p]; ; $i = NEXT($i)) {
        # check_PREV_NEXT($i);
        abort "internal error: non-terminating loop on printing trackpoints of part $p (all == $all)" if $safetynet_count-- < 0;
        if ($i != $PART_END_INDEX[$p]) {
            if (defined $min_dist && $i != $PART_START_INDEX[$p] && !$SEG[$i]) {
                # debug_log "$TIM[$i] recent_dist $recent_dist += $DIS[$i]";
                $recent_dist += $DIS[$i];
                if ($recent_dist < $min_dist && !($FLG[$i] & STATISTICS_WPT)) {
                    debug_log "$TIM[$i] removed due to recent_dist $recent_dist";
                    $num_removed_min_dist++;
                    next;
                }
                $recent_dist = 0;
            }
            if (defined $min_angle && $i != $PART_START_INDEX[$p] && !$SEG[$i]) {
                if (defined $ANG[$i]) {
                    # debug_log "$TIM[$i] recent_angle $recent_angle += abs($ANG[$i])";
                    $recent_angle += abs($ANG[$i]);
                    if ($recent_angle < $min_angle && !($FLG[$i] & STATISTICS_WPT)) {
                        debug_log "$TIM[$i] removed due to recent_angle $recent_angle";
                        $num_removed_min_angle++;
                        next;
                    }
                }
                $recent_angle = 0;
            }
        }
        if ($TXT[$i] # points to be inserted virtually for bridging time gap
            && ($p == $all || $i != $PART_START_INDEX[$p])) {
            print $out $TXT[$i];
            $num_points_output += scalar(() = $TXT[$i] =~ /<trkpt /g);
        }
        print_point($out, $p == $all && defined PREV($i) && $SEG[$i], $LAT[$i], $LON[$i],
                    correct_ele(-1, $i), $TIM[$i], $CMT[$i], $SRC[$i], $FIX[$i],$EXT[$i]);
        $num_points_output++;
        last if $i == $PART_END_INDEX[$p];
    }
    print $out "</trkseg>\n</trk>\n</gpx>\n";
    close $out;
}

my $all_wpts = "";
if ($split) {
    my $j = 0; # index for waypoints, which are assumed to be in same (time-wise) order as trackpoints and parts
    for (my $p = $all + 1; $p <= $part_last; $p++) {
        my $part_wpts = "";
        my $p_start_sec = $SEC[$PART_START_INDEX[$p]];
        while ($j <= $#WTIM && (!defined $WSEC[$j] || !defined $p_start_sec || $WSEC[$j] < $p_start_sec)) {
            last if $WIDX[$j] == $PART_START_INDEX[$p]; # fallback independent of availability of time info
            $j++;
        }
        my $p_end_sec = $SEC[$PART_END_INDEX[$p]];
        my $found_same_idx = 0;
        while ($j <= $#WTIM && (!defined $WSEC[$j] || !defined $p_end_sec || $WSEC[$j] <= $p_end_sec)) {
            if ($WIDX[$j] == $PART_END_INDEX[$p]) { # fallback independent of availability of time info
                $found_same_idx = 1;
            } else {
                last if $found_same_idx;
            }
            my $wpt_str = wpt($j);
            $all_wpts  .= $wpt_str; # assuming that the parts cover all data
            $part_wpts .= $wpt_str if $WIDX[$j] != -1;
            $j++;
        }
        my $part_head = $HEAD;
        my $part_name = part_number('0', $p);
        $part_head =~ s|(<trk>)|$1<name><![CDATA[$split]]></name>|s unless $part_head =~ m|<trk>\s*<name>|s;
        $part_head =~ s|(<trk>.*?)(\]*>?</name>)|$1 $part_name$2|s;

        my $file_name = $split."_$part_name.gpx";
        my $out;
        open($out, "> $file_name") || abort("FATAL: cannot open spilt GPX output file $file_name: $.");
        print_part($out, $p, $part_head, $part_wpts);
    }
} else {
    for (my $j = 0; $j <= $#WTIM; $j++) {
        $all_wpts .= wpt($j);
    }
}

my $all_head = $HEAD;
print_line(""); # empty line before stats
print_part($out, $all, $all_head, $all_wpts);
info("output ".points($num_points_output)
     .(defined $min_dist ? ", removing ".points($num_removed_min_dist)." closer than $min_dist m" : "")
     .(defined $min_angle ? ", removing ".points($num_removed_min_angle)." changing course less than $min_angle°" : ""));

exit 0;