function  deltae2k, lab2, lab1, k

; delta E 2000 code created to match the spreadsheet from Klaus Witt
; DWyble 06Aug02

; If you have some time, make this work for 3xN lab vectors
; As it stands, you have to pass in two 3x1 lab triplets and an optional triplet of k's

    if (n_params() lt 3) then begin
		k = [1, 1, 1]
    endif

    meanL = (lab1(0)+lab2(0))/2.
    
    Cab1 = (lab1(1)^2.+lab1(2)^2.)^.5
    Cab2 = (lab2(1)^2.+lab2(2)^2.)^.5
    meanC = (Cab1 + Cab2)/2.
    
    G = 0.5*(1-((meanC^7.)/(meanC^7. + 25^7.))^0.5)
        
    aprime1=lab1(1)*(1+G)
    aprime2=lab2(1)*(1+G)
    
    CabPrime1 = (aprime1^2.+lab1(2)^2.)^.5
    CabPrime2 = (aprime2^2.+lab2(2)^2.)^.5
    meanCprime = (CabPrime1+CabPrime2)/2
    
    hprime1 = atan(lab1(2), aprime1)*180/!pi
    hprime2 = atan(lab2(2), aprime2)*180/!pi
    if hprime1 LT 0 then hprime1 = hprime1 + 360
    if hprime2 LT 0 then hprime2 = hprime2 + 360

    if ABS(hprime1-hprime2) GT 180 then $
	meanH = (hprime1+hprime2+360)/2. $
	else meanH = (hprime1+hprime2)/2.

    T = 1.-0.17*cos((meanH-30)*!pi/180) + $
		0.24*COS(2*meanH*!pi/180) + $
		0.32*COS((3.*meanH+6)*!pi/180) - $
		0.2*COS((4.*meanH-63)*!pi/180)

    if abs(hprime2 - hprime1) LE 180 then delhp = hprime2 - hprime1 $
    else begin
	if hprime2 EQ min([hprime1,hprime2]) then $
		delhp = hprime2 - hprime1+ 360 $
	else delhp = hprime2 - hprime1 - 360
    endelse
    dellp = lab2(0) - lab1(0)
    delcp = CabPrime2 - CabPrime1
    delCapHp=2*(CabPrime2*CabPrime1)^0.5*sin((!pi/180.)*delhp/2)

    Sl = 1.+(0.015*(meanL-50.)^2.)/((20+(meanL-50.)^2.)^0.5)
    Sc = 1.+0.045*meanCprime
    Sh = 1.+0.015*meanCprime*T

    delTheta = 30.*exp(-1.*((meanH-275.)/25.)^2)
    Rc = ((meanCprime^7.)/(meanCprime^7.+25.^7.))^0.5
    Rt = -Rc*2.*SIN(2.*delTheta*!pi/180.)

    lterm = (dellp/(k[0]*Sl))^2.
	cterm = (delcp/(k[1]*Sc))^2.
	hterm = (delCapHp/(k[2]*Sh))^2.
	rotate = Rt*(delcp/(k[1]*Sc))*(delCapHp/(k[2]*Sh))
	
    de2000 = ( lterm + cterm + hterm + rotate) ^ .5

    return, de2000
end