diff --git a/DESCRIPTION b/DESCRIPTION
index c867845..387d3fb 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: prospect
 Title: PROSPECT leaf radiative transfer model and inversion routines
-Version: 1.6.2
+Version: 1.7.0
 Authors@R: c(person(given = "Jean-Baptiste", 
 					    family = "Feret", 
 					    email = "jb.feret@teledetection.fr", 
diff --git a/NEWS.md b/NEWS.md
index f38b726..5e4eaa0 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+# prospect v1.7.0
+## fix
+- corrected inaccurate accounting for N when computing bNm1
+
 # prospect v1.6.2
 ## fix
 - fixed wrong merge for correction of JOSS paper
diff --git a/R/Lib_PROSPECT.R b/R/Lib_PROSPECT.R
index 2f76eec..069ffc8 100644
--- a/R/Lib_PROSPECT.R
+++ b/R/Lib_PROSPECT.R
@@ -113,7 +113,7 @@ PROSPECT <- function(SpecPROSPECT = NULL, Input_PROSPECT = NULL,
   a <- (1 + rq - tq + D) / (2 * r)
   b <- (1 - rq + tq + D) / (2 * t)
 
-  bNm1 <- b**(N - 1)
+  bNm1 <- b**(Input_PROSPECT$N - 1)
   bN2 <- bNm1**2
   a2 <- a**2
   denom <- a2 * bN2 - 1
diff --git a/R/Lib_PROSPECT_Calibration.R b/R/Lib_PROSPECT_Calibration.R
new file mode 100644
index 0000000..3f2cbeb
--- /dev/null
+++ b/R/Lib_PROSPECT_Calibration.R
@@ -0,0 +1,444 @@
+# ==============================================================================
+# prospect
+# Lib_PROSPECT_Calibration.R
+# ==============================================================================
+# PROGRAMMERS:
+# Jean-Baptiste FERET <jb.feret@teledetection.fr>
+# Florian de Boissieu <florian.deboissieu@inrae.fr>
+# Copyright 2024/02 Jean-Baptiste FERET
+# ==============================================================================
+# This Library includes functions dedicated to PROSPECT calibration
+# ==============================================================================
+
+#'
+#'
+#' @param SpecPROSPECT list. Includes spectral constants derived from
+#' SpecPROSPECT_FullRange: refractive index, specific absorption coefficients
+#' and corresponding spectral bands
+#' @param
+#' @param
+#' @param
+#' @param
+#' @param
+#' @param
+#' @param
+#' @param
+#' @param
+#' @param
+#' @param
+#' @param
+#'
+#' @return
+#' @importFrom
+#' @export
+PROSPECT_cal <- function(SpecPROSPECT = NULL, Input_PROSPECT = NULL,
+                     N = 1.5, CHL = 40.0, CAR = 8.0, ANT = 0.0, BROWN = 0.0,
+                     EWT = 0.01, LMA = NULL, PROT = 0, CBC = 0, alpha = 40.0,
+                     check = TRUE) {
+
+
+
+
+  # define PROSPECT input in a dataframe
+  Input_PROSPECT <- define_Input_PROSPECT(Input_PROSPECT, CHL, CAR,
+                                          ANT, BROWN, EWT, LMA, PROT,
+                                          CBC, N, alpha)
+  # if (check) Input_PROSPECT <- define_Input_PROSPECT(Input_PROSPECT, CHL, CAR,
+  #                                                    ANT, BROWN, EWT, LMA, PROT,
+  #                                                    CBC, N, alpha)
+  # default: simulates leaf optics using full spectral range
+  if (is.null(SpecPROSPECT)) SpecPROSPECT <- prospect::SpecPROSPECT_FullRange
+  # compute total absorption corresponding to each homogeneous layer
+  Kall <- (Input_PROSPECT$CHL * SpecPROSPECT$SAC_CHL +
+             Input_PROSPECT$CAR * SpecPROSPECT$SAC_CAR +
+             Input_PROSPECT$ANT * SpecPROSPECT$SAC_ANT +
+             Input_PROSPECT$BROWN * SpecPROSPECT$SAC_BROWN +
+             Input_PROSPECT$EWT * SpecPROSPECT$SAC_EWT +
+             Input_PROSPECT$LMA * SpecPROSPECT$SAC_LMA +
+             Input_PROSPECT$PROT * SpecPROSPECT$SAC_PROT +
+             Input_PROSPECT$CBC * SpecPROSPECT$SAC_CBC) / Input_PROSPECT$N
+
+  # Non-conservative scattering (normal case) when Kall > 0
+  j <- which(Kall <= 0)
+  t1 <- (1 - Kall) * exp(-Kall)
+  t2 <- (Kall * Kall) * expint(Kall)
+  tau <- t1 + t2
+  if (length(j)>0) tau[j] <- 1
+  # tau <- matrix(1, ncol = 1, nrow = length(t1))
+  # tau[j] <- t1[j] + t2[j]
+
+  # ***********************************************************************
+  # reflectance and transmittance of one layer
+  # ***********************************************************************
+  # Allen W.A., Gausman H.W., Richardson A.J., Thomas J.R. (1969),
+  # Interaction of isotropic ligth with a compact plant leaf, J. Opt.
+  # Soc. Am., 59(10):1376-1379.
+  # ***********************************************************************
+  # reflectivity and transmissivity at the interface
+  # ***********************************************************************
+  if (Input_PROSPECT$alpha == 40) {
+    talf <- SpecPROSPECT$calctav_40
+  } else {
+    talf <- calctav(Input_PROSPECT$alpha, SpecPROSPECT$nrefrac)
+  }
+  ralf <- 1 - talf
+  t12 <- SpecPROSPECT$calctav_90
+  r12 <- 1 - t12
+  t21 <- t12 / (SpecPROSPECT$nrefrac**2)
+  r21 <- 1 - t21
+
+  # top surface side
+  denom <- 1 - (r21 * r21 * (tau**2))
+  Ta <- (talf * tau * t21) / denom
+  Ra <- ralf + (r21 * tau * Ta)
+  # bottom surface side
+  t <- t12 * tau * t21 / denom
+  r <- r12 + (r21 * tau * t)
+
+  # ***********************************************************************
+  # reflectance and transmittance of N layers
+  # Stokes equations to compute properties of next N-1 layers (N real)
+  # Normal case
+  # ***********************************************************************
+  # Stokes G.G. (1862), On the intensity of the light reflected from
+  # or transmitted through a pile of plates, Proc. Roy. Soc. Lond.,
+  # 11:545-556.
+  # ***********************************************************************
+  D <- sqrt((1 + r + t) * (1 + r - t) * (1 - r + t) * (1 - r - t))
+  rq <- r**2
+  tq <- t**2
+  a <- (1 + rq - tq + D) / (2 * r)
+  b <- (1 - rq + tq + D) / (2 * t)
+
+  bNm1 <- b**(N - 1)
+  bN2 <- bNm1**2
+  a2 <- a**2
+  denom <- a2 * bN2 - 1
+  Rsub <- a * (bN2 - 1) / denom
+  Tsub <- bNm1 * (a2 - 1) / denom
+
+  # 	Case of zero absorption
+  j <- which(r + t >= 1)
+  Tsub[j] <- t[j] / (t[j] + (1 - t[j]) * (N - 1))
+  Rsub[j] <- 1 - Tsub[j]
+
+  # leaf reflectance and transmittance : combine top layer with next N-1 layers
+  denom <- 1 - Rsub * r
+  tran <- Ta * Tsub / denom
+  refl <- Ra + (Ta * Rsub * t) / denom
+  return(data.frame('wvl' = SpecPROSPECT$lambda,
+                    'Reflectance' = refl,
+                    'Transmittance' = tran))
+}
+
+#' computation of transmissivity of a dielectric plane surface,
+#' averaged over all directions of incidence and over all polarizations.
+#'
+#' @param alpha numeric. max incidence angle of solid angle of incident light
+#' @param nr numeric. refractive index
+#'
+#' @return numeric. Transmissivity of a dielectric plane surface
+#' @export
+calctav <- function(alpha, nr) {
+  # Stern F. (1964), Transmission of isotropic radiation across an
+  # interface between two dielectrics, Appl. Opt., 3(1):111-113.
+  # Allen W.A. (1973), Transmission of isotropic light across a
+  # dielectric surface in two and three dimensions, J. Opt. Soc. Am.,
+  # 63(6):664-666.
+  # ***********************************************************************
+
+  rd <- pi / 180
+  n2 <- nr**2
+  np <- n2 + 1
+  nm <- n2 - 1
+  a <- (nr + 1) * (nr + 1) / 2
+  k <- -(n2 - 1) * (n2 - 1) / 4
+  sa <- sin(alpha * rd)
+
+  b2 <- (sa**2) - (np / 2)
+  if (alpha == 90) {
+    b1 <- 0 * b2
+  } else {
+    b1 <- sqrt((b2**2) + k)
+  }
+  b <- b1 - b2
+  b3 <- b**3
+  a3 <- a**3
+  ts <- ((k**2) / (6 * b3) + (k / b) - b / 2) - ((k**2) / (6 * a3) + (k / a) - (a / 2))
+
+  tp1 <- -2 * n2 * (b - a) / (np**2)
+  tp2 <- -2 * n2 * np * log(b / a) / (nm**2)
+  tp3 <- n2 * ((1 / b) - (1 / a)) / 2
+  tp4 <- 16 * n2**2 * ((n2**2) + 1) * log(((2 * np * b) - (nm**2)) / ((2 * np * a) - (nm**2))) / ((np**3) * (nm**2))
+  tp5 <- 16 * (n2**3) * (1 / ((2 * np * b) - (nm**2)) - (1 / (2 * np * a - (nm**2)))) / (np**3)
+  tp <- tp1 + tp2 + tp3 + tp4 + tp5
+  tav <- (ts + tp) / (2 * (sa**2))
+  return(tav)
+}
+
+#' This function checks if the input parameters are defined as expected
+#' to run either PROSPECT-D or PROSPECT-PRO
+#' @param LMA numeric. content corresponding to LMA
+#' @param PROT numeric. content corresponding to protein content
+#' @param CBC numeric. content corresponding to carbon based constituents
+#'
+#' @return list. updated LMA, PROT and CBC
+#' @export
+
+check_version_prospect <- function(LMA, PROT, CBC){
+  # PROSPECT-D as default value
+  if (is.null(LMA) & PROT == 0 & CBC == 0) LMA <- 0.008
+  # PROSPECT-PRO if PROT or CBC are not NULL
+  if (is.null(LMA) & (PROT > 0 | CBC > 0)) LMA <- 0
+  # if calling PROSPECT-PRO (protein content or CBC defined by user)
+  # then set LMA to 0 in any case
+  if (!LMA==0 & (PROT > 0 | CBC > 0)) {
+    print_msg('version_PROSPECT')
+    LMA <- 0
+  }
+  return(list('LMA' = LMA, 'PROT' = PROT, 'CBC' = CBC))
+}
+
+
+#' This function produces a data frame from all prospect input variables if not
+#' defined already
+#'
+#'
+
+#' @param Input_PROSPECT numeric. content corresponding to LMA
+#' @param CHL numeric. Chlorophyll content (microg.cm-2)
+#' @param CAR numeric. Carotenoid content (microg.cm-2)
+#' @param ANT numeric. Anthocyanin content (microg.cm-2)
+#' @param BROWN numeric. Brown pigment content (Arbitrary units)
+#' @param EWT numeric. Equivalent Water Thickness (g.cm-2)
+#' @param LMA numeric. Leaf Mass per Area (g.cm-2)
+#' @param PROT numeric. protein content  (g.cm-2)
+#' @param CBC numeric. NonProt Carbon-based constituent content (g.cm-2)
+#' @param N numeric. Leaf structure parameter
+#' @param alpha numeric. Solid angle for incident light at surface of leaf
+#'
+#' @return list. updated LMA, PROT and CBC
+#' @export
+
+define_Input_PROSPECT <- function(Input_PROSPECT, CHL = NULL, CAR = NULL,
+                                  ANT = NULL, BROWN = NULL, EWT = NULL,
+                                  LMA = NULL, PROT = NULL, CBC = NULL,
+                                  N = NULL, alpha = NULL){
+
+  default_PROSPECT <- data.frame('CHL' = 40.0, 'CAR' = 8.0, 'ANT' = 0.0,
+                                 'BROWN' = 0.0, 'EWT' = 0.01, 'LMA' = 0.0,
+                                 'PROT'= 0.0, 'CBC' = 0.0, 'N' = 1.5,
+                                 'alpha' = 40.0)
+  if (is.null(Input_PROSPECT)){
+    dm_val <- prospect::check_version_prospect(LMA, PROT, CBC)
+    Input_PROSPECT <- data.frame('CHL' = CHL, 'CAR' = CAR, 'ANT' = ANT,
+                                 'BROWN' = BROWN, 'EWT' = EWT, 'LMA' = dm_val$LMA,
+                                 'PROT'= dm_val$PROT, 'CBC' = dm_val$CBC,
+                                 'N' = N, 'alpha' = alpha)
+  } else if (!is.null(Input_PROSPECT)){
+    missing <- which(!names(default_PROSPECT)%in%names(Input_PROSPECT))
+    if (length(missing)>0) Input_PROSPECT <- cbind(Input_PROSPECT, default_PROSPECT[missing])
+    dm_val <- prospect::check_version_prospect(Input_PROSPECT$LMA,
+                                               Input_PROSPECT$PROT,
+                                               Input_PROSPECT$CBC)
+    Input_PROSPECT$LMA <- dm_val$LMA
+    Input_PROSPECT$PROT <- dm_val$PROT
+    Input_PROSPECT$CBC <- dm_val$CBC
+  }
+  return(Input_PROSPECT)
+}
+
+
+#' This function adapts SpecPROSPECT accordingly to experimental data
+#' or to a spectral domain defined by UserDomain
+#' @param lambda numeric. Spectral bands corresponding to experimental data
+#' @param SpecPROSPECT list. Includes optical constants: refractive index,
+#' specific absorption coefficients and corresponding spectral bands
+#' @param Refl numeric. Measured reflectance data
+#' @param Tran numeric. Measured Transmittance data
+#' @param UserDomain numeric. either Lower and upper bounds for domain of
+#' interest (optional) or list of spectral bands of interest
+#' @param UL_Bounds boolean. set to TRUE if UserDomain only includes lower and
+#' upper band, set to FALSE if UserDomain is a list of spectral bands (in nm)
+#'
+#' @return list including spectral properties at the new resolution
+#' @import dplyr
+#' @export
+
+FitSpectralData <- function(lambda, SpecPROSPECT = NULL,
+                            Refl = NULL, Tran = NULL,
+                            UserDomain = NULL, UL_Bounds = FALSE) {
+  # default: simulates leaf optics using full spectral range
+  if (is.null(SpecPROSPECT)) SpecPROSPECT <- prospect::SpecPROSPECT_FullRange
+  # convert Refl and Tran into dataframe if needed
+  if (inherits(Refl, what = c('numeric', 'matrix'))) Refl <- data.frame(Refl)
+  if (inherits(Tran, what = c('numeric', 'matrix'))) Tran <- data.frame(Tran)
+  # if (class(Refl)[1]%in%c('numeric', 'matrix')) Refl <- data.frame(Refl)
+  # if (class(Tran)[1]%in%c('numeric', 'matrix')) Tran <- data.frame(Tran)
+  # Adjust LOP: check common spectral domain between PROSPECT and leaf optics
+  if (!is.null(Refl)) Refl <- Refl %>% filter(lambda%in%SpecPROSPECT$lambda)
+  if (!is.null(Tran)) Tran <- Tran %>% filter(lambda%in%SpecPROSPECT$lambda)
+  lambda <- lambda[lambda%in%SpecPROSPECT$lambda]
+  # Adjust PROSPECT
+  lb <- lambda
+  SpecPROSPECT  <- SpecPROSPECT %>% filter(SpecPROSPECT$lambda%in%lb)
+  # if UserDomain is defined
+  if (is.null(UserDomain)) UserDomain <- lambda
+  if (UL_Bounds==TRUE) UserDomain <- seq(min(UserDomain), max(UserDomain))
+  if (!is.null(Refl)) Refl <- Refl %>% filter(lambda%in%UserDomain)
+  if (!is.null(Tran)) Tran <- Tran %>% filter(lambda%in%UserDomain)
+  lambda <- lambda[lambda%in%UserDomain]
+  # Adjust PROSPECT
+  SpecPROSPECT  <- SpecPROSPECT %>% filter(SpecPROSPECT$lambda%in%UserDomain)
+  if (any(!UserDomain%in%lambda)){
+    message('leaf optics out of range defined by UserDomain')
+  }
+  RT <- reshape_lop4inversion(Refl = Refl,
+                              Tran = Tran,
+                              SpecPROSPECT = SpecPROSPECT)
+  return(list("SpecPROSPECT" = SpecPROSPECT, "lambda" = lambda,
+              "Refl" = RT$Refl, "Tran" = RT$Tran, "nbSamples" = RT$nbSamples))
+}
+
+#' computation of a LUT of leaf optical properties using a set of
+#' leaf chemical & structural parameters
+#'
+#' @param Input_PROSPECT dataframe. list of PROSPECT input parameters.
+#' @param SpecPROSPECT list. spectral constants
+#' refractive index, specific absorption coefficients & spectral bands
+#'
+#' @return list. LUT including leaf reflectance and transmittance
+#' @export
+PROSPECT_LUT <- function(Input_PROSPECT, SpecPROSPECT = NULL) {
+
+  if (is.null(SpecPROSPECT)) SpecPROSPECT <- prospect::SpecPROSPECT_FullRange
+  # expected PROSPECT input parameters
+  ExpectedParms <- data.frame('CHL' = 0, 'CAR' = 0, 'ANT' = 0, 'BROWN' = 0,
+                              'EWT' = 0, 'LMA' = 0, 'PROT' = 0, 'CBC' = 0,
+                              'N' = 1.5, 'alpha' = 40)
+  # parameters provided
+  inOK <- names(Input_PROSPECT)
+  # identify missing elements
+  Parm2Add <- which(!names(ExpectedParms) %in% inOK)
+  # check if all parameters are included.
+  if (length(Parm2Add) > 0) print_msg(cause = 'Missing_Input',
+                                      args = list('Input' = inOK,
+                                                  'Expected' = ExpectedParms))
+
+  # re-order missing elements the end of the list using default value
+  Input_PROSPECT <- Complete_Input_PROSPECT(Input_PROSPECT = Input_PROSPECT,
+                                            Parm2Add = Parm2Add,
+                                            ExpectedParms = ExpectedParms)
+  # print number of samples to be simulated
+  nbSamples <- nrow(Input_PROSPECT)
+  messageLUT <- paste('A LUT with', nbSamples, 'samples will be produced')
+  cat(colour_to_ansi('green'), messageLUT, "\033[0m\n")
+
+  # run PROSPECT for nbSamples
+  run_list_PROSPECT <- function(Input_PROSPECT, SpecPROSPECT){
+    LUT_tmp <- do.call(PROSPECT,
+                       c(list(SpecPROSPECT = SpecPROSPECT),
+                         Input_PROSPECT))
+    return(LUT_tmp)
+  }
+  indiv_leaves <- split(Input_PROSPECT,
+                        factor(seq(1:nbSamples)))
+  LUT_tmp <- lapply(X = indiv_leaves,
+                    FUN = run_list_PROSPECT,
+                    SpecPROSPECT = SpecPROSPECT)
+  LUT_Refl <- as.data.frame(lapply(LUT_tmp,'[[', 'Reflectance'))
+  LUT_Tran <- as.data.frame(lapply(LUT_tmp,'[[', 'Transmittance'))
+  names(LUT_Refl) <- names(LUT_Tran) <- paste0('sample_',seq(1,nbSamples))
+  return(list('Reflectance' = LUT_Refl,
+              'Transmittance' = LUT_Tran,
+              'Input_PROSPECT' = Input_PROSPECT))
+}
+
+#' Complete the list of PROSPECT parameters with default values
+#'
+#' @param urldb character. URL for online repository where to download data
+#' @param dbName character. name of the database available online
+#'
+#' @return list. Includes leaf chemistry, refl, tran & number of samples
+#' @export
+
+download_LeafDB <- function(urldb = NULL,
+                            dbName = 'ANGERS'){
+  # repository where data are stored
+  if (is.null(urldb)) urldb <- 'https://gitlab.com/jbferet/myshareddata/raw/master/LOP/'
+  # download leaf chemistry and optical properties
+  DataBioch <- data.table::fread(file.path(urldb,dbName,'DataBioch.txt'))
+  Refl <- data.table::fread(file.path(urldb,dbName,'ReflectanceData.txt'))
+  Tran <- data.table::fread(file.path(urldb,dbName,'TransmittanceData.txt'))
+  # Get wavelengths corresponding to the reflectance & transmittance measurements
+  lambda <- Refl$wavelength
+  Refl$wavelength <- Tran$wavelength <- NULL
+  # Get the number of samples
+  nbSamples <- ncol(Refl)
+  return(list('DataBioch' = DataBioch,
+              'lambda' = lambda,
+              'Refl' = Refl,
+              'Tran' = Tran,
+              'nbSamples' = nbSamples))
+}
+
+
+#' Complete the list of PROSPECT parameters with default values
+#'
+#' @param Input_PROSPECT input parameters sent to PROSPECT by user
+#' @param Parm2Add Parameters to be added to input parameters
+#' @param ExpectedParms full set of parameters expected to run PROSPECT
+#'
+#' @return Input_PROSPECT
+#' @export
+
+Complete_Input_PROSPECT <- function(Input_PROSPECT, Parm2Add, ExpectedParms) {
+  ii <- 0
+  nbSamples <- length(Input_PROSPECT[[1]])
+  nbInputs <- length(Input_PROSPECT)
+  for (i in Parm2Add) {
+    ii <- ii + 1
+    nbInputs <- nbInputs + 1
+    Input_PROSPECT[[nbInputs]] <- matrix(ExpectedParms[[i]],
+                                         ncol = 1, nrow = nbSamples)
+    names(Input_PROSPECT)[[nbInputs]] <- names(ExpectedParms)[[i]]
+  }
+  return(data.frame('CHL' = matrix(Input_PROSPECT$CHL, ncol = 1),
+                    'CAR' = matrix(Input_PROSPECT$CAR, ncol = 1),
+                    'ANT' = matrix(Input_PROSPECT$ANT, ncol = 1),
+                    'BROWN' = matrix(Input_PROSPECT$BROWN, ncol = 1),
+                    'EWT' = matrix(Input_PROSPECT$EWT, ncol = 1),
+                    'LMA' = matrix(Input_PROSPECT$LMA, ncol = 1),
+                    'PROT' = matrix(Input_PROSPECT$PROT, ncol = 1),
+                    'CBC' = matrix(Input_PROSPECT$CBC, ncol = 1),
+                    'N' = matrix(Input_PROSPECT$N, ncol = 1),
+                    'alpha' = matrix(Input_PROSPECT$alpha, ncol = 1)))
+}
+
+
+#' Convert plain text colour to ANSI code
+#'
+#' @param colour colour in plain text ("red", "green", etc.) to convert to ANSI
+#'
+#' @return string representing provided colour as ANSI encoding
+#'
+#' @examples
+#' colour_to_ansi("red") # gives: "\033[31m"
+#' @export
+
+colour_to_ansi <- function(colour) {
+  # Note ANSI colour codes
+  colour_codes <- list("black" = 30,
+                       "red" = 31,
+                       "green" = 32,
+                       "yellow" = 33,
+                       "blue" = 34,
+                       "magenta" = 35,
+                       "cyan" = 36,
+                       "white" = 37)
+
+  # Create ANSI version of colour
+  ansi_colour <- paste0("\033[", colour_codes[[colour]], "m")
+  return(ansi_colour)
+}
+