lf2mpr_nrt.pdl

#!/homef/nabil/perl5/perlbrew/perls/perl-5.42.0/bin/perl
use v5.38;
use FindBin qw($RealBin);
use lib "$RealBin/lib";
use AIALimbfit::DrmsRuntime    qw(configure_drms_environment validate_drms_runtime show_info_lines);
use AIALimbfit::Inventory      qw(limb_path);
use AIALimbfit::LimbfitCommand qw(validate_limbfit_args);
use AIALimbfit::Reducer        qw(kwd_lines masterpoint_filename reduce_limb_points);
use AIALimbfit::Slot           qw(tstr2ts);
use Getopt::Long;
use File::Path qw(make_path);
use PDL::IO::Misc;
use Scalar::Util qw(looks_like_number);

my $config_file = $ENV{AIA_LIMBFIT_CONFIG} // "$RealBin/config.pl";
my $cfg         = do $config_file or die "Cannot load $config_file: " . ( $@ || $! );
my ( $yr, $mo, $da, $hr );
my $inpdir = $cfg->{fits_root};
my $outdir = $cfg->{pointing_dir};
my ( $previous, $next );
my @requested_wavelengths;

GetOptions(
  'year=i'                 => \$yr,
  'month=i'                => \$mo,
  'day=i'                  => \$da,
  'hour=i'                 => \$hr,
  'inpdir=s'               => \$inpdir,
  'outdir=s'               => \$outdir,
  'wavel=i@'               => \@requested_wavelengths,
  'interpolate-previous=s' => \$previous,
  'interpolate-next=s'     => \$next,
) or die "Invalid options\n";
die "--year, --month, --day, and --hour are required\n"
  unless defined $yr && defined $mo && defined $da && defined $hr;
validate_limbfit_args( year => $yr, month => $mo, day => $da, hour => $hr );
validate_limbfit_args(
  year        => $yr,
  month       => $mo,
  day         => $da,
  hour        => $hr,
  wavelength  => $_,
  wavelengths => $cfg->{wl}
) for @requested_wavelengths;
die "Hour must be on the $cfg->{cadence_h}-hour cadence\n" if $hr % $cfg->{cadence_h};
die "Input and output directories must not be empty\n"
  unless defined $inpdir && length($inpdir) && defined $outdir && length($outdir);
die "Interpolation requires both --interpolate-previous and --interpolate-next\n"
  if defined($previous) != defined($next);

umask 0002;
my %seen;
my @wl = @requested_wavelengths ? grep { !$seen{$_}++ } @requested_wavelengths : @{ $cfg->{wl} };
my $outbase = masterpoint_filename( $yr, $mo, $da, $hr, $cfg->{cadence_h} );
my $outnam  = "$outdir/$outbase";
my $slot    = sprintf '%04d-%02d-%02dT%02d:00:00Z', $yr, $mo, $da, $hr;
unlink $outnam or die "Cannot remove stale '$outnam': $!\n" if -e $outnam;

my @kwd;
if ( defined $previous ) {
  configure_drms_environment($cfg);
  validate_drms_runtime( $cfg->{show_info} );
  my ( $target_epoch, $previous_epoch, $next_epoch ) =
    ( tstr2ts($slot), tstr2ts($previous), tstr2ts($next) );
  die "Invalid interpolation timestamp\n"
    unless defined $target_epoch && defined $previous_epoch && defined $next_epoch;
  die "Interpolation target must be strictly between its sources\n"
    if $previous_epoch >= $target_epoch || $target_epoch >= $next_epoch;

  my @keys = map { ( sprintf( 'A_%.3d_X0', $_ ), sprintf( 'A_%.3d_Y0', $_ ) ) } @wl;
  my @source;
  for my $time ( $previous, $next ) {
    my @lines = show_info_lines(
      $cfg->{show_info}, '-q',
      'key=' . join( q{,}, @keys ),
      sprintf( '%s[%s]', $cfg->{mpt_series}, $time )
    );
    chomp @lines;
    @lines = grep { /\S/ } @lines;
    die "Expected one pointing record at $time, found " . scalar(@lines) . "\n" unless @lines == 1;
    my @values = split /\s+/, $lines[0];
    die "Incomplete pointing record at $time\n" unless @values == @keys;
    for my $i ( 0 .. $#values ) {
      die "Cannot interpolate $keys[$i]: invalid value '$values[$i]' at $time\n"
        if !looks_like_number( $values[$i] )
        || $values[$i] =~ /(?:nan|inf)/i
        || ( defined $cfg->{nan_sentinel} && $values[$i] == $cfg->{nan_sentinel} );
    }
    push @source, \@values;
  }
  my $fraction = ( $target_epoch - $previous_epoch ) / ( $next_epoch - $previous_epoch );
  for my $i ( 0 .. $#wl ) {
    push @kwd,
      kwd_lines(
      $wl[$i],
      $source[0][ $i * 2 ] + $fraction * ( $source[1][ $i * 2 ] - $source[0][ $i * 2 ] ),
      $source[0][ $i * 2 + 1 ] + $fraction * ( $source[1][ $i * 2 + 1 ] - $source[0][ $i * 2 + 1 ] )
      );
  }
  warn "INTERPOLATION OK slot=$slot previous=$previous next=$next fraction=$fraction\n";
}
else {
  my @missing = grep { !-s limb_path( $inpdir, $yr, $mo, $da, $hr, $_ ) } @wl;
  die "Missing or empty limb files: " . join( ', ', map { "${_}A" } @missing ) . "\n" if @missing;

  my @failed;
  for my $w (@wl) {
    my $path = limb_path( $inpdir, $yr, $mo, $da, $hr, $w );
    my ( $xavg, $yavg );
    my $reduced = eval {
      my ( $x0, $y0 ) = rcols $path, 0, 1;
      ( undef, undef, $xavg, $yavg ) = reduce_limb_points( $x0, $y0, $w, $cfg );
      1;
    };
    my $failure = $@;
    $failure = "non-finite reducer result\n"
      if $reduced && ( "$xavg" =~ /(?:nan|inf)/i || "$yavg" =~ /(?:nan|inf)/i );
    if ( !$reduced || $failure ) {
      chomp $failure;
      warn "REDUCER FAIL slot=$slot wavelength=${w}A reason=$failure\n";
      push @failed, "${w}A: $failure";
      next;
    }
    push @kwd, kwd_lines( $w, $xavg, $yavg );
  }

  die "Reducer failed for $slot: " . join( '; ', @failed ) . "\n" if @failed;
}
make_path( $outdir, { chmod => oct('755') } ) unless -d $outdir;
my $temporary = "$outnam.$$";
my $written   = eval {
  open my $fh, '>', $temporary or die "Can't open '$temporary' for writing: $!\n";
  print {$fh} @kwd;
  close $fh or die "Can't close '$temporary': $!\n";
  rename $temporary, $outnam or die "Can't replace '$outnam': $!\n";
  1;
};
if ( !$written ) {
  my $error = $@ || "Unknown write failure\n";
  unlink $temporary if -e $temporary;
  die $error;
}
my $partial_note =
  @requested_wavelengths ? q{ partial=} . join( q{,}, map { $_ . q{A} } @wl ) : q{};
warn q{MASTERPOINT OK slot=} . $slot . q{ output=} . $outnam . $partial_note . "\n";