#!/usr/bin/env ruby

###############################
# MM5 file reader for Ruby    #
# written by Shigenori OTSUKA #
# from 2004.5.22              #
# last midified 2005.3.21     #
###############################

require "narray"

#CLASS FOR MM5===========================================================
class MM5Data
  class << self
    alias open new
  end

  def initialize(file,mode=nil)
    @mm5file = nil
    @bh = Big_header.new
    @sh = Sub_header.new
    
    @data = nil
    @current_t = 1
    @xtime = nil
    @variables = []
    @pos_time = []
    
    #file open
    @mm5file = File.open(file,mode="r")
  end
  
  #file close
  def close
    @mm5file.close
  end
  
  #open a variable
  def var(varname)
    return MM5Var.new(self, varname)
  end

  #read iheader
  def iheader
    if @mm5file.eof?
      @mm5file.rewind 
      @current_t = 1
    end
    @mm5file.seek(4,whence=IO::SEEK_CUR)
    ihd = @mm5file.read(4).unpack("N")[0]
    @mm5file.seek(4,whence=IO::SEEK_CUR)
    return ihd
  end

  #seek top of a time "t"
  def seek_time(t)
    if ! @pos_time[t]
      if !(t>@current_t)
        @mm5file.rewind
        @current_t = 1
      end
      begin
        case iheader
        when 0
          @bh.skip(@mm5file)
          if t==1
            @pos_time[t]=@mm5file.pos
            return
          end
        when 1
          @sh.skip(@mm5file)
          skip_field
        when 2
          @current_t += 1
          if !(t>@current_t)
            @pos_time[t]=@mm5file.pos
            return
          end
        else
          print "error\n"
          break
        end
      end until @mm5file.eof?
      @mm5file.rewind
      @current_t = 1
    else
      @mm5file.seek(@pos_time[t],whence=IO::SEEK_SET)
    end
  end
  
  #seek top of a field "var"
  def seek_var(varname,t)
    seek_time(t)
    @pos_var={} if !defined?(@pos_var)
    @pos_var[varname]=[] if @pos_var[varname]==nil
    if @pos_var[varname][t]!=nil
      @mm5file.seek(@pos_var[varname][t],whence=IO::SEEK_SET)
      @sh.read(@mm5file)
      return
    else
      begin
        case iheader
        when 0
          @bh.skip(@mm5file)
        when 1
          pos = @mm5file.pos
          @sh.read(@mm5file)
          if ((varname+" "*(@sh.name.size-varname.size))==@sh.name) 
            @current_var = varname
            @pos_var[varname][t] = pos
            return
          end
          skip_field
        when 2
          @current_t += 1
          print "#{varname}(#{t}): not found\n" 
          return nil
        else
          print "error\n"
          break
        end
      end until @mm5file.eof?
      @mm5file.rewind
      @current_t = 1
      return nil
    end
  end
  
  #read data "var" of time "t"
  def read(t,varname)
    seek_var(varname,t)
    idim = @sh.end_index[0...(@sh.ndim)] - @sh.start_index[0...(@sh.ndim)] + 1
    idim = [idim] if idim.is_a?(Fixnum) 
    if @sh.ordering[0..1]=="YX"
      if @sh.ndim == 2 then idim = [idim[1],idim[0]] end
      if @sh.ndim == 3 then idim = [idim[1],idim[0],idim[2]] end
    end
    data = NArray.sfloat(*idim)
    read_field(data)
    if (@sh.staggering == "C   ") and (@sh.ndim > 1)
      case @sh.ndim
      when 2
        return data[0..-2,0..-2]
      when 3
        return data[0..-2,0..-2,true]
      end
    else
      return data
    end
  end
  
  def read_partial(t,varname,sx,ex,sy,ey,sz,ez,data)
    #seek_time(t)
    #seek_var(varname)
    seek_var(varname,t)
    data = NArray.sfloat(ex-sx+1,ey-sy+1,ez-sz+1)
    read_field_partial(sx,sy,sz,@sh,data)
  end
  
  def file
    @mm5file
  end
  
  def filename
    @mm5file.path
  end
  
  def ntime # return a number of time (of "varname")
    @mm5file.rewind
    @current_t = 1
    t = 0
    if ! @ntime 
      begin
        case iheader
        when 0
          @bh.skip(@mm5file)
          @pos_time[@current_t]=@mm5file.pos if @pos_time[@current_t]==nil
        when 1
          @sh.skip(@mm5file)
          skip_field
        when 2
          @current_t += 1
        else
          print "error\n"
          break
        end
      end until @mm5file.eof?
      @ntime = @current_t - 1 
      @mm5file.rewind
      @current_t = 1
    end
    return @ntime
  end

  def times_var(varname)
    @times_var = {} if ! defined?(@times_var)
    if ! @times_var[varname]
      @times_var[varname]=[]
      @mm5file.rewind
      @current_t = 1
      t = 0
      begin
        case iheader
        when 0
          @bh.skip(@mm5file)
        when 1
          pos = @mm5file.pos
          @sh.read(@mm5file)
          if ((varname+" "*(@sh.name.size-varname.size))==@sh.name) 
            t += 1
            @times_var[varname].push @current_t
            @pos_var={} if !defined?(@pos_var)
            @pos_var[varname]=[] if @pos_var[varname]==nil
            @pos_var[varname][@current_t]=pos if @pos_var[varname][@current_t]==nil
          end
          skip_field
        when 2
          @current_t += 1
        else
          print "error\n"
          break
        end
      end until @mm5file.eof?
      @mm5file.rewind
      @current_t = 1
    end
    return @times_var[varname]
  end

  def ntime_var(varname)
    times_var(varname).length
  end

  def xtimes # read all xtimes
    if ! @xtimes
      @xtimes = []
      @mm5file.rewind
      @current_t = 1
      begin
        case iheader
        when 0
          @bh.skip(@mm5file)
        when 1
          @sh.read(@mm5file)
          skip_field
        when 2
          @xtimes[@current_t - 1] = @sh.xtime
          @current_t += 1
        else
          print "error\n"
          break
        end
      end until @mm5file.eof?
      t = @current_t - 1
      @mm5file.rewind
      @current_t = 1
    end
    return @xtimes
  end
  
  #return names of available variables at time "t"
  def variables(t)
    if ! @variables[t]
      n = 0
      var = []
      seek_time(t)
      begin
        case iheader
        when 0
          @bh.skip(@mm5file)
        when 1
          pos = @mm5file.pos
          @sh.read(@mm5file)
          #the variables with staggering CA are temporary unavailable
          varname = sh.name.strip!
          case sh.ordering
          when "CA  ","S   "
          else
            var[n] = varname
            n += 1
          end
          @pos_var={} if !defined?(@pos_var)
          @pos_var[varname]=[] if @pos_var[varname]==nil
          @pos_var[varname][t]=pos if @pos_var[varname][t]==nil
          skip_field
        when 2
          @current_t += 1
          @variables[t] = var
          return var
        else
          print "error\n"
          break
        end
      end until @mm5file.eof?
      @mm5file.rewind
      @current_t = 1
    else
      return @variables[t]
    end
  end

  #return names of available variables in this file
  def var_names
    if ! defined?(@var_names)
      nt = ntime
      names = []
      for i in 0...nt
        names[i] = variables(i+1)
      end
      @var_names = names.flatten.uniq
    end
    return @var_names
  end
  alias varnames var_names

  def exist?(name)
    var_names.each{|var|
      return true if var==name
    }
    return false
  end

  def each_att
    p "each_att not supported"
  end

  #read big header of this file
  def bh_read
    if ! @bhreadflag
      @mm5file.rewind
      @current_time = 1
	begin
          case iheader
          when 0
            @bh.read(@mm5file)
            return
          when 1
            @sh.skip(@mm5file)
            skip_field
          when 2
            @current_t += 1
          else
            print "error\n"
            break
          end
	end until @mm5file.eof?
      @mm5file.rewind
      @current_time = 1
      @bhreadflag = true
    end
  end
  
  def bh
    if ! @bhreadflag ; bh_read ; end
    @bh
  end
  
  def sh
    @sh
  end

  def current_var
    @current_var
  end

  #read field data of current filepoint
  def read_field(data)
    @mm5file.seek(4,whence=IO::SEEK_CUR)
    case data.dim
    when 1
      data[true] = @mm5file.read(data.size*4).unpack("g*")
    when 2
      if @sh.ordering[0..1]=="YX"
        data[0,0] = NArray.to_na(@mm5file.read(data.size*4),NArray::SFLOAT).ntoh.reshape!(*data.sizes[0..1].reverse).transpose(1,0)
      else
        data[0,0] = NArray.to_na(@mm5file.read(data.size*4),NArray::SFLOAT).ntoh.reshape!(*data.sizes[0..1])
      end
    when 3
      if @sh.ordering[0..1]=="YX"
        data[0,0,0] = NArray.to_na(@mm5file.read(data.size*4),NArray::SFLOAT).ntoh.reshape!(data.sizes[1],data.sizes[0],data.sizes[2]).transpose(1,0,2)
      else
        data[0,0,0] = NArray.to_na(@mm5file.read(data.size*4),NArray::SFLOAT).ntoh.reshape!(*data.sizes[0..2])
      end
    end
    @mm5file.seek(4,whence=IO::SEEK_CUR)	
  end
  
  def read_field_partial(sx,sy,sz,sh,data)
    @mm5file.seek(4,whence=IO::SEEK_CUR)
    @mm5file.seek( (sz-1)*sh.end_index[1]*sh.end_index[0]*4 ,whence=IO::SEEK_CUR)
    for k in 0..(data.sizes[2]-1)
      @mm5file.seek( (sx-1)*sh.end_index[0]*4 ,whence=IO::SEEK_CUR)
      for j in 0..(data.sizes[0]-1)
        @mm5file.seek( (sy-1)*4 ,whence=IO::SEEK_CUR)
        data[j,0..-1,k] = @mm5file.read(data.sizes[1]*4).unpack("g*")
        @mm5file.seek( (sh.end_index[0]-sy-data.size[1]+1)*4 ,whence=IO::SEEK_CUR)
      end
      @mm5file.seek( (sh.end_index[1]-sx-data.size[0]+1)*sh.end_index[0]*4 ,whence=IO::SEEK_CUR)
    end
    @mm5file.seek( (sh.end_index[2]-sz-data.size[2]+1)*sh.end_index[1]*sh.end_index[0]*4 ,whence=IO::SEEK_CUR)
    @mm5file.seek(4,whence=IO::SEEK_CUR)
  end
  
  def skip_field
    size = file.read(4).unpack("N")[0]
    @mm5file.seek(size+4,whence=IO::SEEK_CUR)
  end
