pro get_color,n_col,spec_dlvs,chro_file,white_file,rel_waves,rel_resp,wave_cal,$
  blue_ch,green_ch,red_ch,img_ch,spec_rwc

; v2 1/11/0 implements the XYZ chromaticity formalism using the 1931 CIE x,y,z weighting functions
; v3 1/11/0 employs artificial XYZ weighting functions 

del_red=6.0 ;nm
del_green=6.0 ;nm
del_blue=3.0 ;nm

lc_red=600.0 ;nm
lc_green=550.0 ;nm
lc_blue=500.0 ;nm

;white_factor=1.00
white_factor=1.75

img_rel_resp=rel_resp/total(rel_resp)

;mi=[[1.73009,-0.482301,-0.261062],[-0.813527,1.65171,-0.0233881],[0.0834387,-0.169406,1.28445]]
;mi=[[3.247,-1.540,-0.499],[-0.972,1.875,-0.042],[0.057,-0.205,1.060]]
mi=[[1.000,0.000,0.000],[0.000,1.000,0.000],[0.000,0.000,1.000]]


s_p=size(spec_dlvs)
spec_rwc=fltarr(s_p(1),s_p(2),s_p(3))

green_ch=fltarr(s_p(1),n_col)
blue_ch=fltarr(s_p(1),n_col)
red_ch=fltarr(s_p(1),n_col)
img_ch=fltarr(s_p(1),n_col)

case n_col of
  10: begin
        spec_wv=fltarr(10,200)
        for j=0,9 do begin
          spec_wv(j,*)=0.5*(wave_cal(2*j,*)+wave_cal(2*j+1,*))
        endfor
      end
else: begin
        print,'n_col = 20'
        spec_wv=wave_cal
      end
endcase

; Input the response to a white card
wvc=fltarr(200)
wc=fltarr(200)
rwc=fltarr(200)
openr,olw1,white_file,/get_lun
for i=0,199 do begin
  readf,olw1,duma,dumb,dumc
  wvc(i)=duma
  wc(i)=dumb
  rwc(i)=dumc
endfor
wc=wc*white_factor
wc(195:199)=0.0
spec_dlvs(*,*,195:199)=0.0
blue_correct,spec_wv,spec_dlvs

; Input the equal-energy chromaticity standards
chr_wv=fltarr(81)
x_red=fltarr(81)
y_green=fltarr(81)
z_blue=fltarr(81)
openr,olo2,chro_file,/get_lun
for i=0,80 do begin
  readf,olo2,duma,dumb,dumc,dumd
  chr_wv(i)=duma
  x_red(i)=dumb
  y_green(i)=dumc
  z_blue(i)=dumd
endfor

; Create the artificial weighting functions x_red,y_green and z_blue
chr_wv_red=fltarr(25)
x_red=fltarr(25)
for i=0,24 do begin 
  chr_wv_red(i)=lc_red+0.150*del_red*(i-12)
  x_red(i)=1/del_red*sqrt(2./!pi)*exp(-2.*((chr_wv_red(i)-lc_red)/del_red)^2)
endfor
x_red=x_red/total(x_red)

chr_wv_green=fltarr(25)
y_green=fltarr(25)
for i=0,24 do begin 
  chr_wv_green(i)=lc_green+0.150*del_green*(i-12)
  y_green(i)=1/del_green*sqrt(2./!pi)*exp(-2.*((chr_wv_green(i)-lc_green)/del_green)^2)
endfor
y_green=y_green/total(y_green)

chr_wv_blue=fltarr(25)
z_blue=fltarr(25)
for i=0,24 do begin 
  chr_wv_blue(i)=lc_blue+0.150*del_blue*(i-12)
  z_blue(i)=1/del_blue*sqrt(2./!pi)*exp(-2.*((chr_wv_blue(i)-lc_blue)/del_blue)^2)
endfor
z_blue=z_blue/total(z_blue)

wht_blue=0.0
wht_red=0.0
wht_green=0.0
wht_ir=0.0

