#########################
# GS2 CodeRunner module
# Data Analysis
#########################
#

class NumRu::NetCDF
  aliold :var
  def var(*args)
    if args[0].kind_of? Symbol
      args[0] = args[0].to_s
    end
    return old_var(*args)
  end
end

class CodeRunner
class Gs2

eval(File.read(File.dirname(__FILE__) + '/gsl_tools.rb'), GLOBAL_BINDING, File.dirname(__FILE__) + '/gsl_tools.rb')

# def gsl_vector(name, options={})
#   if options[:t_index] or options[:frame_index] or not [:Failed, :Complete].include? status
#     return get_gsl_vector(name, options)
#   else
#     return cache[[:gsl_vector, name, options]] ||= get_gsl_vector(name, options)
#   end
# end

def netcdf_file
  #if @runner.cache[:runs] and (open = @runner.cache[:runs].keys.find_all{|id| @runner.cache[:runs][id][:netcdf_file]}).size > 200
  #ep "my id", id
  if (open = @runner.run_list.keys.find_all{|id|  @runner.run_list[id].cache[:netcdf_file]}).size > 200
    open = open.sort_by{|id| @runner.run_list[id].cache[:netcdf_file_otime]}
    @runner.run_list[open[0]].ncclose
  end

  if cache[:netcdf_file] and not [:Complete, :Failed].include? @status
    ncclose
  end
  cache[:netcdf_file_otime] = Time.now.to_i
  cache[:netcdf_file] ||= NumRu::NetCDF.open(netcdf_filename)
  cache[:netcdf_file].sync
  cache[:netcdf_file]
end

def netcdf_filename
  @directory + '/' +  @run_name + '.out.nc'
end

def ncclose
  cache[:netcdf_file].close
  cache.delete(:netcdf_file)
end

module FixNormOption
#class << self
    def fix_norm_action(options)
      case options[:set_norm_option]
      when "t_over_m", "bd"
        if ["t_over_m", "bd"].include? @norm_option
          return :none
        else
          return :from_root_2
        end
      else
        #eputs "else", norm_option
        if ["t_over_m", "bd"].include? @norm_option
          #eputs "norm old"
          return :to_root_2
        else
          return :none
        end
      end
    end

    def fix_heat_flux_norm(tensor, options={})
      case fix_norm_action(options)
      when :none
        #eputs "none"
        return tensor
      when :to_root_2
        eputs "to_root_2"
        return tensor / 2.0**1.5
      when :from_root_2
        return tensor * 2.0**1.5
      end
    end

    # Return tensor normalised according to options[:set_norm_option]
    # (which may be "t_over_m", "bd" or "with_root_2", "mtk" etc)
    # regardless of the original normalisation.
    #
    # <tt>power</tt> should be the power to which the reference thermal
    # velocity is raised in the normalising quantity. For example,
    # <tt>t</tt> is normalised to a / v_thr, so for times, power should
    # be set equal to -1.

    def fix_norm(tensor, power, options={})
      case fix_norm_action(options)
      when :none
        #eputs "none"
        return tensor
      when :to_root_2
        eputs "to_root_2"
        return tensor / 2.0**(0.5 * power)
      when :from_root_2
        return tensor * 2.0**(0.5 * power)
      end
    end

#end # class << self
end # module FixNormOption

include FixNormOption

def gsl_vector(name, options={})
  Dir.chdir(@directory) do
    options[:t_index_window] ||= @scan_index_window
    options.setup_time_window
    if [:ky, :kx].include? name.to_sym
      vec = fix_norm(
        GSL::Vector.alloc(netcdf_file.var(name.to_s).get.to_a.sort),
        -1, options
      ) # ky, ky are normalised to 1 / rho_i
      if i = options[:interpolate_ + name.to_s.sub(/k/, '').to_sym]
        if name.to_sym == :ky
          s = (vec.size - 1)*i + 1
          #return vec.connect(GSL::Vector.alloc((vec.size-1)*(i-1)) * 0.0)
          return (0...s).map{|k| k.to_f * vec[1]}.to_gslv
        else
          size = vec.size
          #vec = vec.to_box_order
          raise "Hmmm, kx.size should be odd" unless size%2 == 1
          s = (size-1)/2 * i
          return (-s..s).to_a.map{|ii| ii.to_f * vec.to_box_order[1]}.to_gslv
          #new_vec = GSL::Vector.alloc((s-1)*i + 1)
          #new_vec *= 0.0
          #for j in 0...((s-1)/2+1)
            #new_vec[j] = vec[j]
          #end
          #for j in 0...((s-1)/2)
            #new_vec[-j-1] = vec[-j-1]
          #end
          #return new_vec.from_box_order
        end


      else
        return vec
      end
    elsif [:theta].include? name.to_sym
      #ep options; gets
      #vec = GSL::Vector.alloc(netcdf_file.var(name.to_s).get({'start' => [options[:thetamin]||0], 'end' => [options[:thetamax]||-1]}).to_a)
      vec = GSL::Vector.alloc(netcdf_file.var(name.to_s).get.to_a)
      if gryfx? and options[:periodic]
        #vec = vec.connect([2.0*vec[-1] - vec[-2]].to_gslv)
        vec = vec.connect([-vec[0]].to_gslv)
      end
      if ith = options[:interpolate_theta]
        osize = vec.size
        newsize = (osize-1)*ith+1
        newvec = GSL::Vector.alloc(newsize)
        newvec[newsize-1] = vec[osize-1]# * ith.to_f
        for i in 0...(newsize-1)
          im = i%ith
          frac = im.to_f/ith.to_f
          #iold = (i-im)/(new_shape[-1]-1)*(shape[-1]-1)
          iold = (i-im)/ith
          newvec[i] =  (vec[iold] * (1.0-frac) + vec[iold+1] * frac)
        end
        vec = newvec
      end
      start = options[:thetamin]||0
      endv = options[:thetamax]||vec.size-1
      #ep ['options', options, 'vec.size', vec.size]
      vec = vec.subvector(start, (endv-start+1)).dup
      return vec
    elsif name.to_sym == :t
      #options.setup_time_window
      t = GSL::Vector.alloc(netcdf_file.var(name.to_s).get('start' => [options[:begin_element]], 'end' => [options[:end_element]]).to_a)
      t = t - t[0] if options[:sync_time]
      return fix_norm(t, -1, options) # t is normalised to a/v_thi
    end
    options = eval(options) if options.class == String
    if options[:saturated_time_average] or options[:sta]
      raise "Not Saturated" unless @saturation_time_index
      tmax = list(:t).keys.max
      return ((@saturation_time_index..tmax).to_a.map do |t_index|
        gsl_vector(name, options.dup.absorb({t_index: t_index, saturated_time_average: nil, sta: nil}))
      end).sum / (list(:t).values.max - list(:t)[@saturation_time_index])
    elsif options[:time_average] or options[:ta]
      tmax = list(:t).keys.max
      start_t = 2
      return ((start_t..tmax).to_a.map do |t_index|
        gsl_vector(name, options.dup.absorb({t_index: t_index, time_average: nil, ta: nil}))
      end).sum / (list(:t).values.max - list(:t)[start_t])
    end
    if method = self.class.instance_methods.find{|meth| (name + '_gsl_vector').to_sym == meth}
      options[:graphkit_name] = name
      return send(method, options)
    end
  end
  raise "GSL Vector #{name} not found"
end