end

#CLASS FOR BIG HEADER=============================================================
class Big_header
  def initialize
    @bhi = NArray.int(50,20)
    @bhr = NArray.sfloat(20,20)
  end
  
  def read(file)
    file.seek(4,whence=IO::SEEK_CUR)
    # read BHI
    @bhi[true,true] = NArray.to_na(file.read(4*50*20),NArray::INT).ntoh.reshape!(50,20)
    # read BHR
    @bhr[true,true] = NArray.to_na(file.read(4*20*20),NArray::SFLOAT).ntoh.reshape!(20,20)
    file.seek(50*20*80,whence=IO::SEEK_CUR) # for BHIC
    file.seek(20*20*80,whence=IO::SEEK_CUR) # for BHRC
    file.seek(4,whence=IO::SEEK_CUR)
  end
  
  def skip(file)
    size = file.read(4).unpack("N")[0]
    file.seek(size+4,whence=IO::SEEK_CUR)
  end
  
  def bhi
    @bhi
  end
  
  def bhr
    @bhr
  end
end

# CLASS FOR SUB HEADER===================================================
class Sub_header
  def initialize
    @ndim = 0
    @start_index = NArray.int(4)
    @end_index = NArray.int(4)
    @xtime = 0.0
    @staggering = nil
    @ordering = nil
    @current_date = nil
    @name = nil
    @unit = nil
    @description = nil
  end
  
  def read(file)
    file.seek(4,whence=IO::SEEK_CUR)
    # read sub header
    @ndim = file.read(4).unpack("N")[0]
    for i in 0..3
      @start_index[i] = file.read(4).unpack("N")[0]
    end
    for i in 0..3
      @end_index[i] = file.read(4).unpack("N")[0]
    end
    @xtime = file.read(4).unpack("g")[0]
    @staggering = file.read(4)
    @ordering = file.read(4)
    @current_time = file.read(24)
    @name = file.read(9)
    @unit = file.read(25)
    @description = file.read(46)
    file.seek(4,whence=IO::SEEK_CUR)	
  end
  
  def skip(file)
    size = file.read(4).unpack("N")[0]
    file.seek(size+4,whence=IO::SEEK_CUR)
  end
  
  def ndim
    @ndim
  end
  
  def start_index
    @start_index
  end
  
  def end_index
    @end_index
  end
  
  def xtime
    @xtime
  end
  
  def staggering
    @staggering
  end
  
  def ordering
    @ordering
  end
  
  def name
    @name
  end
