-- Ada95 binding to the netcdf library
-- Copyright (c) 2002-2003 Alexandru Dan Corlan (http://dan.corlan.net)
--
-- This program is free software; you can redistribute it and/or modify
-- it under the terms of the GNU Lesser General Public License as published by
-- the Free Software Foundation; either version 2 of the License, or
-- (at your option) any later version.
--
-- As I understand the licence you must distribute this source file,
-- unchanged, with any binary distribution which uses it, or provide
-- a reference to it as specified in the LGPL licence.
--
-- This program is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-- GNU Lesser General Public License for more details.
--
-- You could have received a copy of the GNU General Public License
-- along with this program; if not, write to the Free Software
-- Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
-- or visit http://www.gnu.org
--
-- HISTORY:
-- jul 2003 v0.1 -- initial release, basic thin binding for reading and inquiring
-- dec 2003 v0.2 -- added `varead' for single call extraction of vectors and
--                  matrices from

with Interfaces.C; use Interfaces.C;
with Ada.Text_Io; use Ada.Text_Io;

package Netcdf is

   type Dataset is new Integer;
   type Dim_Id is new Integer;
   type Dim_Set is array(Integer range <>) of Dim_Id;
   type Var_Id is new Integer;
   type Index is array(Integer range <>) of Natural;
   type Attribute_Nr is new Integer;

   type Nc_Type is (Nc_Nat, Nc_Byte, Nc_Char, Nc_Short, Nc_Int, Nc_Float, Nc_Double);
   for Nc_Type use ( Nc_Nat => 0, Nc_Byte => 1, Nc_Char => 2, Nc_Short => 3,
                     Nc_Int => 4, Nc_Float => 5, Nc_Double => 6);
   for Nc_Type'Size use 32;

   type Float_Vector is array(Integer range <>) of Float;
   type Float_Matrix is array(Integer range <>, Integer range <>) of Float;

   function Inq_Libvers return String;

    procedure Create(Path: String;
                     Clobber: Boolean:= True;
                     Share: Boolean:= False;
                     Ncid: out Dataset);

    procedure Open(Path: String;
                   Write: Boolean:= False;
                   Share: Boolean:= False;
                   Ncid: out Dataset);

--    procedure Redef(Ncid: Dataset);

--    procedure Enddef(Ncid: Dataset);

    procedure Close(Ncid: Dataset);

--    procedure Inq(Ncid: Dataset;
--                  Ndims: out Integer;
--                  Nvars: out Integer;
--                  Natts: out Integer;
--                  Unlimdimid: out Dim_Id);

    function Inq_Ndims(Ncid: Dataset) return Integer;
    function Inq_Nvars(Ncid: Dataset) return Integer;
    function Inq_Natts(Ncid: Dataset) return Integer;
    function Inq_Unlimdim(Ncid: Dataset) return Dim_Id;

--    procedure Sync(Ncid: Dataset);

--    procedure Def_Dim(Ncid: Dataset;
--                      Name: String;
--                      Len: Natural;
--                      Id: out Dim_Id);

--    function Dimid(Ncid: Dataset; Name: String) return Dim_Id;

    function Dimname(Ncid: Dataset; Id: Dim_Id) return String;
    function Dimlen(Ncid: Dataset; Id: Dim_Id) return Natural;

--    procedure Rename_Dim(Ncid: Dataset; Id: Dim_Id; New_Name: String);

--    procedure Def_Var(Ncid: Dataset;
--                      Name: String;
--                      Xtype: Nc_Type;
--                      Dims: Dim_Set;
--                      Id: out Var_Id);

    function Varid(Ncid: Dataset; Name: String) return Var_Id;
    function Varname(Ncid: Dataset; Var: Var_Id) return String;
    function Vartype(Ncid: Dataset; Var: Var_Id) return Nc_Type;
    function Varndims(Ncid: Dataset; Var: Var_Id) return Natural;
    function Vardims(Ncid: Dataset; Var: Var_Id) return Dim_Set;
    function Varnatts(Ncid: Dataset; Var: Var_Id) return Natural;

--    function Rename_Var(Ncid: Dataset; Var: Var_Id; New_Name: String);

--    procedure Put(Ncid: Dataset; Var: Var_Id; Ix: Index; C: Character);
--    procedure Put(Ncid: Dataset; Var: Var_Id; Ix: Index; C: Integer);
--    procedure Put(Ncid: Dataset; Var: Var_Id; Ix: Index; C: Long_Long_Integer);
--    procedure Put(Ncid: Dataset; Var: Var_Id; Ix: Index; C: Float);
--    procedure Put(Ncid: Dataset; Var: Var_Id; Ix: Index; C: Double);

--    function Get(Ncid: Dataset; Var: Var_Id; Ix: Index) return Character;
      function Get(Ncid: Dataset; Var: Var_Id; Ix: Index) return Integer;
--    function Get(Ncid: Dataset; Var: Var_Id; Ix: Index) return Long_Long_Integer;
      function Get(Ncid: Dataset; Var: Var_Id; Ix: Index) return Float;
--    function Get(Ncid: Dataset; Var: Var_Id; Ix: Index) return Double;

      function Varead(File_Name: String; Var_Name: String) return Float_Matrix;
      function Varead(File_Name: String; Var_Name: String) return Float_Vector;
      Wrong_Aggregate_Dimensionality: exception;

end Netcdf;


with Interfaces.C; use Interfaces.C;
with Interfaces.C.Strings; use Interfaces.C.Strings;
with System; use System;

package body Netcdf is

   type Modbits is mod 65536;

   type Dataset_Ptr is access all Dataset;

   -- the following constants are taken from netcdf.h, V2.74, 2000/10/13
   -- distributed with Debian Linux 3.0 (Woody) for the i386 architecture
   -- they are unlikely to change in future versions, but it is always worth a check

   Nc_Nowrite:   constant Modbits:= 16#0000#;
   Nc_Write:     constant Modbits:= 16#0001#;
   Nc_Clobber:   constant Modbits:= 16#0000#;
   Nc_Noclobber: constant Modbits:= 16#0004#;
   Nc_Fill:      constant Modbits:= 16#0000#;
   Nc_Nofill:    constant Modbits:= 16#0100#;
   Nc_Lock:      constant Modbits:= 16#0400#;
   Nc_Share:     constant Modbits:= 16#0800#;

   Nc_Max_Name:  constant Integer:= 128;

   pragma Linker_Options("-lnetcdf");

   function Nc_Inq_Libvers return Chars_Ptr;
   pragma Import(C,Nc_Inq_Libvers,"nc_inq_libvers");

   function Inq_Libvers return String is

   begin
      return Value(Nc_Inq_Libvers);
   end Inq_Libvers;

   procedure Nc_Create(Path: System.Address; Cmode: Modbits; Ncid: System.Address);
   pragma Import(C,Nc_Create,"nc_create");

   procedure Create(Path: String;
                    Clobber: Boolean:= True;
                    Share: Boolean:= False;
                    Ncid: out Dataset) is

      Tonc: aliased Dataset;
      Flgdata: Modbits:= 0;
      Path_0: String:= Path & Ascii.Nul;

   begin
      if Clobber then
         Flgdata := Nc_Clobber;
      end if;
      if Share then
         Flgdata := Flgdata or Nc_Share;
      end if;
      Nc_Create(Path_0'Address, Flgdata, Tonc'Address);
      Ncid := Tonc;
   end Create;

   procedure Nc_Open(Path: System.Address; Cmode: Modbits; Ncid: System.Address);
   pragma Import(C,Nc_Open,"nc_open");

   procedure Open(Path: String;
                  Write: Boolean:= False;
                  Share: Boolean:= False;
                  Ncid: out Dataset) is

      Tonc: aliased Dataset;
      Flgdata: Modbits:= 0;
      Path_0: String:= Path & Ascii.Nul;

   begin
      if Write then
         Flgdata := Nc_Write;
      end if;
      if Share then
         Flgdata := Flgdata or Nc_Share;
      end if;
      Nc_Open(Path_0'Address, Flgdata, Tonc'Address);
      Ncid := Tonc;
   end Open;

   procedure Nc_Close(Ncid: Dataset);
   pragma Import(C,Nc_Close,"nc_close");

   procedure Close(Ncid: Dataset) is

   begin
      Nc_Close(Ncid);
   end Close;

   procedure Nc_Inq_Nvars(Ncid: Dataset; Nvars: System.Address);
   pragma Import(C,Nc_Inq_Nvars,"nc_inq_nvars");

   function Inq_Nvars(Ncid: Dataset) return Integer is

      R: Integer;

   begin
      Nc_Inq_Nvars(Ncid,R'Address);
      return R;
   end Inq_Nvars;

   procedure Nc_Inq_Ndims(Ncid: Dataset; Ndims: System.Address);
   pragma Import(C,Nc_Inq_Ndims,"nc_inq_ndims");

   function Inq_Ndims(Ncid: Dataset) return Integer is

      R: Integer;

   begin
      Nc_Inq_Ndims(Ncid,R'Address);
      return R;
   end Inq_Ndims;

   procedure Nc_Inq_Natts(Ncid: Dataset; Natts: System.Address);
   pragma Import(C,Nc_Inq_Natts,"nc_inq_natts");

   function Inq_Natts(Ncid: Dataset) return Integer is

      R: Integer;

   begin
      Nc_Inq_Natts(Ncid,R'Address);
      return R;
   end Inq_Natts;

   procedure Nc_Inq_Unlimdim(Ncid: Dataset; Dim: System.Address);
   pragma Import(C,Nc_Inq_Unlimdim,"nc_inq_unlimdim");

   function Inq_Unlimdim(Ncid: Dataset) return Dim_Id is

      R: Dim_Id;

   begin
      Nc_Inq_Unlimdim(Ncid,R'Address);
      return R;
   end Inq_Unlimdim;

   procedure Nc_Inq_Dimname(Ncid: Dataset; Dimid: Dim_Id; Name: System.Address);
   pragma Import(C,Nc_Inq_Dimname,"nc_inq_dimname");

   function Dimname(Ncid: Dataset; Id: Dim_Id) return String is

      Nm: Char_Array(1..Size_T(Nc_Max_Name));

   begin
      Nc_Inq_Dimname(Ncid,Id,Nm'Address);
      return To_Ada(Nm, Trim_Nul => True);
   end Dimname;

   procedure Nc_Inq_Dimlen(Ncid: Dataset; Dim: Dim_Id; Len: System.Address);
   pragma Import(C,Nc_Inq_Dimlen,"nc_inq_dimlen");

   function Dimlen(Ncid: Dataset; Id: Dim_Id) return Natural is

      R: Natural;

   begin
      Nc_Inq_Dimlen(Ncid,Id,R'Address);
      return R;
   end Dimlen;

   procedure Nc_Inq_Varid(Ncid: Dataset; Name: System.Address; Id: System.Address);
   pragma Import(C,Nc_Inq_Varid,"nc_inq_varid");

   function Varid(Ncid: Dataset; Name: String) return Var_Id is

      R: Var_Id;
      Varname: String:= Name & Ascii.Nul;

   begin
      Nc_Inq_Varid(Ncid,Varname'Address,R'Address);
      return R;
   end Varid;

   procedure Nc_Inq_Varname(Ncid: Dataset; Varid: Var_Id; Name: System.Address);
   pragma Import(C,Nc_Inq_Varname,"nc_inq_varname");

   function Varname(Ncid: Dataset; Var: Var_Id) return String is

      Nm: Char_Array(1..Size_T(Nc_Max_Name));

   begin
      Nc_Inq_Varname(Ncid,Var,Nm'Address);
      return To_Ada(Nm, Trim_Nul => True);
   end Varname;

   procedure Nc_Inq_Varndims(Ncid: Dataset; V: Var_Id; Len: System.Address);
   pragma Import(C,Nc_Inq_Varndims,"nc_inq_varndims");

   function Varndims(Ncid: Dataset; Var: Var_Id) return Natural is

      R: Natural;

   begin
      Nc_Inq_Varndims(Ncid,Var,R'Address);
      return R;
   end Varndims;

   procedure Nc_Inq_Vardimid(Ncid: Dataset; V: Var_Id; Dimids: System.Address);
   pragma Import(C,Nc_Inq_Vardimid,"nc_inq_vardimid");

   function Vardims(Ncid: Dataset; Var: Var_Id) return Dim_Set is

      R: Dim_Set(1..Varndims(Ncid,Var));

   begin
      Nc_Inq_Vardimid(Ncid,Var,R'Address);
      return R;
   end Vardims;

   procedure Nc_Inq_Vartype(Ncid: Dataset; V: Var_Id; Nct: System.Address);
   pragma Import(C,Nc_Inq_Vartype,"nc_inq_vartype");

   function Vartype(Ncid: Dataset; Var: Var_Id) return Nc_Type is

      R: Nc_Type;

   begin
      Nc_Inq_Vartype(Ncid,Var,R'Address);
      return R;
   end Vartype;

   procedure Nc_Inq_Varnatts(Ncid: Dataset; V: Var_Id; Len: System.Address);
   pragma Import(C,Nc_Inq_Varnatts,"nc_inq_varnatts");

   function Varnatts(Ncid: Dataset; Var: Var_Id) return Natural is

      R: Natural;

   begin
      Nc_Inq_Varnatts(Ncid,Var,R'Address);
      return R;
   end Varnatts;

   procedure Nc_Get_Var1_Int(Ncid: Dataset; Var: Var_Id;
                                           Indx: System.Address; Into: System.Address);
   pragma Import(C,Nc_Get_Var1_Int,"nc_get_var1_int");

   function Get(Ncid: Dataset; Var: Var_Id; Ix: Index) return Integer is

      R: Integer;

   begin
      Nc_Get_Var1_Int(Ncid,Var,Ix'Address,R'Address);
      return R;
   end Get;


   procedure Nc_Get_Var1_Float(Ncid: Dataset; Var: Var_Id;
                                           Indx: System.Address; Into: System.Address);
   pragma Import(C,Nc_Get_Var1_Float,"nc_get_var1_float");

   function Get(Ncid: Dataset; Var: Var_Id; Ix: Index) return Float is

      R: Float;

   begin
      Nc_Get_Var1_Float(Ncid,Var,Ix'Address,R'Address);
      return R;
   end Get;

   function Varead(File_Name: String; Var_Name: String) return Float_Matrix is

      Vid: Var_Id;
      Id: Dataset;

   begin
      Open(File_Name,Write => False, Ncid => Id);
      Vid := Varid(Id,Var_Name);

      declare
         Ds: Dim_Set:= Vardims(Id,Vid);
         Xsz, Ysz: Integer;
      begin

         if Ds'Length/=2 then
            Close(Id);
            raise Wrong_Aggregate_Dimensionality;
         end if;
         Xsz := Dimlen(Id,Ds(1));
         Ysz := Dimlen(Id,Ds(2));

         declare
            R: Float_Matrix(1..Xsz,1..Ysz);
            Ff: Float:= -221.0;
         begin
            for X in 0..(Xsz-1) loop
               for Y in 0..(Ysz-1) loop
                  Ff := Get(Id,Vid,(X, Y));
                  R(X+1,Y+1) := Ff;
               end loop;
            end loop;
            Close(Id);
            return R;
         end;
      end;
   end Varead;

   function Varead(File_Name: String; Var_Name: String) return Float_Vector is

      Vid: Var_Id;
      Id: Dataset;

   begin
      Open(File_Name,Write => False, Ncid => Id);
      Vid := Varid(Id,Var_Name);

      declare
         Ds: Dim_Set:= Vardims(Id,Vid);
         Xsz: Integer;
      begin
         if Ds'Length/=1 then
            Close(Id);
            raise Wrong_Aggregate_Dimensionality;
         end if;
         Xsz := Dimlen(Id,Ds(1));

         declare
            R: Float_Vector(1..Xsz);
         begin
            for X in 0..(Xsz-1) loop
               R(X+1) := Get(Id,Vid,(1 => X));
            end loop;
            Close(Id);
            return R;
         end;
      end;
   end Varead;

end Netcdf;

with Netcdf; use Netcdf;
with Text_Io; use Text_Io;

procedure Try_Netcdf is

   Id: Dataset;
   Ndims: Integer;
   Nvars: Integer;

begin
   Put_Line(Inq_Libvers);
   Open("tilt02.cdf",Write => False, Ncid => Id);
   Put_Line("Opened tilt02.cdf");
   Nvars := Inq_Nvars(Id);
   Put("Nr of variables  = " & Integer'Image(Nvars) & ": ");
   for I in 1..Nvars loop
      Put(Varname(Id,Var_Id(I-1))& " ");
      Put("(" & Nc_Type'Image(Vartype(Id,Var_Id(I-1))) & ", ");
      declare
         Ds: Dim_Set:= Vardims(Id,Var_Id(I-1));
      begin
         for I in Ds'Range loop
            Put(Dimname(Id,Ds(I)) & "=" & Natural'Image(Dimlen(Id,Ds(I))));
            if I<Ds'Last then
               Put(")");
            end if;
         end loop;
      end;
      Put(")");
   end loop;
   Put_Line("");
   Ndims := Inq_Ndims(Id);
   Put("Nr of dimensions = " & Integer'Image(Ndims) & ": ");
   for I in 0..(Ndims-1) loop
      Put(Dimname(Id,Dim_Id(I)) & " ");
   end loop;
   Put_Line("");
   Put_Line("First 10 data elements, as integers:");
   for I in 1..10 loop
      for V in 1..Nvars loop
         Put(Integer'Image(Get(Id,Var_Id(V-1),(1 => I))) & " ");
      end loop;
      Put_Line("");
   end loop;
   Close(Id);
end Try_Netcdf;