#!/usr/bin/ruby

require "narray"
require "numru/dcl"
include NumRu
include Math


def wl2rgb( wl )
  # Unit of wl is nm.

  # original code is http://www.physics.sfasu.edu/astro/color/spectra.html

  gamma = 0.8
  nwl = wl.shape[0]
  rgb = NArray.sfloat(nwl,3)

  for i in 0..nwl-1
    if (wl[i] < 380) then
      r = -1.0
      g = -1.0
      b = -1.0
    end
    if (wl[i] >= 380) && (wl[i] <= 440) then
      r = -1.0*(wl[i]-440)/(440-380)
      g = 0.0
      b = 1.0
    end
    if (wl[i] >= 440) && (wl[i] <= 490) then
      r = 0.0
      g = (wl[i]-440)/(490-440)
      b = 1.0
    end
    if (wl[i] >= 490) && (wl[i] <= 510) then
      r = 0.0
      g = 1.0
      b = -1.0*(wl[i]-510)/(510-490)
    end
    if (wl[i] >= 510) && (wl[i] <= 580) then
      r = (wl[i]-510)/(580-510)
      g = 1.0
      b = 0.0
    end
    if (wl[i] >= 580) && (wl[i] <= 645) then
      r = 1.0
      g = -1.0*(wl[i]-645)/(645-580)
      b = 0.0
    end
    if (wl[i] >= 645) && (wl[i] <= 780) then
      r = 1.0
      g = 0.0
      b = 0.0
    end

    #      LET THE INTENSITY SSS FALL OFF NEAR THE VISION LIMITS
    if wl[i] > 700 then
      sss=0.3+0.7* (780-wl[i])/(780-700)
    elsif wl[i] < 420 then
      sss=0.3+0.7*(wl[i]-380)/(420-380)
    else
      sss=1.0
    end

    #      GAMMA ADJUST AND WRITE IMAGE TO AN ARRAY
    if (wl[i] >= 380) && (wl[i] <= 780) then
      rgb[i,0]=(sss*r)**gamma
      rgb[i,1]=(sss*g)**gamma
      rgb[i,2]=(sss*b)**gamma
    else
      rgb[i,0]=-1.0
      rgb[i,1]=-1.0
      rgb[i,2]=-1.0
    end
  end

  rgb = rgb * 255 / rgb.max

end

def wl2rgb_old( wl )

  # see http://www.techmind.org/colour/spectra.html

  rgb = NArray.sfloat(wl.shape[0],3)

  #for i in 0..nx-1
  #  # Simple , Single-lobe Fit
  #  # CIE 1931
  #  # see http://jcgt.org/published/0002/02/01/paper.pdf
  #  x = 1.065 * exp(-0.5*((wl-595.8)/33.33)**2) + 0.366 * exp(-0.5*((wl-446.8)/19.44)**2)
  #  y = 1.014 * exp(-0.5*((Math::log(wl/556.3))/0.075)**2)
  #  z = 1.839 * exp(-0.5*((Math::log(wl/449.8))/0.051)**2)
  #  r =   2.3706743 * x - 0.9000405 * y - 0.4706338 * z
  #  g = - 0.5138850 * x + 1.4253036 * y + 0.0885814 * z
  #  b =   0.0052982 * x - 0.0146949 * y + 1.0093968 * z
  #  rgb[i,0] = r ; rgb[i,1] = g ; rgb[i,2] = b
  ##  print wl, " ", y, "\n"
  ##  print wl, " ", x, " ", y, " ", z, "\n"
  ##  print wl, " ", r, " ", g, " ", b, "\n"
  ##  DCL::sgtnxu(vpx,vpy,999,DCL::isgrgb(255.0/(nx-1)*i,0,0))
  #end

  # Simple , Single-lobe Fit
  # CIE 1931
  # see http://jcgt.org/published/0002/02/01/paper.pdf
  x = 1.065 * NMath::exp(-0.5*((wl-595.8)/33.33)**2) + 0.366 * NMath::exp(-0.5*((wl-446.8)/19.44)**2)
  y = 1.014 * NMath::exp(-0.5*((NMath::log(wl/556.3))/0.075)**2)
  z = 1.839 * NMath::exp(-0.5*((NMath::log(wl/449.8))/0.051)**2)
  r =   2.3706743 * x - 0.9000405 * y - 0.4706338 * z
  g = - 0.5138850 * x + 1.4253036 * y + 0.0885814 * z
  b =   0.0052982 * x - 0.0146949 * y + 1.0093968 * z
  rgb[true,0] = r ; rgb[true,1] = g ; rgb[true,2] = b


  # offset in order to avoid negative green and blue (red has negative value)
  rgb = rgb + ( -rgb[true,1..2].min )

  # index for negative rgb
  inegative = (rgb.lt 0.0).where
  ipositive = (rgb.gt 0.0).where

  rgbnegative = NArray.sfloat(rgb.shape[0],3)
  rgbnegative[] = rgb[]
  rgbnegative[ipositive] = 0.0

  rgbpositive = NArray.sfloat(rgb.shape[0],3)
  rgbpositive[] = rgb[]
  rgbpositive[inegative] = 0.0

  # cliping negative values