module GSLVectors

  # The square of the potential summed over all wave numbers, indexed by time, normalised to (e/T)(rho_1/a).

  def phi2tot_over_time_gsl_vector(options)

    Dir.chdir(@directory) do      #Necessary options: ky
      #log 'about to open netcdf file'
      #options.setup_time_window
      phis = netcdf_file.var('phi2').get('start'=>[options[:begin_element]], 'end'=>[options[:end_element]] ).to_a
      log 'about to allocate gsl vector'
      vec = GSL::Vector.alloc(phis)
      log 'finished'
      return fix_norm(vec, 1, options)
    end
  end

  def apar2_over_time_gsl_vector(options)

    Dir.chdir(@directory) do      #Necessary options: ky
      #log 'about to open netcdf file'
      #options.setup_time_window
      phis = netcdf_file.var('apar2').get('start'=>[options[:begin_element]], 'end'=>[options[:end_element]] ).to_a
      log 'about to allocate gsl vector'
      vec = GSL::Vector.alloc(phis)
      log 'finished'
      return fix_norm(vec, 1, options)
    end
  end

  def transient_es_heat_flux_amplification_over_kx_gsl_vector(options)
    options[:direction] = :kx
    transient_es_heat_flux_amplification_over_kxy_gsl_vector(options)
  end

  def transient_es_heat_flux_amplification_over_ky_gsl_vector(options)
    options[:direction] = :ky
    transient_es_heat_flux_amplification_over_kxy_gsl_vector(options)
  end

  def transient_es_heat_flux_amplification_over_kxy_gsl_vector(options)
    Dir.chdir(@directory) do      # i.e. phi2_by_ky_vs_time or phi2_by_kx_vs_time
      kxy = options[:direction].to_sym

#       ep :growth_rate_at_ + kxy
      p send(:transient_es_heat_flux_amplification_at_species_at_ + kxy)
      return GSL::Vector.alloc(send(:transient_es_heat_flux_amplification_at_species_at_ + kxy)[options[:species_index]-1].values)

    end
  end

  def transient_amplification_over_kx_gsl_vector(options)
    options[:direction] = :kx
    transient_amplification_over_kxy_gsl_vector(options)
  end
  
  def transient_amplification_over_ky_gsl_vector(options)
    options[:direction] = :ky
    transient_amplification_over_kxy_gsl_vector(options)
  end
  
  def transient_amplification_over_kxy_gsl_vector(options)
    Dir.chdir(@directory) do      # i.e. phi2_by_ky_vs_time or phi2_by_kx_vs_time
      kxy = options[:direction]
#       ep :growth_rate_at_ + kxy
      return GSL::Vector.alloc(send(:transient_amplification_at_ + kxy).values)

    end
  end
  private :transient_amplification_over_kxy_gsl_vector

  # The growth rate of the fluctuations, calculated from the potential, indexed by time and normalised to vth_1/a.
  # :kx or :kx_index must be specified in options
  def growth_rate_by_kx_over_time_gsl_vector(options)
    options[:direction] = :kx
    growth_rate_by_kxy_over_time_gsl_vector(options)
  end

  # The growth rate of the fluctuations, calculated from the potential, indexed by time and normalised to vth_1/a.
  # :ky or :ky_index must be specified in options
  def growth_rate_by_ky_over_time_gsl_vector(options)
    options[:direction] = :ky
    growth_rate_by_kxy_over_time_gsl_vector(options)
  end
    
  def growth_rate_by_kxy_over_time_gsl_vector(options)
    # i.e. time_dependent_gr_by_ky_vs_time or phi2_by_kx_vs_time

    kxy = options[:direction]

    phi = gsl_vector("phi2_by_#{kxy}_over_time", options).log / 2.0

    size = phi.size
    dphi = phi.subvector(1, size - 1) - phi.subvector(0, size-1)
    # NB dt already has norm fixed, dphi is dimensionless
    return fix_norm(dphi/gsl_vector('dt'), 0, options)
  end

  # The real frequency of the fluctuations, read from the .out file, indexed by time and normalised to vth_1/a.
  # :ky_index or :kx_index must be specified in options.
  def frequency_by_kx_over_time_gsl_vector(options)
    options[:direction] = :kx
    frequency_by_kxy_over_time_gsl_vector(options)
  end

  def frequency_by_ky_over_time_gsl_vector(options)
    options[:direction] = :ky
    frequency_by_kxy_over_time_gsl_vector(options)
  end

  def frequency_by_kxy_over_time_gsl_vector(options)
    kxy = options[:direction]
    kxy_index = kxy + :_index
    kxys = get_list_of(kxy)
    desired_kxy = kxys[options[kxy_index]]
    raise "No k found at the desired index" if desired_kxy.nil?

    omega_reals = []
    File.open(@run_name+".out",'r') do |fileHandle|
      fileHandle.each_line do |fileLine|
        if fileLine.include?('aky=')  # Only examine the lines of the .out file that contain frequency information.

          index = fileLine.index('akx=')
          raise "akx wasn't found where it was expected in the .out file." if index.nil?
          akx = fileLine[(index+4)..-1].to_f

          index = fileLine.index('aky=')
          raise "aky wasn't found where it was expected in the .out file." if index.nil?
          aky = fileLine[(index+4)..-1].to_f

          index = fileLine.index('om=')
          raise "om wasn't found where it was expected in the .out file." if index.nil?
          omr = fileLine[(index+3)..-1].to_f
          if kxy == :kx
            # You need to be careful when testing equality of the desired k with the k in the .out file
            # since the .out file is only written to ~ 5 significant digits:
            omega_reals << omr if ((desired_kxy - akx).abs/(desired_kxy.abs + 1e-7) < 1e-4)
          else
            omega_reals << omr if ((desired_kxy - aky).abs/(desired_kxy.abs + 1e-7) < 1e-4)
          end
        end
      end
    end
    raise "No real frequencies found in the .out file for the desired k" if (omega_reals.size==0)
    GSL::Vector.alloc(omega_reals)
  end

  # The size of each time step,  indexed by time, normalised to a/v_th1.
  def dt_gsl_vector(options)
    t = gsl_vector('t', options)
    size = t.size
    # NB t already has norm fixed
    return t.subvector(1, size - 1) - t.subvector(0, size-1)
  end

  # The growth rate, calculated from the potential, indexed by kx. Only makes sense in linear calculations.
  def growth_rate_over_kx_gsl_vector(options)
    options[:direction] = :kx
    growth_rate_over_kxy_gsl_vector(options)
  end
  # The growth rate, calculated from the potential, indexed by ky. Only makes sense in linear calculations.
  def growth_rate_over_ky_gsl_vector(options)
    options[:direction] = :ky
    growth_rate_over_kxy_gsl_vector(options)
  end

  def growth_rate_over_kxy_gsl_vector(options)
    Dir.chdir(@directory) do      # i.e. phi2_by_ky_vs_time or phi2_by_kx_vs_time
      kxy = options[:direction]
      return GSL::Vector.alloc(send(:growth_rate_at_ + kxy).values)
    end
  end
  private :growth_rate_over_kxy_gsl_vector

  # The growth rate, calculated from the potential, indexed by kx. Only makes sense in linear calculations.
  def growth_rate_over_kx_slice_gsl_vector(options)
    Dir.chdir(@directory) do
      slice_of_growth_rates = send(:growth_rate_at_ky_at_kx)[options[:ky]].values
      raise "Something went wrong: slice of growth rates seems empty" if slice_of_growth_rates.nil?
      return GSL::Vector.alloc(slice_of_growth_rates)
      #return GSL::Vector.alloc(send(:growth_rate_at_ky_at_kx[ky]).values)
    end
  end

  # The growth rate, calculated from the potential, indexed by ky. Only makes sense in linear calculations.
  def growth_rate_over_ky_slice_gsl_vector(options)
    Dir.chdir(@directory) do
      slice_of_growth_rates = send(:growth_rate_at_ky_at_kx).values.map{|h| h[options[:kx]]}
      raise "Something went wrong: slice of growth rates seems empty" if slice_of_growth_rates.nil?
      return GSL::Vector.alloc(slice_of_growth_rates)
    end
  end

  def es_heat_flux_by_kx_over_time_gsl_vector(options)
    options[:direction] = :kx
    es_heat_flux_by_kxy_over_time_gsl_vector(options)
  end

  def es_heat_flux_by_ky_over_time_gsl_vector(options)
    options[:direction] = :ky
    es_heat_flux_by_kxy_over_time_gsl_vector(options)
  end

  def es_heat_flux_by_kxy_over_time_gsl_vector(options)
    raise "The write_fluxes_by_mode flag was not set." unless @write_fluxes_by_mode
    Dir.chdir(@directory) do
      kxy = options[:direction]
      #kxy_index = kxy + :_index
      options.convert_to_index(self, kxy)
      raise "Please provide species_index " unless options[:species_index]
      if kxy==:ky
        lkx = list(:kx)
        es_heat_av = (lkx.keys.map do |kx_index|
          es_heat =  netcdf_file.var('es_heat_flux_by_mode').get({'start' => [kx_index-1,options[:ky_index]-1,options[:species_index]-1, 0], 'end' => [kx_index-1,options[:ky_index]-1,options[:species_index]-1, -1]})
          #ep phi.shape
          es_heat.reshape(*es_heat.shape.values_at(3))
        end).sum / lkx.size
        return es_heat_av.to_gslv
      else
        lky = list(:ky)
        es_heat_av = (lky.keys.map do |ky_index|
          es_heat =  netcdf_file.var('es_heat_flux_by_mode').get({'start' => [options[:kx_index]-1,ky_index-1,options[:species_index]-1, 0], 'end' => [options[:kx_index]-1,ky_index-1,options[:species_index]-1, -1]})
          #ep phi.shape
          es_heat.reshape(*es_heat.shape.values_at(3))
        end).sum / lky.size
        return es_heat_av.to_gslv
      end

    end
  end
  private :es_heat_flux_by_kxy_over_time_gsl_vector

  def es_heat_flux_over_kx_gsl_vector(options)
    options[:direction] = :kx
    es_heat_flux_over_kxy_gsl_vector(options)
  end

  def es_heat_flux_over_ky_gsl_vector(options)
    options[:direction] = :ky
    es_heat_flux_over_kxy_gsl_vector(options)
  end

  #This function will output the heat flux as a function of kx or ky.
  #Default behaviour will be to average the heat flux over the time domain.
  def es_heat_flux_over_kxy_gsl_vector(options)
    raise "The write_fluxes_by_mode flag was not set." unless @write_fluxes_by_mode
    Dir.chdir(@directory) do
      kxy = options[:direction]
      raise "Please provide species_index " unless options[:species_index]
      if kxy==:ky
        es_heat = (netcdf_file.var('es_heat_flux_by_mode').get({'start' => [0,0,options[:species_index]-1, 0], 'end' => [-1,-1,options[:species_index]-1, -1]})) #index = [kx,ky,spec,t]
        #Need to average over time and sum over kx
        shape = es_heat.shape
        es_heat_av = []; temp = [];
        for iy in 0...shape[1]
          for ix in 0...shape[0]
            temp[ix] = es_heat[ix,iy,0,0..-1].sum / shape[3]
          end
          es_heat_av[iy] = temp.sum
        end
        return es_heat_av.to_gslv
      else
        es_heat = (netcdf_file.var('es_heat_flux_by_mode').get({'start' => [0,0,options[:species_index]-1, 0], 'end' => [-1,-1,options[:species_index]-1, -1]})) #index = [kx,ky,spec,t]
        shape = es_heat.shape
        es_heat_av = []; temp = [];
        for ix in 0...shape[0]
          for iy in 0...shape[1]
            temp[iy] = es_heat[ix,iy,0,0..-1].sum / shape[3]
          end
          es_heat_av[ix] = temp.sum
        end
        return es_heat_av.to_gslv.from_box_order
      end
    end
  end

  def phi2_by_kx_over_time_gsl_vector(options)
    options[:direction] = :kx
    phi2_by_kxy_over_time_gsl_vector(options)
  end
  
  def phi2_by_ky_over_time_gsl_vector(options)
    options[:direction] = :ky
    phi2_by_kxy_over_time_gsl_vector(options)
  end
  
  def phi2_by_kxy_over_time_gsl_vector(options)
    Dir.chdir(@directory) do
      # i.e. phi2_by_ky_vs_time or phi2_by_kx_vs_time

      kxy = options[:direction]
      if list(kxy).size == 1
        return phi2tot_over_time_gsl_vector(options)
      end
      kxy_index = kxy + :_index

      #Necessary options: :ky or :kx
      #Optional options: :t_index_window
      phi_t_array=nil
      if @grid_option == "single"
        phi_t_array = netcdf_file.var('phi2').get('start' => [options[:begin_element]], 'end' => [options[:end_element]]).to_a.flatten
      else
        options.convert_to_index(self, kxy)
        phi_t_array = netcdf_file.var("phi2_by_#{kxy}").get('start' => [options[kxy_index] - 1, options[:begin_element]], 'end' => [options[kxy_index] - 1, options[:end_element]]).to_a.flatten
      end
      return GSL::Vector.alloc(phi_t_array)
    end
  end
  private :phi2_by_kxy_over_time_gsl_vector

  def phi2_by_mode_over_time_gsl_vector(options)
    Dir.chdir(@directory) do      #Necessary options: :ky and :kx
      #Optional options: :t_index_window
      #     eputs "got here"
      #options[:begin_element], options[:end_element] = (options[:t_index_window] ? options[:t_index_window].map{|ind| ind -1} : [0, -1])
      options.setup_time_window
      phi_t_array=nil
      if @grid_option == "single"
        phi_t_array = netcdf_file.var('phi2').get('start' => [options[:begin_element]], 'end' => [options[:end_element]]).to_a.flatten
      else
