function values = getcdf(file)

%
% Modified 6/92 to use the c version of mexcdf. J. Allen
%
%--------------------------------------------------------------------
%     Copyright (C) J. V. Mansbridge, CSIRO, 
% CHANGE   1.5 92/03/29
%
%  function values = getcdf(file)
%
% DESCRIPTION:
%  getcdf is an interactive function that gets a NetCDF hyperslab in
%  a user-friendly manner.  The user is asked questions to determine
%  information necessary for the data retrieval.  If the file is in
%  compressed form then the user will be given an option for
%  automatic uncompression.
%
% INPUT:
%  file is the name of a netCDF file with or without the .cdf extent.
%
% OUTPUT:
%  values is a scalar, vector or matrix of data that is read in from
%  the NetCDF file
%
% EXAMPLE:
%  velocity = getcdf('fred')
%
% CALLER:   general purpose
% CALLEE:   fill_attnames.m, fill_varnames.m, check_cdf.m, check_string.m,
%           return_val.m, mexcdf.mex
%
% AUTHOR:   J. V. Mansbridge, CSIRO
%---------------------------------------------------------------------

% @(#)getcdf.m   1.5   92/03/29
% 
%--------------------------------------------------------------------

% Set some constants.

blank = abs(' ');

% Check that there are the correct number of arguments.

if nargin ~= 1
   help getcdf
   return
end

% Check that the file is accessible and optionally uncompress it.

nf=length(file);
if(file(nf-3:nf)~='.cdf'),
  cdf = [ file '.cdf' ];
else
  cdf = file;
end
err = check_cdf(cdf);
if err ~= 0
    return;
end

% Open the netcdf file.

[cdfid, rcode] = mexcdf('OPEN', cdf, 'NC_NOWRITE');

% Suppress all error messages from netCDF 

[rcode] = mexcdf('setopts', 0);

% Collect information about the cdf file.

[ndims, nvars, ngatts, recdim, rcode] =  mexcdf('ncinquire', cdfid);

varstring = fill_varnames(cdfid, nvars);

% Prompt the user for the name of the variable containing the hyperslab.

k = -1;
while k <1 | k > nvars
   disp(' ')
   s = [ '----- Choose the variable -----'];
   disp(s)
   disp(' ')
   for i = 0:3:nvars-1
      stri = int2str(i+1);
      if length(stri) == 1
         stri = [ ' ' stri];
      end
      [varnam, vartyp, nvdims, vdims, nvatts, rcode] = ...
      mexcdf('ncvarinq', cdfid, i);
      s = [ '  ' stri ') ' varnam ];
      addit = 26 - length(s);
      for j =1:addit
         s = [ s ' '];
      end

      if i < nvars-1
         stri = int2str(i+2);
         if length(stri) == 1
            stri = [ ' ' stri];
         end
         [varnam, vartyp, nvdims, vdims, nvatts, rcode] = ...
         mexcdf('ncvarinq', cdfid, i+1);
         s = [ s '  ' stri ') ' varnam ];
         addit = 52 - length(s);
         for j =1:addit
            s = [ s ' '];
         end
      end 

      if i < nvars - 2
         stri = int2str(i+3);
         if length(stri) == 1
            stri = [ ' ' stri];
         end
         [varnam, vartyp, nvdims, vdims, nvatts, rcode] = ...
         mexcdf('ncvarinq', cdfid, i+2);
         s = [ s '  ' stri ') ' varnam ];
      end 
      disp(s)
   end
   disp(' ')
   s = [ 'Select a menu number: '];
%   k = input(s);
   k = return_val(s, 0);
end

% try to get information about the variable

varid = k-1;
[varnam, vartyp, nvdims, vdims, nvatts, rcode] = ...
        mexcdf('ncvarinq', cdfid, varid);
attstring = fill_attnames(cdfid, varid, nvatts);

% Print out information about the dimensions of the variable and calculate
% some things needed for the subsequent questions.

disp(' ')
disp([varnam ' has ' int2str(nvdims) ' dimensions'])
nadim = '';
lcount = 1;
message = '';
mcount = 1;
for i = 1:nvdims
    dimid = vdims(i);
    [name, size, rcode] = mexcdf('ncdiminq', cdfid, dimid);

% If the variable is one dimensional then check whether it has length 1.
% If it does then find and return its value and exit.

   if size == 1 & nvdims == 1
      [values, rcode] = mexcdf('ncvarget1', cdfid, varid, [0]);
      return
   end
   llow(i) = lcount;
   lcount = lcount + length(name);
   lup(i) = lcount - 1;
   nadim = [nadim name];
   ledim(i) = size-1;

% Test that the dimension name is also a variable name.  If it is then
% store information about its initial and final values in the string s.

   rhid = check_string(name, varstring, nvars)-1;
   if rhid >= 0
      [namejunk, dvartyp, dnvdims, vdimsjunk, nvattsjunk, rcode] = ...
                mexcdf('ncvarinq', cdfid, rhid);
      if size <= 6
         [temp, rcode] = mexcdf('ncvarget', cdfid, rhid, [0], [size]);
         s = ' : Elements';
         for j = 1:size
            s = [ s ' ' num2str(temp(j)) ];
         end
      else
         [temp1, rcode] = mexcdf('ncvarget', cdfid, rhid, [0], [3]);
         s = ' : Elements';
         for j = 1:3
            s = [ s ' ' num2str(temp1(j)) ];
         end
         [temp2, rcode] = mexcdf('ncvarget', cdfid, rhid, [size-3], [3]);
         s= [ s ' ...' ];
         for j = 1:3
            s = [ s ' ' num2str(temp2(j)) ];
         end
      end
   else
      s = [ ' '];
   end
   s = [ '   ' int2str(i) ')  ' name ' : Length ' int2str(size) s ];
   mlow(i) = mcount;
   mcount = mcount + length(s);
   mup(i) = mcount - 1;
   message = [ message s ];
   disp(s)
end

% Ask the user for the number of dimensions and check the answer.

if nvdims == 1
   s = 'What sort of hyperslab do you want?';
   s1 = 'a scalar';
   s2 = 'a vector';
   k = -1;
   while any(k == [1 2]) == 0
      k = menu(s, s1, s2);
      if any(k == [1 2]) == 0
         disp(' ')
         disp('You have asked for a non-existent option - try again')
      end
   end
else
   s = 'What sort of hyperslab do you want?';
   s1 = 'a scalar';
   s2 = 'a vector';
   s3 = 'a matrix';
   k = -1;
   while any(k == [1 2 3]) == 0
      k = menu(s, s1, s2, s3);
      if any(k == [1 2 3]) == 0
         disp(' ')
         disp('You have asked for a non-existent option - try again')
      end
   end
end
num_edge = k - 1;

% ask for the corners and edges according to the number of dimensions.

if num_edge == 0

% ask for the corners in order to retrieve a scalar.

   for i = 1:nvdims
      name = [nadim(llow(i):lup(i))];
      corner(i) = -1;
      while corner(i) < 0 | corner(i) > ledim(i)
         s = [ message(mlow(i)+2:mup(i)) ];
         disp(' ')
         disp(s)
         s=[ '    ' name ' : Index (between 0 and ' ...
            int2str(ledim(i)) ')  ' ];
%         corner(i) = input(s);
         corner(i) = return_val(s, 0);
      end
   end

% Get the scalar.

   [values, rcode] = mexcdf('ncvarget1', cdfid, varid, corner);

else

% identify which dimensions are on the edge of the hyperslab in 2 parts.

% first get the necessary information from the user.

   y_count = 0;
   n_count = 0;
   n_lim = nvdims - num_edge;
   for i = 1:nvdims
      choice(i) = ' ';
   end
   i = 1;
   while y_count < num_edge & n_count < n_lim
      name = [nadim(llow(i):lup(i))];
      st = [ 'Should the hyperslab vary with ... ' name '? (y/n) [y] '];
      reply = input(st, 's');
      if isempty(reply)
         choice(i) = 'y';
         y_count = y_count + 1;
      elseif reply == 'y' | reply == 'Y'
         choice(i) = 'y';
         y_count = y_count + 1;
      elseif reply == 'n' | reply == 'N'
         choice(i) = 'n';
         n_count = n_count + 1;
      end
      if i == nvdims
         i = 1;
      else
         i = i + 1;
      end
   end
   disp(' ')

% second fill in the remainder of the choice vector

   if y_count == num_edge
      for i = 1:nvdims
         if choice(i) ~= 'y';
            choice(i) = 'n';
         end
      end
   else
      for i = 1:nvdims
         if choice(i) ~= 'n';
            choice(i) = 'y';
         end
      end
   end



% ask for the index at a point or the corners and edges in order to retrieve
% a hyperslab.

   for i = 1:nvdims

      if choice(i) == 'y' | choice(i) == 'Y'

% first get the starting point

         name = [nadim(llow(i):lup(i))];
         corner(i) = -1;
         while corner(i) < 0 | corner(i) > ledim(i)
            s = [ message(mlow(i)+2:mup(i)) ];
            disp(' ')
            disp(s)
            s=[ '    ' name ' : Starting index (between 0 and ' ...
            int2str(ledim(i)) ')  (cr for all indices)  ' ];
            clear xtemp;
            xtemp = input(s);
            if isempty(xtemp)
               corner(i) = 0;
               edge(i) = ledim(i)+1;
               notdone = 0;
            else
               corner(i) = xtemp;
               notdone = 1;
            end
         end

% next, get the finishing point if this is required.

         if notdone
            end_point = -1;
            ste = [];
            for ii = 1:length(name)
               ste = [ ste ' ' ];
            end
            while end_point < corner(i) | end_point > ledim(i)
               s=[ ste '      finishing index (between ' ...
                  int2str(corner(i)) ' and ' int2str(ledim(i)) ')  '];
               end_point = return_val(s, end_point);
%               clear xtemp;
%               xtemp = input(s);
%               if isempty(xtemp) == 0
%                  end_point = xtemp;
%               end
            end

% Calculate the edge length

            edge(i) = end_point - corner(i) + 1;
         end
      else

% get the index for this dimension

         name = [nadim(llow(i):lup(i))];
         corner(i) = -1;
         while corner(i) < 0 | corner(i) > ledim(i)
            s = [ message(mlow(i)+2:mup(i)) ];
            disp(' ')
            disp(s)
            s=[ '    ' name ' : Index (between 0 and ' ...
               int2str(ledim(i)) ')  ' ];
%            clear xtemp;
%            xtemp = input(s);
%            if isempty(xtemp) == 0
%               corner(i) = xtemp;
%            end
           corner(i) = return_val(s, 0);
         end

% set the edge length to 1

         edge(i) = 1;
      end
   end

   if num_edge == 1

% get the row vector

       lenstr = prod(edge);
       [values, rcode] = mexcdf('ncvarget', cdfid, varid, corner, edge);
% Do not transpose      values = values.';

   else

% If we have wish to return a matrix then some manipulations must be done.

% First identify the dimensions and their order in the matrix.  Note that the
% order of the dimensions is from slowest changing to fastest as we
% are making C calls (not FORTRAN).

      count = 0;
      for i = 1:nvdims
         if choice(i) == 'y' | choice(i) == 'Y'
            if count == 0
               name1 = [nadim(llow(i):lup(i))];
               length1 = edge(i);
               count = 1;
            else
               name2 = [nadim(llow(i):lup(i))];
               length2 = edge(i);
               count = 2;
            end
         end
      end

% ask about the order of the indicies.

      order = 0;
      s = 'In which order do you want the indicies?';
      s1 = [ varnam '(' name2 ',' name1 ')' ];
      s2 = [ varnam '(' name1 ',' name2 ')' ];
      while any(order == [1 2]) == 0
         order = menu(s, s1, s2);
      end

% create the appropriate 2-d matrix according to the value of order.  If
% order == 2 then the 2-d matrix must be transposed.

       lenstr = prod(edge);
       [values, rcode] = mexcdf('ncvarget', cdfid, varid, corner, edge);
	length1
	length2
%      values = reshape(values, length1, length2);
       values = reshape(values, length2, length1);
      if order == 2
         values = values.';
      end

   end
end

% check for missing values.  Note that a
% missing value is taken to be one less than valid_min, greater than
% valid_max or 'close to' _FillValue.
% Note 1: valid_min and valid_max may be specified by the attribute
%   valid_range and if valid_range exists than the existence of
%   valid_min and valid_max is not checked.
% Note 2: a missing value must be OUTSIDE the valid range to be
%   recognised.
% Note 3: a range does not make sense for character arrays.
% Note 4: By 'close to' _FillValue I mean that an integer or character
%   must equal _FillValue and a real must be in the range
%   0.99999*_FillValue tp 1.00001*_FillValue.  This allows real*8 
%   rounding errors in moving the data from the netcdf file to matlab;
%   these errors do occur although I don't know why given that matlab
%   works in double precision.

% The indices of the data points containing missing value indicators
% will be stored separately in index_miss_low, index_miss_up and
% index__FillValue

index_miss_low = [];
index_miss_up = [];
index__FillValue = [];

% First find the indices of the data points that are outside the valid
% range.

pos_vr = check_string('valid_range', attstring, nvatts);
if pos_vr > 0
   [attype, attlen, rcode] = mexcdf('ncattinq', cdfid, varid, 'valid_range');
   [ miss, rcode] = mexcdf('ncattget', cdfid, varid, 'valid_range');
   miss_low = miss(1);
   miss_up = miss(2);
   index_miss_low = find ( values < miss_low );
   index_miss_up = find ( values > miss_up );
else
   pos_min = check_string('valid_min', attstring, nvatts);
   if pos_min > 0
      [attype, attlen, rcode] = mexcdf('ncattinq', cdfid, varid, 'valid_min');
      [miss_low, rcode] = mexcdf('ncattget', cdfid, varid, 'valid_min');
      index_miss_low = find ( values < miss_low );
   end

   pos_max = check_string('valid_max', attstring, nvatts);
   if pos_max > 0
      [attype, attlen, rcode] = mexcdf('ncattinq', cdfid, varid, 'valid_max');
      [miss_up, rcode] = mexcdf('ncattget', cdfid, varid, 'valid_max');
      index_miss_up = find ( values > miss_up );
   end
end

% Now find the indices of the data points that are 'close to'
% _FillValue.  Note that 'close to' is different according to the
% data type.

pos_missv = check_string('_FillValue', attstring, nvatts);
if pos_missv > 0
   [attype, attlen, rcode] = mexcdf('ncattinq', cdfid, varid, '_FillValue');
   [miss_val, rcode] = mexcdf('ncattget', cdfid, varid, '_FillValue');
   if attype == 1 | attype == 2
      index__FillValue = find ( values == miss_val );
   elseif attype == 3 | attype == 4
      need_index_m = 1;
      if pos_vr > 0 | pos_min > 0
         if miss_val < miss_low
            need_index_m = 0;
         end
      end
      if pos_vr > 0 | pos_max > 0
         if miss_val > miss_up
            need_index_m = 0;
         end
      end
      if need_index_m
         index__FillValue = find ( values == miss_val );
      end
   elseif attype == 5 | attype == 6
      need_index_m = 1;
      if miss_val < 0
         miss_val_low = 1.00001*miss_val;
         miss_val_up = 0.99999*miss_val;
      else
         miss_val_low = 0.99999*miss_val;
         miss_val_up = 1.00001*miss_val;
      end

      if pos_vr > 0 | pos_min > 0
         if miss_val_up < miss_low
            need_index_m = 0;
         end
      end
      if pos_vr > 0 | pos_max > 0
         if miss_val_low > miss_up
            need_index_m = 0;
         end
      end
      if need_index_m
         index__FillValue = find ( miss_val_low <= values & ...
                                      values <= miss_val_up );
      end
   end
end

% Combine the arrays of missing value indices into one ordered array.
% Note that for real numbers the range of the _FillValue may intersect
% both the valid and invalid range there may be some indices that
% appear twice; this does not cause any inaccuracy, although it will
% result in some inefficiency.  In particular, rescaling is done on
% the set of indices NOT in index_miss and so is not affected.

index_miss = [ index_miss_low; index__FillValue; index_miss_up ];
index_miss = sort(index_miss);
len_index_miss = length(index_miss);

% If there are any missing values then offer to change them to a
% more convenient value.

if len_index_miss > 0
   s = [ varnam ' contains missing values:  Choose an action' ];
   s1 = 'Leave the missing value unchanged';
   s2 = 'Replace the missing value with NaN';
   s3 = 'Replace the missing value with a new value';
   k = -1;
   while any(k == [1 2 3]) == 0
      k = menu(s, s1, s2, s3);
      if k == 1
      elseif k == 2
         if length(index_miss) == 1
           values(index_miss) = NaN;
         else
            values(index_miss) = NaN*ones(index_miss);
         end
      elseif k == 3
         s = '   Type in your new missing value marker [0]  ';
         new_miss = return_val(s, 0);
         if length(index_miss) == 1
           values(index_miss) = new_miss;
         else
            values(index_miss) = new_miss*ones(index_miss);
         end
      else
         disp(' ')
         disp('You have asked for a non-existent option - try again')
      end
   end
end

% Close the netcdf file.

[rcode] = mexcdf('ncclose', cdfid);

