Public Function lstar(Y, Yn)
    ratio = Y / Yn
    If ratio > 0.008856 Then
        lstar = 116 * (ratio ^ (1 / 3)) - 16
    Else
        lstar = 903.3 * ratio
    End If
End Function

Public Function astar(X, Y, Xn, Yn)
    xratio = X / Xn
    yratio = Y / Yn
    third = 1 / 3
    
    If xratio > 0.008856 Then
        fx = xratio ^ third
    Else
        fx = 7.787 * xratio + (16 / 116)
    End If

    If yratio > 0.008856 Then
        fy = yratio ^ third
    Else
        fy = 7.787 * yratio + (16 / 116)
    End If
    
    astar = 500 * (fx - fy)
    
End Function

Public Function bstar(Y, Z, Yn, Zn)
    
    yratio = Y / Yn
    zratio = Z / Zn
    third = 1 / 3
    
    If zratio > 0.008856 Then
        fz = zratio ^ third
    Else
        fz = 7.787 * zratio + (16 / 116)
    End If

    If yratio > 0.008856 Then
        fy = yratio ^ third
    Else
        fy = 7.787 * yratio + (16 / 116)
    End If
    
    bstar = 200 * (fy - fz)
    
End Function

Public Function deltaE2000(l1, a1, b1, l2, a2, b2)

    ' the usual chroma calculation
    cab1 = (a1 ^ 2 + b1 ^ 2) ^ 0.5
    cab2 = (a2 ^ 2 + b2 ^ 2) ^ 0.5

    meanc = (cab1 + cab2) / 2#
    
    ' this bit of nonsense is required because VB seems to choke on the big numbers
    logmeanc = Application.Log10(meanc)
    log25 = Application.Log10(25)
    cratio = (10 ^ (logmeanc * 7 - Application.Log10((10 ^ (logmeanc * 7) + 10 ^ (log25 * 7))))) ^ 0.5
    g = 0.5 * (1 - cratio)

    ' adjust a* values
    a1 = (1 + g) * a1
    a2 = (1 + g) * a2

    cab1 = (a1 ^ 2 + b1 ^ 2#) ^ 0.5
    cab2 = (a2 ^ 2 + b2 ^ 2#) ^ 0.5

    deltac = cab1 - cab2
    
    ' calculate delta L*
    deltal = l1 - l2

    ';calculate hue angles
    h1 = Application.Degrees(Application.Atan2(a1, b1))
    h2 = Application.Degrees(Application.Atan2(a2, b2))
    If h1 < 0 Then
        h1 = h1 + 360
        End If
    If h2 < 0 Then
        h2 = h2 + 360
    End If
    
    ' calculate delta h and mean h
    deltah = h1 - h2
    meanh = h1 + h2

    If deltah > 180 Then
        deltah = deltah - 360
        meanh = meanh - 360
    End If
    
    If deltah < -180 Then
        deltah = deltah + 360
        meanh = meanh - 360
    End If

    delth = 2 * ((cab1 * cab2) ^ 0.5) * Sin(Application.Radians(deltah / 2#))

    deltal = l1 - l2
    deltaE = (deltal ^ 2 + deltac ^ 2 + delth ^ 2) ^ 0.5

    deltaha = (deltaE ^ 2 - deltal ^ 2 - deltac ^ 2)
    If deltaha < 0 Then
        deltaha = 0
    End If
    
    delth = deltaha ^ 0.5
    delth = 2 * ((cab1 * cab2) ^ 0.5) * Sin(Application.Radians(deltah / 2))
    
    meanl = (l1 + l2) / 2#
    meanc = (cab1 + cab2) / 2#
    meanh = meanh / 2#


    ' Computes the weighting functions.

    sl = 1 + ((0.015 * (meanl - 50) ^ 2) / (20 + (meanl - 50) ^ 2) ^ 0.5)
    sc = 1 + 0.045 * meanc

    ' find the hue angle goodness
    ' split these up because I don't know how to split lines in vb
    t1 = 1 - 0.17 * Cos(Application.Radians(meanh - 30))
    t2 = 0.24 * Cos(Application.Radians(2 * meanh))
    t3 = 0.32 * Cos(Application.Radians(3 * meanh + 6))
    t4 = -0.2 * Cos(Application.Radians(4 * meanh - 63))
    t = t1 + t2 + t3 + t4
    
    sh = 1 + 0.015 * meanc * t
    sr = sc * sh
    
    ' same problem as above with the big numbers
    logmeanc = Application.Log10(meanc)
    log25 = Application.Log10(25)
    rc_numerator = 10 ^ (logmeanc * 7)
    rc_denominator = rc_numerator + 10 ^ (log25 * 7)
    rc = 2 * (rc_numerator / rc_denominator) ^ 0.5

    dthet = 30 * Exp(-1 * ((meanh - 275) / 25) ^ 2)
    rt = -Sin(2 * Application.Radians(dthet)) * rc

    ' Now computes deltaV
    vl = (deltal / sl) ^ 2
    vc = (deltac / sc) ^ 2
    vh = (delth / sh) ^ 2

    deltaE2000 = (vl + vc + vh + (rt * deltac * delth / sr)) ^ 0.5

End Function