#! /usr/bin/env perl
#
#    This is an example program which reads some 4D pressure and
#    temperatures. The companion program pres_temp_4D_wr.pl shows how
#    to write the netCDF data file read by this program.
#
# Heiko Klein <Heiko.Klein@met.no> 2007-02-06
#

use strict;
use warnings;

use PDL::Lite;
use PDL::NetCDF;
use Fcntl;

my $NLAT = 6;
my $NLON = 12;
my $SAMPLE_PRESSURE = 900.0;
my $SAMPLE_TEMP = 9.0;
my $START_LAT = 25.0;
my $START_LON = -125.0;
my $NREC = 2;
my $NLEVEL = 2;

my $latitude = (PDL::Basic::sequence($NLAT) * 5 + $START_LAT)->float;
my $longitude = (PDL::Basic::sequence($NLON) * 5 + $START_LON)->float;
my $pressure = PDL::Core::zeroes($NLON, $NLAT, $NLEVEL)->float;
my $temp = PDL::Core::zeroes($NLON, $NLAT, $NLEVEL)->float;
for (my $l = 0; $l < $NLEVEL; $l++) {
    for (my $i = 0; $i < $NLAT; $i++) {
	for (my $j = 0; $j < $NLON; $j++) {
	    $pressure->mslice($j, $i, $l) .= $SAMPLE_PRESSURE + $i*$NLON + $j + $l*$NLON*$NLAT; # no change in t here
	    $temp->mslice($j, $i, $l) .= $SAMPLE_TEMP + $i*$NLON + $j + $l*$NLAT*$NLON; # no change in t here
	}
    }
}

eval {
    my $ncfile = new PDL::NetCDF("pres_temp_4D.nc", {REVERSE_DIMS => 1,
						     MODE => O_RDONLY});

    # checking data
    if (not equalPdls($latitude, $ncfile->get('latitude'))) {
	die "different latitude values";
    }
    if (not equalPdls($longitude, $ncfile->get('longitude'))) {
	die "different longitude values";
    }
    for (my $t = 0; $t < $NREC; $t++) {
	if (not equalPdls($pressure, $ncfile->get('pressure', [0,0,0,$t], [$NLON, $NLAT, $NLEVEL, 1]))) {
	    die "different pressure values, time $t";
	}
	if (not equalPdls($temp, $ncfile->get('temperature'))) {
	    die "different temperature values, time $t";
	}
    }

    # checking attributes
    if ('degrees_north' ne $ncfile->getatt('units', 'latitude')) {
	die "wrong units for latitude";
    }
    if ('degrees_east' ne $ncfile->getatt('units', 'longitude')) {
	die "wrong units for longtude";
    }
    if ('hPa' ne $ncfile->getatt('units', 'pressure')) {
	die "wrong units for pressure";
    }
    if ('celcius' ne $ncfile->getatt('units', 'temperature')) {
	die "wrong units for temperature";
    }


$ncfile->close;
};
if ($@) {
    print STDERR "$@\n";
} else {
    print "SUCCESS reading file\n";
}


sub equalPdls {
    my ($pdlA, $pdlB) = @_;
    return $pdlA->where($pdlA != $pdlB)->isempty;
}