for id=0,24 do begin  
  wht_blue=wht_blue+z_blue(id)*interpol(wc,wvc,chr_wv_blue(id))
  wht_green=wht_green+y_green(id)*interpol(wc,wvc,chr_wv_green(id))
  wht_red=wht_red+x_red(id)*interpol(wc,wvc,chr_wv_red(id))
endfor
for ii=0,48 do begin
  wht_ir=wht_ir+img_rel_resp(id)*interpol(wc,wvc,rel_waves(id))
endfor

g_tot=0
gw_tot=0
go_tot=0

print,'white RGB'
print,wht_red,wht_green,wht_blue
;print,'RGB'

for i=0,s_p(1)-1 do begin
  for j=0,n_col-1 do begin

    num_red=0.0
    num_green=0.0
    num_blue=0.0
    num_ir=0.0
    
    for iwv=0,24 do begin
      signal_red=interpol(spec_dlvs(i,j,*),spec_wv(j,*),chr_wv_red(iwv))
      signal_green=interpol(spec_dlvs(i,j,*),spec_wv(j,*),chr_wv_green(iwv))
      signal_blue=interpol(spec_dlvs(i,j,*),spec_wv(j,*),chr_wv_blue(iwv))
      num_red=num_red+x_red(iwv)*signal_red
      num_green=num_green+y_green(iwv)*signal_green
      num_blue=num_blue+z_blue(iwv)*signal_blue
    endfor

    for jwv=0,48 do begin
      signal_ir=interpol(spec_dlvs(i,j,*),spec_wv(j,*),rel_waves(jwv))
      num_ir=num_ir+img_rel_resp(jwv)*signal_ir
    endfor

;   print,i,j,num_red,num_green,num_blue,format='(i4,i4,f10.3,f10.3,f10.3)'
    
    arg_rgb=255.0*(mi##[num_red,num_green,num_blue])/(mi##[wht_red,wht_green,wht_blue])+0.5

    gv_index=where(arg_rgb gt 255.0,gv_cnt)
    if gv_cnt ne 0 then arg_rgb(gv_index)=255.0
    go_index=where(arg_rgb lt 0.0,go_cnt)
    if go_cnt ne 0 then arg_rgb(go_index)=0.0
    print,i,j,num_red/wht_red*255,num_green/wht_green*255,num_blue/wht_blue*255,$
      format='(i4,i4,f10.3,f10.3,f10.3)'
    print,arg_rgb(0),arg_rgb(1),arg_rgb(2),format='(f10.3,f10.3,f10.3)'
    rgb=byte(arg_rgb)

    arg_ir=255.0*num_ir/wht_ir
    if arg_ir(0) gt 255.0 then arg_ir(0)=255.0
    if arg_ir(0) lt 1.0 then arg_ir(0)=1.0

    g_tot=g_tot+gv_cnt

    red_ch(i,j)=rgb(0)
    green_ch(i,j)=rgb(1)
    blue_ch(i,j)=rgb(2)
    img_ch(i,j)=arg_ir

    arg_rgb_wht=255.0*(mi##[wht_red,wht_green,wht_blue])/(mi##[wht_red,wht_green,wht_blue])
    gw_index=where(arg_rgb_wht gt 255.0,gw_cnt)
    if gw_cnt ne 0 then arg_rgb_wht(gw_index)=255.0
    rgb_wht=byte(arg_rgb_wht)    
    gw_tot=gw_tot+gw_cnt
    go_tot=go_tot+go_cnt

    spec_rwc(i,j,*)=spec_dlvs(i,j,*)/rwc(*)
;   if i eq 20 and j eq 3 then stop
;   if i eq 20 and j eq 10 then stop
;   if i eq 20 and j eq 14 then stop
  endfor
  print,i,'   Max = ',max(spec_dlvs(i,*,*))
endfor

print,'g_tot = ',g_tot
print,'gw_tot = ',gw_tot
print,'go_tot = ',go_tot

stop
return
end