#  rgb = rgbpositive

  # 
  rgb = rgbpositive - ( rgbnegative[true,0] + rgbnegative[true,1] + rgbnegative[true,2] )


  rgb = rgb * 255 / rgb.max

  return rgb
end


def wlnm2rgb( wl )

  rgb = NArray.sfloat(wl.shape,3)

  #for i in 0..nx-1
  #  # Simple , Single-lobe Fit
  #  # CIE 1931
  #  # see http://jcgt.org/published/0002/02/01/paper.pdf
  #  x = 1.065 * exp(-0.5*((wl-595.8)/33.33)**2) + 0.366 * exp(-0.5*((wl-446.8)/19.44)**2)
  #  y = 1.014 * exp(-0.5*((Math::log(wl/556.3))/0.075)**2)
  #  z = 1.839 * exp(-0.5*((Math::log(wl/449.8))/0.051)**2)
  #  r =   2.3706743 * x - 0.9000405 * y - 0.4706338 * z
  #  g = - 0.5138850 * x + 1.4253036 * y + 0.0885814 * z
  #  b =   0.0052982 * x - 0.0146949 * y + 1.0093968 * z
  #  rgb[i,0] = r ; rgb[i,1] = g ; rgb[i,2] = b
  ##  print wl, " ", y, "\n"
  ##  print wl, " ", x, " ", y, " ", z, "\n"
  ##  print wl, " ", r, " ", g, " ", b, "\n"
  ##  DCL::sgtnxu(vpx,vpy,999,DCL::isgrgb(255.0/(nx-1)*i,0,0))
  #end

  # Simple , Single-lobe Fit
  # CIE 1931
  # see http://jcgt.org/published/0002/02/01/paper.pdf
  x = 1.065 * NMath::exp(-0.5*((wl-595.8)/33.33)**2) + 0.366 * NMath::exp(-0.5*((wl-446.8)/19.44)**2)
  y = 1.014 * NMath::exp(-0.5*((NMath::log(wl/556.3))/0.075)**2)
  z = 1.839 * NMath::exp(-0.5*((NMath::log(wl/449.8))/0.051)**2)
  r =   2.3706743 * x - 0.9000405 * y - 0.4706338 * z
  g = - 0.5138850 * x + 1.4253036 * y + 0.0885814 * z
  b =   0.0052982 * x - 0.0146949 * y + 1.0093968 * z
  rgb[true,0] = r ; rgb[true,1] = g ; rgb[true,2] = b


  # offset in order to avoid negative green and blue (red has negative value)
  rgb = rgb + ( -rgb[true,1..2].min )

  # index for negative rgb
  inegative = (rgb.lt 0.0).where
  ipositive = (rgb.gt 0.0).where

  rgbnegative = NArray.sfloat(rgb.shape[0],3)
  rgbnegative[] = rgb[]
  rgbnegative[ipositive] = 0.0

  rgbpositive = NArray.sfloat(rgb.shape[0],3)
  rgbpositive[] = rgb[]
  rgbpositive[inegative] = 0.0

  # cliping negative values
  rgb = rgbpositive

  # 
  #rgb = rgbpositive - ( rgbnegative[true,0] + rgbnegative[true,1] + rgbnegative[true,2] )


  rgb = rgb * 255 / rgb.max

  return rgb
end