require "numru/ggraph"
include NumRu

def axis_conv_wn( gp, flag_fluxspe )
  gp = gp.copy

  if ( flag_fluxspe ) then
    # conversion from Wm-2 (m-1)-1 to Wm-2 (cm-1)-1
    gp = gp * 1e-2
#    gp.units = 'Wm|-2"(cm|-1")|-1"'
    gp.units = 'W m-2 cm'
  end

  # conversion from m-1 to cm-1
  wn = gp.coord("WaveNum")
  wn = wn.convert_units( Units['cm-1'] )
  # set axis
  gp.axis("WaveNum").set_pos(wn)

  return gp
end

def axis_conv_wl( gp, flag_fluxspe )
  gp = gp.copy

  wn = gp.coord("WaveNum")  # wavenumber
  wl = 1.0/wn               # wavelength

  if ( flag_fluxspe ) then
    # conversion from Wm-2 (m-1)-1 to Wm-2 (um)-1

#    gp = gp / wl**2 * 1e-6
#    gp *= 1.0 / wl.reshape(1,wl.val.size)**2 * 1e-6
    case gp.rank
    when 1 then
      tmp_wn = wn
    when 2 then
      if ( gp.coord(0).name == 'WaveNum' ) then
        tmp_wn = wn.reshape(wn.val.size,1)
      else
        tmp_wn = wn.reshape(1,wn.val.size)
      end
    else
      print 'Unexpected rank of gp.', "\n"
      exit
    end
    gp *= tmp_wn**2 * 1e-6

#    gp.units = 'Wm|-2"(um)|-1"'
    gp.units = 'W m-2 um-1'
#    gp.units = 'Wm|-2"('+DCL::csgi(163)+'m)|-1"'
#    gp.units = 'W m-2 '+DCL::csgi(163)+'m-1'
  end

  # conversion from m-1 to um
  wl = wl.convert_units( Units['um'] )
  # set axis
  gp.axis("WaveNum").set_pos(wl)
  wl.long_name = "wavelength"

  return gp
end

def wn_running_mean( gp, nave )
  # gp   : input GPhys object
  # nave : number of data for averaging

  # pickup values of wavenumber axis
  na_wn = gp.coord('WaveNum').val
  # calculate size of the new wavenumber axis for an output GPhys object
  nwn = na_wn.length / nave
  # create the output GPhys object
  gpout = gp.cut('WaveNum'=>na_wn[0]..na_wn[nwn-1]).copy
  # pickup axis object of the output GPhys object
  ax_wn = gpout.coord("WaveNum").copy
  # calculate running averaged values
  for i in 0..(nwn-1)
    iwns = i*nave          # first index for each averaging
    iwne = iwns+nave-1     # last  index for each averaging
    # calculate averaged wavenumber
    wn_ave = 0
    for j in 0..(nave-1)
      wn_ave += na_wn[iwns+j]
    end
    wn_ave /= nave
    ax_wn[i] = wn_ave
    # calculate averaged values
    gpout[true,i] = gp.cut('WaveNum'=>na_wn[iwns]..na_wn[iwne]).mean('WaveNum')
  end
  # set new wavenumber axis into the output GPhys object
  gpout.axis("WaveNum").set_pos(ax_wn)
  return gpout
end