#         value = options[:ky]
#         eputs value
#         get_list_of(:ky)
#         index = @ky_list.find{|index,val| (val-value).abs < Float::EPSILON}[0]
        options.convert_to_index(self, :kx, :ky)
#         p options
        phi_t_array = netcdf_file.var("phi2_by_mode").get('start' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:begin_element]], 'end' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:end_element]]).to_a.flatten
#         eputs 'phi_t_array.size', phi_t_array.size
      end
      return GSL::Vector.alloc(phi_t_array)

    end
  end

  def tpar2_by_mode_over_time_gsl_vector(options)
    Dir.chdir(@directory) do      #Necessary options: :ky and :kx
      #Optional options: :t_index_window
      #     eputs "got here"
      #options[:begin_element], options[:end_element] = (options[:t_index_window] ? options[:t_index_window].map{|ind| ind -1} : [0, -1])
      options.setup_time_window
      tpar_t_array=nil
      if @grid_option == "single"
        tpar_t_array = netcdf_file.var('tpar2').get('start' => [options[:begin_element]], 'end' => [options[:end_element]]).to_a.flatten
      else
#         value = options[:ky]
#         eputs value
#         get_list_of(:ky)
#         index = @ky_list.find{|index,val| (val-value).abs < Float::EPSILON}[0]
        options.convert_to_index(self, :kx, :ky, :species)
#         p options
        tpar_t_array = netcdf_file.var("tpar2_by_mode").get('start' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:species_index] - 1, options[:begin_element]], 'end' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:species_index] - 1, options[:end_element]]).to_a.flatten
#         eputs 'tpar_t_array.size', tpar_t_array.size
      end
      return GSL::Vector.alloc(tpar_t_array)

    end
  end

  def tperp2_by_mode_over_time_gsl_vector(options)
    Dir.chdir(@directory) do      #Necessary options: :ky and :kx
      #Optional options: :t_index_window
      #     eputs "got here"
      #options[:begin_element], options[:end_element] = (options[:t_index_window] ? options[:t_index_window].map{|ind| ind -1} : [0, -1])
      options.setup_time_window
      tperp_t_array=nil
      if @grid_option == "single"
        tperp_t_array = netcdf_file.var('tperp2').get('start' => [options[:begin_element]], 'end' => [options[:end_element]]).to_a.flatten
      else
#         value = options[:ky]
#         eputs value
#         get_list_of(:ky)
#         index = @ky_list.find{|index,val| (val-value).abs < Float::EPSILON}[0]
        options.convert_to_index(self, :kx, :ky, :species)
#         p options
        tperp_t_array = netcdf_file.var("tperp2_by_mode").get('start' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:species_index] - 1, options[:begin_element]], 'end' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:species_index] - 1, options[:end_element]]).to_a.flatten
#         eputs 'tperp_t_array.size', tperp_t_array.size
      end
      return GSL::Vector.alloc(tperp_t_array)

    end
  end

  def phi0_by_kx_by_ky_over_time_gsl_vector(options)
    Dir.chdir(@directory) do
      options.convert_to_index(self, :kx, :ky)
      phi0_array = netcdf_file.var('phi0').get.to_a.map{|arr| arr[options[:kx_index] - 1][options[:ky_index] - 1][options[:ri]]}
      return GSL::Vector.alloc(phi0_array)
    end
  end

  def linked_kx_elements_gsl_vector(options)
    Dir.chdir(@directory) do
      return GSL::Vector.alloc([0]) if @grid_option == "single" or agk?
      if agk? or (@s_hat_input or @shat).abs < 1.0e-5
        #p 'op1', options
        options.convert_to_index(self, :ky, :kx)
        #p 'op2', options
        #eputs "No Magnetic Shear"