end

#CLASS FOR VARIABLE===============================================================================
class MM5Var
  def initialize(file,varname)
    @varname = varname
    if file.is_a?(String)
      file = MM5Data.open(file)
    elsif ! file.is_a?(MM5Data)
      raise ArgumentError, "1st arg must be a MM5data or a file name"
    end
    @mm5data = file

    # set dimensions of the variable
    @dimensions = []
    @mm5data.seek_time(1)
    if @varname == "t"
      ntime = @mm5data.ntime
      @ordering = ["T"]# @mm5data.sh.ordering.dup
      @staggering = "C"# @mm5data.sh.staggering.dup[0..0]
      @dimensions[0] = ntime
    elsif @varname == "x"
      ntime = 1
      @mm5data.seek_var("LONGICRS",1)
      @start_index = @mm5data.sh.start_index.dup
      @end_index = @mm5data.sh.end_index.dup
      @ordering = @mm5data.sh.ordering.dup
      @staggering = @mm5data.sh.staggering.dup[0..0]
      @dimensions[0] = @end_index[1] - @start_index[1] +1 -1
    elsif @varname == "xdot"
      ntime = 1
      @mm5data.seek_var("LONGICRS",1)
      @start_index = @mm5data.sh.start_index.dup
      @end_index = @mm5data.sh.end_index.dup + 1
      @ordering = @mm5data.sh.ordering.dup
      @staggering = @mm5data.sh.staggering.dup[0..0]
      @dimensions[0] = @end_index[1] - @start_index[1] +1 -1
    elsif @varname == "y"
      ntime = 1
      #p "test"
      #if @mm5data.exist?("LATITCRS")
        @mm5data.seek_var("LATITCRS",1) 
        @start_index = @mm5data.sh.start_index.dup
        @end_index = @mm5data.sh.end_index.dup
        @ordering = @mm5data.sh.ordering.dup
        @staggering = @mm5data.sh.staggering.dup[0..0]
        @dimensions[0] = @end_index[0] - @start_index[0] +1 -1
      #elsif @mm5data.sh.ordering[0..0]=="Y"
      #  @start_index = @mm5data.sh.start_index.dup[0]
      #  @end_index = @mm5data.sh.end_index.dup[0]
      #  @ordering = "Y"
      #  @staggering = @mm5data.sh.staggering.dup[0..0]
      #  @dimensions[0] = @end_index[0] - @start_index[0] +1 -1
      #end
    elsif @varname == "ydot"
      ntime = 1
        #p "test",@mm5data.sh.ordering #[0..0]
        #p @mm5data.current_var
        #p @mm5data.sh.start_index.dup[0]
        #p @mm5data.sh.end_index.dup[0]

        # @start_index = @mm5data.sh.start_index.dup[0]
        # @end_index = @mm5data.sh.end_index.dup[0]
        # @ordering = "Y"
        # @staggering = @mm5data.sh.staggering.dup[0..0]
        # @dimensions = [@end_index[0] - @start_index[0] +1 -1]
      #if @mm5data.exist?("LATITCRS")
        @mm5data.seek_var("LATITCRS",1)
        @start_index = @mm5data.sh.start_index.dup
        @end_index = @mm5data.sh.end_index.dup + 1
        @ordering = @mm5data.sh.ordering.dup
        @staggering = @mm5data.sh.staggering.dup[0..0]
        @dimensions[0] = @end_index[0] - @start_index[0] +1 -1
      #elsif @mm5data.sh.ordering[0..0]=="Y"
        # @start_index = @mm5data.sh.start_index.dup[0]
        # @end_index = @mm5data.sh.end_index.dup[0]
        # @ordering = "Y"
        # @staggering = @mm5data.sh.staggering.dup[0..0]
        # @dimensions = [@end_index[0] - @start_index[0] +1 -1]
      #end
    elsif @varname == "z" or @varname == "SIGMAH"
      ntime = 1
      @mm5data.seek_var("SIGMAH",1)
      @start_index = @mm5data.sh.start_index.dup
      @end_index = @mm5data.sh.end_index.dup
      @dimensions[0] = @end_index[0] - @start_index[0] + 1
      @ordering = @mm5data.sh.ordering.dup
      @staggering = @mm5data.sh.staggering.dup[0..0]
    elsif @varname == "sigmafull"
      ntime = 1
      @mm5data.seek_var("SIGMAH",1)
      @start_index = @mm5data.sh.start_index.dup
      @end_index = @mm5data.sh.end_index.dup + 1
      @ordering = @mm5data.sh.ordering.dup
      @staggering = @mm5data.sh.staggering.dup[0..0]
      @dimensions[0] = @end_index[0] - @start_index[0] + 1
    elsif @varname == "PRESSURE"
      ntime = 1
      @mm5data.seek_var("PRESSURE",1)
      @start_index = @mm5data.sh.start_index.dup
      @end_index = @mm5data.sh.end_index.dup
      @dimensions[0] = @end_index[0] - @start_index[0] + 1
      @ordering = @mm5data.sh.ordering.dup
      @staggering = @mm5data.sh.staggering.dup[0..0]
    elsif @varname == "boundary"
      ntime = 1
      # @mm5data.seek_var("PRESSURE")
      @start_index = [0] # @mm5data.sh.start_index.dup
      @end_index = [4] # @mm5data.sh.end_index.dup
      @dimensions[0] = @end_index[0] - @start_index[0] + 1
      @ordering = "boundary" # @mm5data.sh.ordering.dup
      @staggering = nil # @mm5data.sh.staggering.dup[0..0]
    elsif @varname == "landusecategory" 
      raise "landusecategory is not supported yet."
    elsif @varname == "landuse2" 
      raise "landuse2 is not supported yet."
    else
      @mm5data.seek_var(varname,1)
      @start_index = @mm5data.sh.start_index.dup
      @end_index = @mm5data.sh.end_index.dup
      for i in 0...(@mm5data.sh.ndim)
        @dimensions[i]=@end_index[i] - @start_index[i] + 1
      end
      # if order of axes is Y-X, convert to X-Y
      @ordering = @mm5data.sh.ordering.dup
      if @ordering[0..1]=="YX"
        @dimensions[0..1] = @dimensions[0..1].reverse!
      end
      # if cross point, cut off dummy cells
      @staggering = @mm5data.sh.staggering.dup[0..0]
      if @staggering == "C"
        @dimensions[0] = @dimensions[0] - 1
        @dimensions[1] = @dimensions[1] - 1
      end 
      # Add time axis
      #ntime = @mm5data.ntime(@varname)
      ntime = @mm5data.ntime_var(@varname)
      @dimensions.push(ntime)
    end
    @rank = @dimensions.length
    #p [@varname,@dimensions,@rank]
    # Set Attribute (still not working now)
    @attr = Attribute.new
  end

  def rank
    @rank
  end

  def dim_names
    if ! defined?(@dim_names)
      case @ordering
      when "YXP "
        if @staggering[0..0]=="C"
          @dim_names = ["x","y","PRESSURE","t"]
        else
          @dim_names = ["xdot","ydot","PRESSURE","t"]
        end
      when "YXS "
        if @staggering[0..0]=="C"
          @dim_names = ["x","y","SIGMAH","t"]
        else
          @dim_names = ["xdot","ydot","SIGMAH","t"]
        end
      when "YXW "
        if @staggering[0..0]=="C"
          @dim_names = ["x","y","sigmafull","t"]
        else
          @dim_names = ["xdot","ydot","sigmafull","t"]
        end
      when "YX  "
        if @staggering[0..0] == "C"
          @dim_names = ["x","y","t"]
        else
          @dim_names = ["xdot","ydot","t"]
        end
      when "CA  "
        @dim_names = ["landusecategory","landuse2"]
      when "XSB "
        if @staggering[0..0]=="C"
          @dim_names = ["x","SIGMAH","boundary","t"]
        else
          @dim_names = ["xdot","SIGMAH","boundary","t"]
        end
      when "YSB "
        if @staggering[0..0]=="C"
          @dim_names = ["y","SIGMAH","boundary","t"]
        else
          @dim_names = ["ydot","SIGMAH","boundary","t"]
        end
      when "XWB "
        if @staggering[0..0]=="C"
          @dim_names = ["x","sigmafull","boundary","t"]
        else
          @dim_names = ["xdot","sigmafull","boundary","t"]
        end
      when "YWB "
        if @staggering[0..0]=="C"
          @dim_names = ["y","sigmafull","boundary","t"]
        else
          @dim_names = ["ydot","sigmafull","boundary","t"]
        end
      when "S   "
        @dim_names = ["SIGMAH"]
      when "P   "
        @dim_names = ["PRESSURE"]
      when "boundary"
        @dim_names = ["boundary"]
      else
        @dim_names = ["t"] if @varname == "t"
      end
    end
    return @dim_names
  end

  def name
    @varname
  end

  def name=(name=nil)
    p "name= not supported"
  end

  def vartype
    "sfloat"
  end

  def typecode
    4
  end

  def att(name)
    @attr[name]
  end

  def natts
    @attr.length
  end

  def att_names
    @attr.keys
  end

  def each_att
    p "each_att not supported"
  end

  def attr
    @attr
  end

  def file
    @mm5data.filename
  end

  def staggering
    @staggering
  end

  def ordering
    @ordering
  end

  def shape_ul0
    shape = []
    for i in 0...(@dimensions.length)
      if self.dim_names[i]=="t"
        shape.push(0)
      else
        shape.push(@dimensions[i])
      end
    end
    return shape
  end

  def shape_current
    shape = []
    @dimensions.each{|dim|
      #shape.push(dim[:len])
      shape.push(dim)
    }
    return shape
  end

  def get(*a)
    #p ["a",a]
    ary = []
    @dimensions.each{|d|
      ary.push((0...(d)).to_a)
    }
    if a==nil
    elsif a.is_a?(Array)
      for i in 0...(ary.length)
        if a[i].is_a?(Range)
          ary[i] = NArray.int(@dimensions[i]).indgen![a[i]]
        elsif a[i]==true
        elsif a[i]==false
          break
        elsif a[i].is_a?(Array) or a[i].is_a?(NArray)
          ary[i] = NArray.int(@dimensions[i]).indgen![a[i]].to_a
        elsif a[i].is_a?(Fixnum)
          ary[i] = [NArray.int(@dimensions[i]).indgen![a[i]]]
        end
      end
    elsif a.is_a?(Fixnum)
      ary[i] = [NArray.int(@dimensions[i]).indgen![a]]
    end
    case @varname
    #variables as axes
    when "t"
      if @mm5data.xtimes.uniq.length != @mm5data.xtimes.length
        return NArray.sfloat(@mm5data.xtimes.length).indgen![(ary[0])]
      else 
        return NArray.to_na(@mm5data.xtimes)[(ary[0])]
      end
    when "x","xdot","y","ydot"
      return NArray.sfloat(@dimensions[0]).indgen!(0.0,@mm5data.bh.bhr[8,0])[(ary[0])]
    when "z"
      return NArray.sfloat(@dimensions[0]).indgen![(ary[0])]
    when "SIGMAH"
      return @mm5data.read(1,"SIGMAH")[(ary[0])]
    when "sigmafull"
      tmp = @mm5data.read(1,"SIGMAH")
      return NArray.to_na([0.0,((tmp[0..-2]+tmp[1..-1])/2.0).to_a,1.0].flatten!)[(ary[0])]
    when "PRESSURE"
      return @mm5data.read(1,"PRESSURE")[(ary[0])]
    #ordinary variables
    else
      dim = []
      ary.each{|b| dim.push(b.length)}
      na = NArray.sfloat(*dim)
      for i in 0...(ary[-1]).length
	case @dimensions.length
	when 3
	  na[true,true,i] = @mm5data.read((ary[-1][i])+1, @varname)[*ary[0..-2]]
	when 4
	  na[true,true,true,i] = @mm5data.read((ary[-1][i])+1, @varname)[*ary[0..-2]]
	end
      end
      return na
    end
  end

  def [](*a)
    # temporary routine
    #return self.get[*a]
    return self.get(*a)

    # old routine
    # copied from grads_gridded.rb by S. OTSUKA
    # copied from NetCDFVar class by R. MIZUTA
    #print "[]="
    #first = Array.new
    #last = Array.new
    #stride = Array.new
    #set_stride = false
    #a.each{|i|
    #  if(i.is_a?(Fixnum))
    #    first.push(i)
    #    last.push(i)
    #    stride.push(1)
    #  elsif(i.is_a?(Range))
    #    first.push(i.first)
    #    last.push(i.exclude_end? ? i.last-1 : i.last)
    #    stride.push(1)
    #  elsif(i.is_a?(Hash))
    #    r = (i.to_a[0])[0]
    #    s = (i.to_a[0])[1]
    #    if ( !( r.is_a?(Range) ) || ! ( s.is_a?(Integer) ) )
    #      raise TypeError, "Hash argument must be {a_Range, step}"
    #    end
    #    first.push(r.first) 
    #    last.push(r.exclude_end? ? r.last-1 : r.last)
    #    stride.push(s)
    #    set_stride = true
    #  elsif(i.is_a?(TrueClass))
    #    first.push(0)
    #    last.push(-1)
    #    stride.push(1)
    #  end
    #}
    #if(set_stride)
    #  na = self.get({"start"=>first, "end"=>last, "stride"=>stride})
    #else
    #  na = self.get({"start"=>first, "end"=>last})
    #end
    #shape = na.shape
    #(a.length-1).downto(0){ |i|
    #  #shape.delete_at(i) if a[i].is_a?(Fixnum)
    #  shape.delete_at(i) if first[i]==last[i]
    #}
    #na.reshape!( *shape )
    #na
  end
  
  def []=(*a)  # copied from NetCDFVar class
    val = a.pop
    first = Array.new
    last = Array.new
    stride = Array.new
    set_stride = false
    a.each{|i|
      if(i.is_a?(Fixnum))
        first.push(i)
        last.push(i)
        stride.push(1)
      elsif(i.is_a?(Range))
        first.push(i.first)
        last.push(i.exclude_end? ? i.last-1 : i.last)
        stride.push(1)
      elsif(i.is_a?(Hash))
        r = (i.to_a[0])[0]
        s = (i.to_a[0])[1]
        if ( !( r.is_a?(Range) ) || ! ( s.is_a?(Integer) ) )
          raise ArgumentError, "Hash argument must be {first..last, step}"
        end
        first.push(r.first) 
        last.push(r.exclude_end? ? r.last-1 : r.last)
        stride.push(s)
        set_stride = true
      elsif(i.is_a?(TrueClass))
        first.push(0)
        last.push(-1)
        stride.push(1)
      end
    }
    
    if(set_stride)
      self.put(val, {"start"=>first, "end"=>last, "stride"=>stride})
    else
      self.put(val, {"start"=>first, "end"=>last})
    end
  end

end