diff --git a/R/munsell2rgb.R b/R/munsell2rgb.R
index e661b946a..e71445b7d 100644
--- a/R/munsell2rgb.R
+++ b/R/munsell2rgb.R
@@ -382,11 +382,16 @@ munsell2rgb <- function(the_hue, the_value, the_chroma, alpha = 1, maxColorValue
   # plyr::join() 2x faster than base::merge
   # data.table::merge() (with conversion to/from) 5x faster than base::merge()
   
+  
   ## TODO: maybe more efficient with keys
   # round-trip through data.table is still faster
   d <- data.table::as.data.table(d)
   munsell <- data.table::as.data.table(munsell)
   
+  
+  ## TODO: optimize by first filtering full lookup table on e.g. hue
+  
+  
   # join
   res <- merge(d, munsell, by = c('hue','value','chroma'), all.x = TRUE, sort = FALSE)
   
diff --git a/data/munsell.rda b/data/munsell.rda
index de8e7a901..7546ffdbe 100644
Binary files a/data/munsell.rda and b/data/munsell.rda differ
diff --git a/misc/utils/Munsell/interpolate-spectra.R b/misc/utils/Munsell/interpolate-spectra.R
index cb76d5daa..fad316286 100644
--- a/misc/utils/Munsell/interpolate-spectra.R
+++ b/misc/utils/Munsell/interpolate-spectra.R
@@ -1,6 +1,6 @@
 library(lattice)
 library(tactile)
-library(pbapply)
+library(purrr)
 library(reshape2)
 
 # load simplified spectra
@@ -47,7 +47,7 @@ interpolateOddChromaSpectra <- function(i) {
   # short circuit: 0 candidates for interpolation
   if(length(s.chroma) < 1)
     return(NULL)
-    
+  
   
   # setup interpolation function: natural splines
   # fit is exact at training points
@@ -89,7 +89,7 @@ interpolateOddChromaSpectra <- function(i) {
 }
 
 # do interpolation
-mm <- pblapply(m, interpolateOddChromaSpectra)
+mm <- map(m, .f = interpolateOddChromaSpectra, .progress = TRUE)
 
 # combine
 mm <- do.call('rbind', mm)
@@ -109,6 +109,16 @@ xyplot(reflectance ~ chroma | factor(wavelength), data=s,
        par.settings = tactile.theme()
 )
 
+idx <- which(m.final$hue %in% c('2.5YR') & m.final$value == 4)
+s <- m.final[idx, ]
+
+xyplot(reflectance ~ chroma | factor(wavelength), data=s, 
+       type='b', as.table=TRUE,
+       scales = list(y = list(tick.number = 10)),
+       par.settings = tactile.theme()
+)
+
+
 
 # check for reflectance <= 0
 m.final[m.final$reflectance <= 0, ]
@@ -133,7 +143,7 @@ m.final$reflectance[idx] <- min(m.final$reflectance[-idx])
 
 
 ## check: OK
-s <- subset(m.final, subset = hue == '2.5YR' & value == 4 & chroma %in% 2:4)
+s <- subset(m.final, subset = hue == '5YR' & value == 4 & chroma %in% 2:4)
 
 xyplot(reflectance ~ wavelength, data = s, 
        groups = munsell, type='b',
@@ -143,6 +153,15 @@ xyplot(reflectance ~ wavelength, data = s,
 )
 
 
+s <- subset(m.final, subset = hue == '2.5Y' & value == 4 & chroma %in% 2:4)
+
+xyplot(reflectance ~ wavelength, data = s, 
+       groups = munsell, type='b',
+       scales = list(y = list(tick.number = 10)),
+       auto.key=list(lines=TRUE, points=FALSE, cex=1, space='top', columns = 3),
+       par.settings = tactile.theme()
+)
+
 
 
 ## interpolate 2.5 value
@@ -195,7 +214,7 @@ interpolateValueSpectra <- function(i) {
 
 
 # do interpolation
-mm <- pblapply(m, interpolateValueSpectra)
+mm <- map(m, .f = interpolateValueSpectra, .progress = TRUE)
 
 # combine
 mm <- do.call('rbind', mm)
@@ -238,6 +257,11 @@ save(munsell.spectra, file = '../../../data/munsell.spectra.rda', compress = 'xz
 save(munsell.spectra.wide, file = '../../../data/munsell.spectra.wide.rda', compress = 'xz')
 
 # cleanup
-unlink(c('interpolated-Munsell-spectra-wide.rds', 'interpolated-Munsell-spectra.rds', 'simplified-Munsell-spectra.rds'))
+unlink(
+  c('interpolated-Munsell-spectra-wide.rds', 
+    'interpolated-Munsell-spectra.rds', 
+    'simplified-Munsell-spectra.rds'
+  )
+)
 
 
diff --git a/misc/utils/Munsell/local-functions.R b/misc/utils/Munsell/local-functions.R
index 1176eeb10..91c74e082 100644
--- a/misc/utils/Munsell/local-functions.R
+++ b/misc/utils/Munsell/local-functions.R
@@ -38,8 +38,47 @@ interpolateChroma <- function(m.i) {
 }
 
 
-## TODO: consider re-writing for entire range, and splinefun() based interpolation
+# 2024-09-26
+# re-write of interpolateValue() -> now safely interpolates all 0.5 values
+interpolateValue2 <- function(m.i) {
+  
+  # can only proceed with >=2 rows
+  # some combinations of hue, value, chroma have 1 row.
+  # there will be other combinations created by split() with 0 rows
+  if(nrow(m.i) < 2) {
+    return(NULL)
+  }
+  
+  # linear interpolation ~ munsell value
+  # x ~ V
+  s.1 <- splinefun(m.i$V, m.i$x)
+  # y ~ V
+  s.2 <- splinefun(m.i$V, m.i$y)
+  # Y ~ V
+  s.3 <- splinefun(m.i$V, m.i$Y)
+  
+  # all odd values
+  new.V <- seq(from = min(m.i$V) + 0.5, to = max(m.i$V) - 0.5, by = 1)
+  
+  # combine interpolated values into data.frame
+  # H, C are constant
+  m.new <- data.frame(
+    H = m.i$H[1],
+    V = new.V,
+    C = m.i$C[1],
+    p1 = s.1(new.V),
+    p2 = s.2(new.V),
+    p3 = s.3(new.V)
+  )
+  
+  names(m.new) <- c('H', 'V', 'C', 'x', 'y', 'Y')
+  
+  # only return new rows
+  return(m.new)
+}
+
 
+## NOTE: this can only interpolate between two integer values
 # 2022-03-29
 # for now only interpolating 2.5 value
 # usually interpolating xyY, 
diff --git a/misc/utils/Munsell/main.R b/misc/utils/Munsell/main.R
index 4ab53c6bb..e8b4ac12a 100644
--- a/misc/utils/Munsell/main.R
+++ b/misc/utils/Munsell/main.R
@@ -1,9 +1,11 @@
 ## Code / Data related to preparation of Munsell color interpretation in aqp
-## 2022-03-29
+## 2024-09-26
 ## D.E. Beaudette, A.G. Brown
 
 # make Munsell and related LUT
-# add neutral chips
+# + neutral chips
+# + odd chroma
+# + 0.5 value 
 # xyY [C] -> XYZ [D65] -> sRGB -> CIELAB
 source('prepare-munsell-LUT.R')
 
diff --git a/misc/utils/Munsell/munsell-LUT-2024-09-25.rds b/misc/utils/Munsell/munsell-LUT-2024-09-25.rds
new file mode 100644
index 000000000..6559ccfbc
Binary files /dev/null and b/misc/utils/Munsell/munsell-LUT-2024-09-25.rds differ
diff --git a/misc/utils/Munsell/prepare-munsell-LUT.R b/misc/utils/Munsell/prepare-munsell-LUT.R
index de7caad00..b661b635e 100644
--- a/misc/utils/Munsell/prepare-munsell-LUT.R
+++ b/misc/utils/Munsell/prepare-munsell-LUT.R
@@ -156,16 +156,17 @@ p1 <- xyplot(y ~ C | factor(V), groups = H, data = m, subset = H %in% c('2.5YR',
 p2 <- xyplot(y ~ C | factor(V), groups = H, data = m.new.chroma, subset = H %in% c('2.5YR', '2.5Y', '5G'), type = 'p', par.settings = tactile.theme(), as.table = TRUE, scales = list(alternating = 1), cex = 0.5, pch = 16)
 
 # good
-p.1 + p.2
+p1 + p2
 
 # verify odd chroma frequencies
 table(m.new.chroma$C)
 
 
+
 ## TODO:
 # * vectorize interpolateValue() or re-factor
 # * do we need multivariate interpolation?
-# 
+# * safely handle impossible interpolation (missing end points ~ some hues)
 
 .n <- length(unique(m.new.chroma$C[m.new.chroma$H %in% c('2.5Y', '2.5YR', '2.5R')]))
 .cols <- hcl.colors(n = .n, palette = 'blues3')
@@ -220,44 +221,59 @@ xyplot(
 
 
 ## interpolate all half-value chips
-# .s <- seq(from = 1.5,  to = 9.5, by = 1)
-# map(m.new.chroma, .f = interpolateValue, .progress = TRUE)
+z <- split(m.new.chroma, list(m.new.chroma$H, m.new.chroma$C))
+zz <- map(z, .f = interpolateValue2, .progress = TRUE)
 
+# remove NULLs
+# these are H/C combinations where interpolation is not possible
+idx <- which(!sapply(zz, is.null))
+zz <- zz[idx]
 
+zz <- do.call('rbind', zz)
+nrow(zz)
 
-## interpolate 2.5 values
-# only need 2 value-slices
-m.sub <- subset(m.new.chroma, V %in% c(2, 3))
-m.sub <- split(m.sub, list(m.sub$H, m.sub$C))
+# stack interpolated values
+m.new.chroma <- rbind(m.new.chroma, zz)
 
-# note: some combinations are missing values 2 AND 3
-table(sapply(m.sub, nrow))
-# 0    1    2 
-# 1102  140  718 
-
-# only process those with 2 records
-idx <- which(sapply(m.sub, nrow) == 2)
-m.sub <- m.sub[idx]
+# sort
+m.new.chroma <- m.new.chroma[order(m.new.chroma$H, m.new.chroma$V, m.new.chroma$C), ]
 
-m.2.5.values <- map(m.sub, .f = interpolateValue, .progress = TRUE)
-m.2.5.values <- do.call('rbind', m.2.5.values)
+str(m.new.chroma)
 
-# 758 rows
-nrow(m.2.5.values)
-head(m.2.5.values)
 
-## stack interpolated 2.5 values
-m.new.chroma <- rbind(m.new.chroma, m.2.5.values)
+# ## interpolate 2.5 values
+# # only need 2 value-slices
+# m.sub <- subset(m.new.chroma, V %in% c(2, 3))
+# m.sub <- split(m.sub, list(m.sub$H, m.sub$C))
+# 
+# # note: some combinations are missing values 2 AND 3
+# table(sapply(m.sub, nrow))
+# # 0    1    2 
+# # 1102  140  718 
+# 
+# # only process those with 2 records
+# idx <- which(sapply(m.sub, nrow) == 2)
+# m.sub <- m.sub[idx]
+# 
+# m.2.5.values <- map(m.sub, .f = interpolateValue, .progress = TRUE)
+# m.2.5.values <- do.call('rbind', m.2.5.values)
+# 
+# # 758 rows
+# nrow(m.2.5.values)
+# head(m.2.5.values)
+# 
+# ## stack interpolated 2.5 values
+# m.new.chroma <- rbind(m.new.chroma, m.2.5.values)
+# 
+# # sort
+# m.new.chroma <- m.new.chroma[order(m.new.chroma$H, m.new.chroma$V, m.new.chroma$C), ]
 
-# sort
-m.new.chroma <- m.new.chroma[order(m.new.chroma$H, m.new.chroma$V, m.new.chroma$C), ]
 
-str(m.new.chroma)
 
 ## graphical check
 g <- make.groups(
   m.new.chroma,
-  m.2.5.values
+  zz
 )
 
 # ok
@@ -275,6 +291,20 @@ xyplot(
   xlim = c(0, 10)
 )
 
+xyplot(
+  y ~ V | factor(C), 
+  groups = which, 
+  data = g, 
+  subset = H %in% c('2.5YR'), 
+  type = 'p', 
+  par.settings = tactile.theme(), 
+  as.table = TRUE, 
+  scales = list(alternating = 1), 
+  cex = 0.5, 
+  pch = 16, 
+  xlim = c(0, 10)
+)
+
 xyplot(
   x ~ V | factor(C), 
   groups = which, 
@@ -290,6 +320,50 @@ xyplot(
 )
 
 
+.n <- length(unique(m.new.chroma$C[m.new.chroma$H %in% c('2.5Y', '2.5YR', '2.5R')]))
+.cols <- hcl.colors(n = .n, palette = 'zissou1')
+
+xyplot(
+  x ~ V | factor(H), 
+  groups = C, 
+  data = m.new.chroma, 
+  subset = H %in% c('2.5Y', '2.5YR', '2.5R'),
+  type = 'b', 
+  par.settings = tactile.theme(
+    background = list(col = 'black'),
+    axis.text = list(col = 'white'),
+    par.xlab.text = list(col = 'white'),
+    par.ylab.text = list(col = 'white'),
+    superpose.symbol = list(col = .cols, pch = 16),
+    superpose.line = list(col = .cols, lwd = 1)
+    ), 
+  as.table = TRUE, 
+  scales = list(alternating = 1, x = list(at = seq(1, 10))),
+  panel = function(...) {
+    panel.grid(-1, -1)
+    panel.xyplot(...)
+  }
+) 
+
+xyplot(
+  x ~ V | factor(H), 
+  groups = C, 
+  data = m.new.chroma, 
+  subset = H %in% c('2.5Y', '2.5YR', '2.5R'),
+  type = 'b', 
+  par.settings = tactile.theme(
+    superpose.symbol = list(col = .cols, pch = 16),
+    superpose.line = list(col = .cols, lwd = 1)
+  ), 
+  as.table = TRUE, 
+  scales = list(alternating = 1, x = list(at = seq(1, 10))),
+  panel = function(...) {
+    panel.grid(-1, -1)
+    panel.xyplot(...)
+  }
+) 
+
+
 
 ##
 ## convert xyY [C] ---> XYZ [D65]
@@ -320,15 +394,20 @@ m.final <- data.frame(m.new.chroma, m.sRGB)
 plot_cols <- rgb(m.final$R, m.final$G, m.final$B)
 
 
-p1 <- xyplot(V ~ C | factor(H, levels=c('2.5Y', '10YR', '7.5YR', '5YR', '2.5YR', '10R')),
-             main="Common Soil Colors", 
-             data=m.final, subset=H %in% c('2.5Y', '10YR', '7.5YR', '5YR', '2.5YR', '10R') & V > 1 & V <= 8 & C <= 8, 
-             as.table=TRUE, subscripts=TRUE, xlab='Chroma', ylab='Value',
-             par.settings = tactile.theme(),
-             panel=function(x, y, subscripts, ...) 
-             {
-               panel.xyplot(x, y, pch=15, cex=2, col=plot_cols[subscripts])
-             }
+p1 <- xyplot(
+  V ~ C | factor(H, levels = c('2.5Y', '10YR', '7.5YR', '5YR', '2.5YR', '10R')),
+  main = "Common Soil Colors", 
+  data = m.final, 
+  subset = H %in% c('2.5Y', '10YR', '7.5YR', '5YR', '2.5YR', '10R') & V <= 8 & C <= 8, 
+  as.table = TRUE, 
+  subscripts = TRUE, 
+  scales = list(alternating = 1, y = list(at = 1:8)),
+  xlab = 'Chroma', 
+  ylab = 'Value',
+  par.settings = tactile.theme(),
+  panel = function(x, y, subscripts, ...) {
+    panel.xyplot(x, y, pch = 15, cex = 2, col = plot_cols[subscripts])
+  }
 )
 
 p1
@@ -374,7 +453,8 @@ n.agg.final <- n.agg.final[order(n.agg.final$V), ]
 # combine
 m.final <- rbind(m.final, n.agg.final)
 
-# 9227
+# 2022: 9,227 (2.5 value chips)
+# 2024: 15,709 (all half-value chips)
 nrow(m.final)
 
 
@@ -392,17 +472,32 @@ row.names(m.final.lab) <- NULL
 
 str(m.final.lab)
 
+
 ##
-## check
+## check differences from previous versions of the LUT
 ##
 
+## 2022:
+## the new N chips are the top differences
+## everything else has dE00 < 4
+## mostly value == 1
+
+
+## 2024:
+## dE00 > 0.4 (but all < 1.5) are 2.5 value chips
+## 
+## likely related to interpolation over full range of V vs. single, linear interpolation 2->2.5<-3
+
+
+
+
 ## make backup copy of old LUT
 # data(munsell)
-# saveRDS(munsell, file = 'munsell-LUT-2022-03-29.rds')
+# saveRDS(munsell, file = 'munsell-LUT-2024-09-25.rds')
 
-z.old <- readRDS('munsell-LUT-2022-03-29.rds')
+z.old <- readRDS('munsell-LUT-2024-09-25.rds')
 
-z <- merge(z.old, m.final.lab, by = c('hue', 'value', 'chroma'), all.x = TRUE)
+z <- merge(z.old, m.final.lab, by = c('hue', 'value', 'chroma'), all.x = TRUE, sort = FALSE)
 
 str(z)
 
@@ -431,6 +526,8 @@ for(i in 1:nrow(z)) {
   )
 }
 
+hist(d)
+
 
 # changes with dE00 > 2
 idx <- which(d > 2)
@@ -444,9 +541,19 @@ table(zz$hue)
 table(zz$value)
 table(zz$chroma)
 
-## N chips are the top differences
-## everything else is dE00 < 4
-## mostly value == 1
+
+# changes with dE00 > 0.4
+idx <- which(d > 0.4)
+zz <- z[idx, ]
+zz$dE00 <- d[idx]
+zz <- zz[order(zz$dE00, decreasing = TRUE), ]
+
+nrow(zz)
+head(zz, 20)
+
+table(zz$hue)
+table(zz$value)
+table(zz$chroma)
 
 
 
@@ -460,6 +567,12 @@ save(munsell, file = '../../../data/munsell.rda', compress = 'xz')
 
 ## install / or reload from source
 
+munsell2rgb('10YR', 9, 2, returnLAB = TRUE)
+munsell2rgb('10YR', 9.5, 2, returnLAB = TRUE)
+
+# dE00 ~ 3
+# colorContrastPlot('10YR 9/2', '10YR 9.5/2')
+
 
 munsell2rgb('10YR', 3.5, 2, returnLAB = TRUE)
 munsell2rgb('10YR', 4, 2, returnLAB = TRUE)
@@ -472,7 +585,7 @@ munsell2rgb('10YR', 5, 1, returnLAB = TRUE)
 
 
 # check neutral
-m <- sprintf('N %s/', 2:9)
+m <- sprintf('N %s/', c(2, 2.5,  3:8))
 cols <- parseMunsell(m)
 soilPalette(cols, lab = m)
 
diff --git a/misc/utils/Munsell/simplified-Munsell-spectra.rds b/misc/utils/Munsell/simplified-Munsell-spectra.rds
new file mode 100644
index 000000000..adf7030d2
Binary files /dev/null and b/misc/utils/Munsell/simplified-Munsell-spectra.rds differ