#         begin
#           options.convert_to_index(:kx)
#         rescue
#           raise "Must specify kx or kx_index if no magnetics shear"
#         end
# #         theta0 = (options[:theta0] || 0)
# #         theta0 += jump(options) if @g_exb

        #theta0 = (options[:kx_index])
        #if @g_exb and @g_exb.abs > 0.0
          #theta0 += jump(options)
          #theta0 = theta0%((list(:kx).size-1)/2) if list(:kx).size > 1
        #end

        return GSL::Vector.alloc([options[:kx_index] - 1])
      end

      options.convert_to_index(self, :ky, :kx)
      nkx = netcdf_file.var('kx').dims[0].length
#       p nkx
      stride = @jtwist * (options[:ky_index] -1 )
      #stride = 3
      nlinks = [(nkx / stride).floor, 1].max
      theta0 = options[:kx_index] % @jtwist  #(options[:theta0] || 0)
      #log 'stride', stride, 'nlinks', nlinks, 'theta0', theta0
      #if @g_exb and @jtwist > 1 #and options[:t_index]
#         kx_shift = list(:ky)[options[:ky_index]]  * @g_exb
#         p list(:t)[options[:t_index]], options[:t_index], kx_shift

#         kx_shift *=  list(:t)[(options[:t_index] or list(:t).keys.max)]
#         jump = (kx_shift / list(:kx)[2]).round
        #theta0  += (@jtwist - jump(options) % @jtwist) % @jtwist

#         else
#           jump = 0
      #end
      #ep 'stride', stride, 'nlinks', nlinks, 'theta0', theta0
      #ep GSL::Vector.indgen(nlinks / 2,  nkx + theta0 - nlinks / 2 * stride, stride).connect(GSL::Vector.indgen(nlinks / 2, theta0, stride)).reverse if nlinks > 1
      #return [7,5,3,1,34].to_gslv
      return GSL::Vector.alloc([theta0 % jtwist]) if nlinks ==1
      return GSL::Vector.indgen(nlinks / 2,  nkx + theta0 - nlinks / 2 * stride, stride).connect(GSL::Vector.indgen(nlinks / 2, theta0, stride)).reverse

    end
  end

  def spectrum_over_kpar_gsl_vector(options)
    Dir.chdir(@directory) do
 # , /kpar_spectrum/
      #ep 'zero?', (@s_hat_input||@shat)==0.0
      unless agk? or (@s_hat_input||@shat||0.0).abs<1.0e-5
        phi = gsl_vector_complex('phi_along_field_line', options)
        phi = phi.subvector(0,phi.size-1)
        #i = 0
        #phi = phi.collect{|re,im|
          #i+=1; GSL::Complex.alloc(Math.sin(0.1*i), Math.cos(0.1*i))+
          #GSL::Complex.alloc(Math.sin(0.4*i), Math.cos(0.4*i))

        #}
        ##GraphKit.quick_create([phi.square]).gnuplot
        phi_k = phi.forward
        phi_kr = phi_k.square
        case phi_kr.size%2
        when 0
          spec = phi_kr.subvector((phi_kr.size+2)/2, (phi_kr.size-2)/2).connect(phi_kr.subvector(0, (phi_kr.size+2)/2))
        when 1
          spec = phi_kr.subvector((phi_kr.size + 1)/2, (phi_kr.size-1)/2).connect(phi_kr.subvector(0, (phi_kr.size+1)/2))
        end
        ##spec = phi_kr
        #ep 'spec.class', spec.class
        return spec
      else

        gm = gsl_matrix('spectrum_over_ky_over_kpar', options)
        vec = GSL::Vector.alloc(gm.shape[1])
        vec.set_all(0.0)
        for ky_element in 0...gm.shape[0]
          vec+= gm.row(ky_element)
        end
        vec = vec/gm.shape[0]
        return vec
      end
    end
  end
 
  def kpar_gsl_vector(options)

    Dir.chdir(@directory) do
      if agk? or (@s_hat_input||@shat).abs  < 1.0e-5
        dk = 1
        phi = list(:theta).values
      else
        kxe = gsl_vector('linked_kx_elements', options)
        dk = 1.0/kxe.size
      phi = gsl_vector_complex('phi_along_field_line', options)
      end
      case phi.size%2
      when 0
        kpar = GSL::Vector.indgen(phi.size-1, -((phi.size-3)/2))*dk
      when 1
        kpar = GSL::Vector.indgen(phi.size-1, -((phi.size-2)/2))*dk
      end
      #ep 'kpar', kpar, 'phi.size', phi.size

      #ep 'kpar.class', kpar.class
      return kpar

    end
  end

  def phi_along_field_line_gsl_vector(options)
    Dir.chdir(@directory) do
      complex_phi_vector= gsl_vector_complex('phi_along_field_line', options)
      case options[:imrc]
      when :im
        phi_vector = complex_phi_vector.imag
      when :mag
        #mag = true
        phi_vector = complex_phi_vector.abs2
      when :corr
        thetas = gsl_vector('theta_along_field_line', options)
        min = thetas.abs.to_a.index(thetas.abs.min)
        at_0 = complex_phi_vector[min]
#         ep at_0.class
        phi_vector = (complex_phi_vector * (at_0 / at_0.mag).conj).real
#         gsl_complex('correcting_phase', options)).real
      when :real
        phi_vector = complex_phi_vector.real
      else
        raise CRError.new("options[:imrc] was: #{options[:irmc]}")
      end
      phi_vector *= -1.0 if options[:flip]
      (phi_vector /= phi_vector.abs.max; phi_vector *= (options[:height] || 1.0)) if options[:norm]
      phi_vector = phi_vector.reverse if options[:rev]
      return phi_vector

    end
  end

  def theta_along_field_line_gsl_vector(options)
    Dir.chdir(@directory) do
      case @grid_option
      when "single", "range"
        theta_vector = gsl_vector(:theta)
      when "box"
        #eputs "Start theta_along_field_line"

        kx_elements = gsl_vector('linked_kx_elements', options).to_a
#         ep list(:kx).keys.max
#         ep kx_elements[0], list(:kx)[(kx_elements[0] + 1).to_i]
#         ep kx_elements[-1], list(:kx)[(kx_elements[-1] + 1).to_i]
        thetas = gsl_vector(:theta)
#         ep thetas
        #eputs "End theta_along_field_line"
        return thetas if agk? or (@s_hat_input or @shat).abs < 1.0e-5
        if gryfx?
          theta_list = ((1..kx_elements.size).to_a.map do |i|
            thetas * i
          end)
          thetas = theta_list.inject{|o,n| o.connect(n)}
          thetas -= Math::PI*(kx_elements.size-1)
          return thetas

        end
        theta_list = (kx_elements.map do |element|

          kx = list(:kx)[(element + 1).to_i]
#                      ep element
                     #ep 'kx', kx, 'shat', (@s_hat_input or @shat), 'ky',   list(:ky)[options[:ky_index]]
          thetas - 1.0 / (@s_hat_input or @shat) / list(:ky)[options[:ky_index]] * kx
        end).inject{|old, new| old.connect(new)}
#         thetas = gsl_vector(:theta) - 1.0 / @shat / list(:ky)[options[:ky_index]] * list(:kx)[(kx_elements[0] + 1).to_i] #- Math::PI*(kx_elements.size  - 1)
#         get_list_of(:ky, :t)
#         if @g_exb #and options[:t_index]

          if options[:moving]
            theta_list = theta_list  -  Math::PI * 2.0 * (jump(options) / @jtwist)
          else
#             ep 'jump % jtwist is!!', jump(options) % @jtwist
            theta_list = theta_list - Math::PI * 2.0 / @nx.to_f * ((jump(options) % @jtwist).to_f / @jtwist.to_f)
          end
#           jump = 0
#         end
#         theta_list = thetas.dup #gsl_vector(:theta) - Math::PI*kx_elements.size
#         (kx_elements.size - 1).times do
#           thetas = thetas + Math::PI * 2.0
#           theta_list = theta_list.connect(thetas)
#         end
#         pp theta_list.to_a.values_at(0, theta_list.size - 1)
#         pp theta_list.to_a.max
        theta_vector = theta_list
      end
#       theta_vector = theta_vector.reverse if options[:rev]
      theta_vector *= (@shat) if options[:z]
      return theta_vector

    end
  end

  def phi_for_eab_movie_gsl_vector(options)
    Dir.chdir(@directory) do      #options required are x_index, y_index and tm_index (Time)
      mvf_name = @run_name + '.movie.nc'
      raise CRError.new("cannot find file #{mvf_name}") unless FileTest.exist? mvf_name
      ncf = NumRu::NetCDF.open(mvf_name)
#       p ncf.var('phi_by_xmode').get.to_a[0][0][0]
      return GSL::Vector.alloc(ncf.var('phi_by_xmode').get.to_a[options[:tm_index] - 1].map{|xy_arr| xy_arr[options[:x_index] - 1][options[:y_index] - 1]})

    end
  end

  def hflux_tot_over_time_gsl_vector(options)
    Dir.chdir(@directory) do
      options.setup_time_window
      narr = netcdf_file.var('hflux_tot').get('start' => [options[:begin_element]], 'end' => [options[:end_element]])
      #eputs 'Got narr'
      #ep 'hflux_tot', hflux
      #eputs "fixing norm"
      return fix_heat_flux_norm(GSL::Vector.alloc(narr.to_a), options)
    end
  end
  alias :hflux_tot_gsl_vector :hflux_tot_over_time_gsl_vector
  
  def es_heat_flux_over_time_gsl_vector(options)
    Dir.chdir(@directory) do

      options.setup_time_window
      return GSL::Vector.alloc(netcdf_file.var('es_heat_flux').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
    end
  end
  
  def es_heat_par_over_time_gsl_vector(options)
    Dir.chdir(@directory) do

      options.setup_time_window
      return GSL::Vector.alloc(netcdf_file.var('es_heat_par').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
    end
  end
  alias :es_heat_par_gsl_vector :es_heat_par_over_time_gsl_vector
  
  def es_heat_perp_over_time_gsl_vector(options)
    Dir.chdir(@directory) do

      options.setup_time_window
      return GSL::Vector.alloc(netcdf_file.var('es_heat_perp').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
    end
  end
  alias :es_heat_perp_gsl_vector :es_heat_perp_over_time_gsl_vector
  
  def es_heat_flux_over_time_gsl_vector(options)
    Dir.chdir(@directory) do

      options.setup_time_window
      return GSL::Vector.alloc(netcdf_file.var('es_heat_flux').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
    end
  end
  
  def es_mom_flux_over_time_gsl_vector(options)
    Dir.chdir(@directory) do
      options.setup_time_window
      return GSL::Vector.alloc(netcdf_file.var('es_mom_flux').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
    end
  end
  
  def es_part_flux_over_time_gsl_vector(options)
    Dir.chdir(@directory) do
      options.setup_time_window
      return GSL::Vector.alloc(netcdf_file.var('es_part_flux').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
    end
  end
  
  # Velocity space diagnostics: fraction of dist func in higher
  # pitch angle harmonics
  def lpc_pitch_angle_gsl_vector(options)
    raise "Velocity space lpc diagnostics not found" unless FileTest.exist? "#@directory/#@run_name.lpc"
    lpc = GSL::Vector.filescan("#@directory/#@run_name.lpc")
    return lpc[1]
  end
  
  # Velocity space diagnostics: fraction of dist func in higher
  # energy harmonics
  def lpc_energy_gsl_vector(options)
    raise "Velocity space lpc diagnostics not found" unless FileTest.exist? "#@directory/#@run_name.lpc"
    lpc = GSL::Vector.filescan("#@directory/#@run_name.lpc")
    return lpc[2]
  end
  
  # Velocity space diagnostics: integral error due to
  # pitch angle resolution
  def vres_pitch_angle_gsl_vector(options)
    raise "Velocity space vres diagnostics not found" unless FileTest.exist? "#@directory/#@run_name.vres"
    vres = GSL::Vector.filescan("#@directory/#@run_name.vres")
    return vres[1]
  end
  
  # Velocity space diagnostics: integral error due to
  # energy resolution
  def vres_energy_gsl_vector(options)
    raise "Velocity space vres diagnostics not found" unless FileTest.exist? "#@directory/#@run_name.vres"
    vres = GSL::Vector.filescan("#@directory/#@run_name.vres")
    return vres[2]
  end
  
  def par_mom_flux_over_time_gsl_vector(options)
  Dir.chdir(@directory) do

    options.setup_time_window
    # This is a hack... one day some one will put it in the NetCDF file (haha).
    momlines = `grep parmom #@run_name.out`
    mom = []
    momlines.scan(Regexp.new("#{LongRegexen::FLOAT.to_s}$")) do
      mom.push $~[:float].to_f
    end
    options[:end_element] = (mom.size + options[:end_element]) if options[:end_element] < 0
#       p options
    return GSL::Vector.alloc(mom).subvector(options[:begin_element], options[:end_element] - options[:begin_element] + 1)
  end
  end

  def perp_mom_flux_over_time_gsl_vector(options)

    Dir.chdir(@directory) do
      options.setup_time_window
      # This is a hack... one day some one will put it in the NetCDF file (haha).
      momlines = `grep perpmom #@run_name.out`
      mom = []
      momlines.scan(Regexp.new("#{LongRegexen::FLOAT.to_s}$")) do
        mom.push $~[:float].to_f
      end
      options[:end_element] = (mom.size + options[:end_element]) if options[:end_element] < 0
#       p options
      return GSL::Vector.alloc(mom).subvector(options[:begin_element], options[:end_element] - options[:begin_element] + 1)
    end
  end

  def scan_parameter_value_gsl_vector(options)
    return GSL::Vector.alloc(netcdf_file.var('scan_parameter_value').get.to_a)
  end
  
  def spectrum_over_kx_gsl_vector(options)
    options[:direction] = :kx
    spectrum_over_kxy_gsl_vector(options)
  end

  def spectrum_over_kx_avg_gsl_vector(options)
    options[:direction] = :kx
    spectrum_over_kxy_avg_gsl_vector(options)
  end
 
  def spectrum_over_ky_gsl_vector(options)
    options[:direction] = :ky
    spectrum_over_kxy_gsl_vector(options)
  end
  
  def spectrum_over_ky_avg_gsl_vector(options)
    options[:direction] = :ky
    spectrum_over_kxy_avg_gsl_vector(options)
  end
  
  def spectrum_over_kxy_gsl_vector(options)
    Dir.chdir(@directory) do
      # i.e. spectrum_over_ky or spectrum_over_kx
      kxy = options[:direction]
#       eputs options[:t_index]
      raise "Spectrum makes no sense for single modes" if @grid_option == "single"

      options.convert_to_index(:t) if options[:t] or options[:t_element]
#       eputs options[:t_index]

      options[:t_index] ||= list(:t).keys.max
#       eputs options[:t_index]
      phi_array = netcdf_file.var("phi2_by_#{kxy}").get('start' => [0, options[:t_index] - 1], 'end' => [-1, options[:t_index] - 1]).to_a.flatten
      v = GSL::Vector.alloc(phi_array)
      v = v.from_box_order if kxy == :kx
      v = v.mul(gsl_vector(kxy).square) unless options[:phi2_only]
      return v
    end
  end

  #spectrum averaged in time
  def spectrum_over_kxy_avg_gsl_vector(options)
    Dir.chdir(@directory) do
      # i.e. spectrum_over_ky or spectrum_over_kx
      kxy = options[:direction]
      raise "Spectrum makes no sense for single modes" if @grid_option == "single"

      phi_array = netcdf_file.var("phi2_by_#{kxy}").get('start' => [0, 0], 'end' => [-1, -1]) #index = [kx or ky, t]

      shape = phi_array.shape
      phi_av = [];
      #average over time for each kx or ky individually
      for i in 0...shape[0]
        phi_av[i] = phi_array[i,0..-1].sum / shape[1]
      end

      v = GSL::Vector.alloc(phi_av)
      v = v.from_box_order if kxy == :kx
      v = v.mul(gsl_vector(kxy).square) unless options[:phi2_only]
      return v
    end
  end

  def x_gsl_vector(options)
    raise "options nakx and interpolate_x are incompatible" if options[:nakx] and options[:interpolate_x]
    kx = gsl_vector(:kx, options)
    lx = 2*Math::PI/kx.to_box_order[1]
    #ep 'lx', lx
    nx = options[:nakx]||kx.size
    GSL::Vector.indgen(nx, 0, lx/nx)
  end
  
  def y_gsl_vector(options)
    raise "options naky and interpolate_y are incompatible" if options[:naky] and options[:interpolate_y]
    ky = gsl_vector(:ky, options)
    ly = 2*Math::PI/ky[1]
    ny = options[:naky]||ky.size
    ysize = ny*2-2+ny%2
    GSL::Vector.indgen(ysize, 0, ly/ysize)
  end

  #This function reads in the 'drhodpsi' variable from the netcdf file.
  def drhodpsi_gsl_vector(options)
    drhodpsi = netcdf_file.var('drhodpsi').get()[0]
    return drhodpsi
  end

  #This function returns the zonal flow velocity as a function of x (the radial coordinate).
  #This is v_ZF = kxfac*IFT(i k_x phi_imag), where kxfac = (qinp/rhoc)*grho(rhoc).
  def zf_velocity_over_x_gsl_vector(options)
    Dir.chdir(@directory) do
      raise CRFatal.new("Need either qinp or pk and epsl specified in order to calculate kxfac.
                        If using numerical equil use the option :kxfac to override calculation.") unless @qinp or (@pk and @eps) or options[:kxfac]

      kx = gsl_vector(:kx).to_box_order
      drhodpsi = gsl_vector('drhodpsi')

      phi = GSL::Vector.alloc(kx.size)
      for it in 0...gsl_vector(:t).size
        options[:t_index] = it
        phi += gsl_vector_complex('phi_zonal', options)
      end
      phi /= gsl_vector(:t).size

      if @qinp
        kxfac = (@qinp/@rhoc)/drhodpsi
      elsif @pk and @epsl
        kxfac = (@epsl/@pk)/drhodpsi
      elsif options[:kxfac]
        kxfac = options[:kxfac]
      end

      vec_zf_vel = GSL::Vector.alloc(kx.size)
      #Take imaginary part since i k_x will lead to imaginary part being real
      vec_zf_vel = 0.5*kxfac*(phi*kx).backward.imag
      return vec_zf_vel
    end
  end

  #This function returns the mean flow velocity as a function of x (the radial coordinate).
  #This is v_g_exb = (x - x(centre))*g_exb. The x-x(centre) ensures that the flow is zero
  #at the middle of the box.
  def mean_flow_velocity_over_x_gsl_vector(options)
    Dir.chdir(@directory) do
      raise CRFatal.new("Need to have g_exb > 0 to have a mean flow.") unless @g_exb
      x = gsl_vector(:x)

      vec_exb_vel = GSL::Vector.alloc(x.size)
      #Take imaginary part since i k_x will lead to imaginary part being real
      vec_exb_vel = (x - x[x.size/2])*@g_exb
      return vec_exb_vel
    end
  end

  def zonal_spectrum_gsl_vector(options)
    Dir.chdir(@directory) do
      gmzf = gsl_matrix('spectrum_over_ky_over_kx',options)
      veczf = GSL::Vector.alloc(gmzf.shape[1])
      gmzf.shape[1].times{|i| veczf[i] = gmzf[0,i]}
      return veczf
    end
  end

end # module GSLVectors
include GSLVectors

def gsl_vector_complex(name, options={})
  options = eval(options) if options.class == String

    if method = self.class.instance_methods.find{|meth| (name + '_gsl_vector_complex').to_sym == meth}
      options[:graphkit_name] = name
      return send(method, options)
    end
end

module GSLVectorComplexes

  def phi_along_field_line_gsl_vector_complex(options)
  Dir.chdir(@directory) do
#     eputs options[:ky]
#     eputs Dir.pwd
      #eputs "Start phi_along_field_line"
      options.convert_to_index(self, :ky)
      if options[:t_index] or options[:t]
        #extra option required is t_index
        raise CRFatal.new("write_phi_over_time is not enabled so this function won't work") unless (@write_phi_over_time and @write_phi_over_time.fortran_true?)

        options.convert_to_index(self, :t)
        case @grid_option
        when "single"
          temp = GSL::Vector.alloc(netcdf_file.var('phi_t').get({'start' => [0,0,0,0, options[:t_index] - 1], 'end' => [-1,-1,0,0, options[:t_index] - 1]}).to_a[0][0][0].flatten)
        when "range"
          a = netcdf_file.var('phi_t').get({'start' => [0, 0, options[:kx_index]-1, options[:ky_index] - 1, options[:t_index] - 1], 'end' => [-1, -1, options[:kx_index]-1, options[:ky_index] - 1, options[:t_index]-1]})
          #temp =  GSL::Vector.alloc(a.to_a[0].values_at(*kx_elements).flatten)
          temp =  GSL::Vector.alloc(a.to_a[0][0].flatten)
        when "box"
          options.convert_to_index(self, :ky, :kx)
          kx_elements = gsl_vector('linked_kx_elements', options).to_a
  #         pp kx_elements
          a = netcdf_file.var('phi_t').get({
            'start' => [0,0,0,options[:ky_index] - 1, options[:t_index] - 1],
            'end' => [-1,-1,-1, options[:ky_index] - 1, options[:t_index] - 1]
          }).to_a[0][0].values_at(*kx_elements).flatten
  #         pp a.index(nil)
  #         temp = GSL::Vector.alloc(netcdf_file.var('phi').get.to_a[options[:ky_index] - 1 ].values_at(*kx_elements).flatten)
          #ep a
          temp = GSL::Vector.alloc(a)
        end

        #eputs "End phi_along_field_line"
        return GSL::Vector::Complex.alloc(temp.subvector_with_stride(0, 2), temp.subvector_with_stride(1, 2))
      else
        case @grid_option
        when "single"
          temp = GSL::Vector.alloc(netcdf_file.var('phi').get({'start' => [0,0, 0, 0], 'end' => [-1,-1,0,0]}).to_a.flatten)
        when "range"
          a = netcdf_file.var('phi').get({'start' => [0, 0, 0, options[:ky_index] - 1], 'end' => [-1, -1, -1, options[:ky_index] - 1]})
          #temp =  GSL::Vector.alloc(a.to_a[0].values_at(*kx_elements).flatten)
          temp =  GSL::Vector.alloc(a.to_a[0][0].flatten)
        when "box"
          a = netcdf_file.var('phi').get({'start' => [0, 0, 0, options[:ky_index] - 1], 'end' => [-1, -1, -1, options[:ky_index] - 1]})
          temp =  GSL::Vector.alloc(a.to_a[0].values_at(*kx_elements).flatten)
        else
          raise "invalid grid option"
        end

        vector = GSL::Vector::Complex.alloc(temp.subvector_with_stride(0, 2), temp.subvector_with_stride(1, 2))
        #ep 'vector', vector.real
        return vector
      end
    end
  #       eputs data; gets
  end

  #This function returns a complex GSL vector of the zonal (ky=0) component of phi/phi_t at a given theta index
  #and time index if write_phi_over_time was enabled during the simulation
  def phi_zonal_gsl_vector_complex(options)
    Dir.chdir(@directory) do
        raise CRFatal.new("write_eigenfunc is not enabled so this function won't work") unless @write_eigenfunc
        options.convert_to_index(self, :t)
        a = netcdf_file.var('phi_igomega_by_mode').get({
          'start' => [0,0,0, options[:t_index] - 1],
          'end' => [-1,-1,0, options[:t_index] - 1]
        })
        vector = GSL::Vector::Complex.alloc(GSL::Vector.alloc(a[0,0..-1,0,0]), GSL::Vector.alloc(a[1,0..-1,0,0]))
        return vector
    end
  end

end
include GSLVectorComplexes

def gsl_matrix(name, options={})
  options = eval(options) if options.class == String
  if options[:saturated_time_average] or options[:sta]
    raise "Not Saturated" unless @saturation_time_index
    tmax = list(:t).keys.max
    return ((@saturation_time_index..tmax).to_a.map do |t_index|
      gsl_matrix(name, options.dup.absorb({t_index: t_index, saturated_time_average: nil, sta: nil}))
    end).sum / (list(:t).values.max - list(:t)[@saturation_time_index])
  end
  if method = self.class.instance_methods.find{|meth| (name + '_gsl_matrix').to_sym == meth}
      options[:graphkit_name] = name
      return send(method, options)
  end
end

module GSLMatrices
  
  def growth_rate_over_ky_over_kx_gsl_matrix(options)
    if @growth_rate_at_ky_at_kx.nil?
       raise("The CodeRunner variable growth_rate_at_ky_at_kx does not seem to have been calculated for this run. This may result when the environment variable GS2_CALCULATE_ALL is not set when the run was analyzed. Try setting GS2_CALCULATE_ALL and then re-analyze the run using, e.g. from the command line,\n $ coderunner rc 'cgrf\' -j #{@id}")
    end
    array = @growth_rate_at_ky_at_kx.values.map{|h| h.values}
    return GSL::Matrix.alloc(array.flatten, array.size, array[0].size)
  end
  
  def transient_amplification_over_ky_over_kx_gsl_matrix(options)
      array = @transient_amplification_at_ky_at_kx.values.map{|h| h.values}
      return GSL::Matrix.alloc(array.flatten, array.size, array[0].size)
  end
  
  def es_heat_flux_over_ky_over_kx_gsl_matrix(options)
    raise "The write_fluxes_by_mode flag was not set." unless @write_fluxes_by_mode
    Dir.chdir(@directory) do
      raise "Heat flux spectrum makes no sense for single modes" if @grid_option == "single"
      options.convert_to_index(:t) if options[:t] or options[:t_element]
      options[:t_index] ||= list(:t).keys.max
      #es_heat_by_k index order (in Fortran) is kx, ky, t
      es_heat_narray = netcdf_file.var("es_heat_flux_by_mode").get('start' => [0, 0, 0, options[:t_index] - 1], 'end' => [-1, -1, 0, options[:t_index] - 1])
      es_heat_narray.reshape!(*es_heat_narray.shape.slice(0..1))

      gm =  es_heat_narray.to_gm.move_cols_from_box_order
      if options[:limit]
        for i in 0...gm.shape[0]
          for j in 0...gm.shape[1]
#             j+= extra if
            gm[i, j] = [[gm[i,j], (options[:limit][0] or gm[i,j])].max, (options[:limit][1] or gm[i,j])].min
#             mat[i, j+extra] = gm[i,-j] unless j==0
          end
        end
      end
      return gm
  end
  end
  
  def spectrum_over_ky_over_kx_gsl_matrix(options)
  Dir.chdir(@directory) do
      raise "Spectrum makes no sense for single modes" if @grid_option == "single"
      options.convert_to_index(:t) if options[:t] or options[:t_element]
      options[:t_index] ||= list(:t).keys.max
      #phi2_by_mode index order (in Fortran) is kx, ky, t
      phi_narray = netcdf_file.var("phi2_by_mode").get('start' => [0, 0, options[:t_index] - 1], 'end' => [-1, -1, options[:t_index] - 1])
      phi_narray.reshape!(*phi_narray.shape.slice(0..1))

      gm =  phi_narray.to_gm.move_cols_from_box_order
      if options[:times_kx4] or options[:times_kx2]
#         puts 'normalising'
        vals = list(:kx).values.sort
        for i in 0...gm.shape[0]
          for j in 0...gm.shape[1]
#             p vals[j]
            gm[i,j] =  gm[i,j] * (vals[j])**4 if options[:times_kx4]
            gm[i,j] =  gm[i,j] * (vals[j])**2 if options[:times_kx2]
          end
        end
      end
      if options[:no_zonal]

        for i in 0...gm.shape[1]
          gm[0,i] = 0.0
        end
      end
      if options[:log]
        gm = gm.log
      end

      return gm
  end
  end
  
  def spectrum_over_ky_over_kpar_gsl_matrix(options)
  Dir.chdir(@directory) do

      #:re, :theta, :kx, :ky
      lkx = list(:kx)

      if options[:t_index] or options[:t]
        #extra option required is t_index
        raise CRFatal.new("write_phi_over_time is not enabled so this function won't work") unless @write_phi_over_time
        options.convert_to_index(self, :t)
      end
       phi_av = (lkx.keys.map do |kx_index|
        if options[:t_index]
          phi =  netcdf_file.var('phi_t').get({'start' => [0,0,kx_index-1,0, options[:t_index] - 1], 'end' => [-1,-2,kx_index-1,-1, options[:t_index] - 1]})
        else
          phi = netcdf_file.var('phi').get({'start' => [0, 0, kx_index - 1, 0], 'end' => [-1, -2, kx_index-1, -1]})
        end
        #ep phi.shape
        phi.reshape(*phi.shape.values_at(0,1,3))
      end).sum / lkx.size

      phi_t = phi_av.to_a #.map{|arr| arr.transpose}.transpose.map{|a| a.transpose}
      #ep 'phi_t', phi_t.size, phi_t[0].size, phi_t[0][0].size
      gvky = gsl_vector('ky')
      gm = GSL::Matrix.alloc(gvky.size, gsl_vector('theta').size-1)
      for ky_element in 0...gm.shape[0]
        #p phi_t[ky_element].transpose[0]
        spectrum = GSL::Vector::Complex.alloc(phi_t[ky_element]).forward.square
        if options[:no_kpar0]
          spectrum[0]=0.0
        end
        #ep spectrum.size
        spectrum = spectrum.from_box_order
        #ep spectrum.shape
        spectrum = spectrum*gvky[ky_element]**2 unless options[:phi2_only]
        #ep gm.size
        #ep spectrum.size
        gm.set_row(ky_element, spectrum)
      end
      if options[:no_zonal]
        gm.row(0).set_all(0.0)
      end
      if options[:log]
        gm = gm.log
      end

      return gm
  end
  end

  def phi0_over_x_over_y_gsl_matrix(options)
  Dir.chdir(@directory) do

      #options.convert_to_index(:t) if options[:t] or options[:t_element]
      options.convert_to_index(self, :t) if options[:t] or options[:t_element]
      options[:t_index] ||= list(:t).keys.max
      phi_re_narray = netcdf_file.var("phi0").get('start' => [0, 0, 0, options[:t_index] - 1], 'end' => [0, -1, -1, options[:t_index] - 1])
      phi_re_narray.reshape!(*phi_re_narray.shape.slice(1..2))
      # The narray has index order ky, kx, but we want kx, ky for historical reasons, hence the transpose.
      gm_re = phi_re_narray.to_gm
      phi_im_narray = netcdf_file.var("phi0").get('start' => [1, 0, 0, options[:t_index] - 1], 'end' => [1, -1, -1, options[:t_index] - 1])
      phi_im_narray.reshape!(*phi_im_narray.shape.slice(1..2))
      # The narray has index order ky, kx, but we want kx, ky for historical imasons, hence the transpose.
      gm_im = phi_im_narray.to_gm
      gm = GSL::Matrix::Complex.re_im(gm_re, gm_im)

      ntheta0_temp = gm.shape[1]
      naky_temp = gm.shape[0]



      # Due to a strange GS2 convention, non zonal modes must be divided by 2:
      for i in 1...naky_temp
        for j in 0...ntheta0_temp
          gm[i,j] = gm[i,j]/2.0
        end
      end

      if options[:no_zonal]
        for i in 0...gm.shape[1]
          gm[0,i] = GSL::Complex.alloc([0,0])
        end
      end
      if xres = (options[:xres] or options[:x_resolution])
         if xres < nx
            puts "Warning: xres should be at least nx. Using nx instead of the xres you specified."
            xres = nx
         end
      else
         xres = nx
      end
      if yres = (options[:yres] or options[:y_resolution])
         if yres < ny
            puts "Warning: yres should be at least ny. Using ny instead of the yres you specified."
            yres = ny
         end
      else
         yres = ny
      end

      # Next, pad with 0's:
      padded = GSL::Matrix::Complex.calloc(yres, xres)
      # Zonal modes first:
      for ix in 0...((ntheta0_temp+1)/2)
          padded[0, ix] = gm[0, ix]
      end
      for ix in ((ntheta0_temp+1)/2)...ntheta0_temp
          padded[0, ix+xres-ntheta0_temp] = gm[0, ix]
      end
      # Now include the non-zonal modes in the padded matrix mat:
      for iy in 1...naky_temp
          for ix in 0...((ntheta0_temp + 1)/2)
             padded[iy, ix] = gm[iy, ix]
          end
          for ix in ((ntheta0_temp+1)/2)...ntheta0_temp
             padded[iy, ix+xres-ntheta0_temp] = gm[iy, ix]
          end
          padded[yres-iy, 0] = gm[iy,0].conj
          for ix in 1...xres
             padded[yres-iy, ix] = padded[iy, xres-ix].conj
          end
      end
      gm = padded

      gm = gm.backward_cols_c2c(false).backward_rows_c2c(false)
      # At this point, gm should be purely real (within machine precision), but let's check to be sure:
      should_be_zero = gm.imag.abs.max
      if should_be_zero > 1.0e-10
         puts "should_be_zero = #{should_be_zero}"
         raise "Something went wrong - reconstructed phi is not purely real."
      end

      gm = gm.real

      if options[:limit]
        for i in 0...gm.shape[0]
          for j in 0...gm.shape[1]
            gm[i, j] = [[gm[i,j], options[:limit][0]].max, options[:limit][1]].min
          end
        end
      end

      return gm
  end
  end
end

include GSLMatrices

def kx_shift(options)
# ep options
  return 0 unless @g_exb and @g_exb.abs > 0.0
  #p options
  return - list(:ky)[options[:ky_index]] * list(:t)[(options[:t_index] or list(:t).keys.max)] * @g_exb
end

def jump(options)
# ep 'kx_shift',  kx_shift(options)
  jump =  ((kx_shift(options) / list(:kx)[2]).round)
  case options[:t_index]
  when 1
    return jump
  else
    if @g_exb and @g_exb.abs > 0
      return jump + 1
    else
      return 0
    end
  end
end

# This function is used in the presence of perpendicular flow shear. It returns the (Eulerian) GS2
# kx_index as a function of the Lagrangian kx,  which is the kx_index of the mode in a shearing
# coordinate system, I.e. if you give it an Lagrangian kx (which is the same as the Eulerian
# kx at t=0) it will tell you where it has now got to. It may have left the box, in which case
# this function will return an error.
#
# A given Lagrangian kx moves through the GS2 box, and thus for such a kx the response matrix varies
# in time. This is done because the effect of flow shear can be reduced by a shearing coordinate
# transformation to become merely a time varying kx.
#
# At each timestep, phi(ikx_indexed(it)) is set equal to phi(ikx_indexed(it - jump(iky))
# kx_indexed is defined in the following way.
#   do it=itmin(1), ntheta0
#     ikx_indexed (it+1-itmin(1)) = it
#   end do
#
#   do it=1,itmin(1)-1
#     ikx_indexed (ntheta0 - itmin(1) + 1 + it)= it
#   end do
#
# In other words, what this means is that akx(ikx_indexed(0)) is the minimum kx,
# and that akx(ikx_indexed(ntheta0)) gives the maximum kx, kx_indexed moves the
# kxs out of box order.
#
# So. remembering that jump is negative, phi(kx) is set equal phi(kx - jump * dkx)
# so the Lagrangian mode has moved to a lower kx. So get the Eulerian index, one
# starts with the Lagrangian index, and adds jump (which is negative!). This, however,
# must be done with indexes that are in the physical (not box) order. So this function
# first moves the indexes out of box order, then adds jump, then moves them back
# into box order so that the index returned will give the correct kx from the GS2
# array.

def eulerian_kx_index(options)
  #eputs "Start eulerian_kx_index"
  lagrangian_kx_index = options[:kx_index]
  phys = physical_kx_index(lagrangian_kx_index)
  #ep 'jump', jump(options)
  index = phys + jump(options)
  raise ArgumentError.new("Lagrangian kx out of range") if index <= 0
  box= box_kx_index(index)
  #eputs "End eulerian_kx_index"
  return box
end

def kx_indexed
  return cache[:kx_indexed] if cache[:kx_indexed]
  #kx = cache[:kx_array] ||= gsl_vector('kx').to_a
  #kxphys = kx.from_box_order
  #min_index = kx.min_index + 1
  #cache[:kx_indexed] ||= kx.size.times.inject({}) do |hash, kx_element|
    #hash[kx_element + 1] = kxphs
  kx = gsl_vector('kx')
  size = kx.size
  box =  GSL::Vector::Int.indgen(size) + 1
  zero_element = kx.abs.min_index
  phys = box.subvector(zero_element, size-zero_element).connect(box.subvector(0, zero_element))
  cache[:kx_indexed] = [phys.to_a, box.to_a].transpose.inject({}){|hash, (physi, boxi)| hash[physi] = boxi; hash}
end

def box_kx_index(physical_kx_index)

  return kx_indexed[physical_kx_index]
end

def physical_kx_index(box_kx_index)
  return kx_indexed.key(box_kx_index)
  #kx = cache[:kx_gslv] ||= gsl_vector('kx')
  #return kx.from_box_order.to_a.index(kx[box_kx_index-1]) + 1
  #kx = cache[:kx_gslv] ||= gsl_vector('kx')
  #index_of_min_kx = cache[:index_of_min_kx] ||= kx.min_index + 1 # kx.min_index returns a 0-based index
  #if box_kx_index < index_of_min_kx
    #box_kx_index + (1 + kx.size - index_of_min_kx)
  #else
    #box_kx_index - (index_of_min_kx - 1)
  #end
end

def gsl_complex(name, options={})
  options = eval(options) if options.class == String
#   p @directory
  Dir.chdir(@directory) do
#     eputs Dir.pwd
    case name
    when /correcting_phase/
#       options.convert_to_index(self, :ky)
#       theta0 = (options[:theta0] or 0)
# #       p 'options[:ky_index]', options[:ky_index]
#       phase_array = NumRu::NetCDF.open("#@directory/#@run_name.out.nc").var('phase').get({"start" => [0, options[:ky_index] - 1, theta0], 'end' => [1, options[:ky_index] - 1, theta0] }).to_a.flatten
#       p 'phase_array', phase_array
#       thetaelement0 = (list(:theta).key(0.0) - 1).to_i
#       p 'list(:theta)[thetaelement0 + 1]', list(:theta)[thetaelement0 + 1]
#       p 'thetaelement0', thetaelement0
#       p 'theta0 - jump(options)', theta0 - jump(options) % @jtwist
#       p 'list(:kx)[2] * (theta0 - jump(options)%@jtwist)', list(:kx)[2] * (theta0 - jump(options)%@jtwist)
#       kx_element = list(:kx).key(list(:kx)[2] * (theta0 - jump(options)%@jtwist)) - 1
#       at_0 = NumRu::NetCDF.open("#@directory/#@run_name.out.nc").var('phi').get({"start" => [0, thetaelement0, kx_element, options[:ky_index] - 1], 'end' => [1, thetaelement0, kx_element, options[:ky_index] - 1] }).to_a.flatten
#       p 'at_0', at_0
#       at_0 = GSL::Complex.alloc(at_0)
#       p 'at_0', at_0
#       return (at_0 / at_0.mag).conj
# #       pp 'theta0', theta0
# #       pp phase_array[5][theta0]
#       return GSL::Complex.alloc(phase_array)
# #       new_options = options.dup
#       new_options[:imrc] = :real
#       thetas = gsl_vector('theta_along_field_line', new_options)
#       at_0 = gsl_vector_complex('phi_along_field_line', new_options)[.to_a.index(0.0)]
#       p at_0
      exit
    else
      raise CRError.new("Unknown gsl_complex requested: #{name}")
    end
  #       eputs data; gets
  end
end

# def gsl_matrix(name, options={})
#   if options[:t_index] or options[:frame_index]
#     return get_gsl_matrix(name, options)
#   else
#     return cache[[:gsl_vector, name, options]] ||= get_gsl_matrix(name, options)
#   end
# end

end
end