diff --git a/Tetmer/DESCRIPTION b/Tetmer/DESCRIPTION
index bbda728..29785a6 100644
--- a/Tetmer/DESCRIPTION
+++ b/Tetmer/DESCRIPTION
@@ -1,11 +1,10 @@
 Package: Tetmer
 Type: Package
 Title: Fitting Population Parameters to k-mer Spectra
-Version: 2.0.0
+Version: 2.1.0
 Author: Hannes Becher
 Maintainer: The package maintainer <hannesbecher@gmx.de>
-Description: More about what it does (maybe more than one line)
-    Use four spaces when indenting paragraphs within the Description.
+Description: This package contains a shiny app for fitting population parameters to k-mer spectra.
 License: GPL (>= 3)
 Depends: 
     R (>= 2.10),
diff --git a/Tetmer/R/data.R b/Tetmer/R/data.R
index 0506342..132e9c1 100644
--- a/Tetmer/R/data.R
+++ b/Tetmer/R/data.R
@@ -6,10 +6,11 @@
 #'
 #' @format A \code{spectrum} object containing a name and data frame:
 #' \describe{
-#'   \item{name}{The name E028}
-#'   \item{data}{A dataframe with coulmns \describe{
-#'   \item{mult}{K-mer multiplicity}
-#'   \item{count}{The number of different k-mer at a given multiplicity}
+#'   \item{`name`}{The name 'E. arctica, E028'}
+#'   \item{`k`}{The k-mer length of 27}
+#'   \item{`data`}{A dataframe with coulmns \describe{
+#'   \item{`mult`}{K-mer multiplicity}
+#'   \item{`count`}{The number of different k-mer at a given multiplicity}
 #'   }
 #'   }
 #' }
@@ -17,6 +18,7 @@
 #'
 #' @examples plot(E028, log="xy")
 #' plot(E028, xlim=c(0, 200), ylim=c(0, 10000000))
+#' tetmer(E028)
 "E028"
 
 
@@ -27,16 +29,18 @@
 #'
 #' @format A \code{spectrum} object containing a name and data frame:
 #' \describe{
-#'   \item{name}{The name E028}
-#'   \item{data}{A dataframe with coulmns \describe{
-#'   \item{mult}{K-mer multiplicity}
-#'   \item{count}{The number of different k-mer at a given multiplicity}
+#'   \item{`name`}{The name 'E. anglica, E030'}
+#'   \item{`k`}{The k-mer length of 21}
+#'   \item{`data`}{A dataframe with coulmns \describe{
+#'   \item{`mult`}{K-mer multiplicity}
+#'   \item{`count`}{The number of different k-mer at a given multiplicity}
 #'   }
 #'   }
 #' }
-#' @source Becher et al. (2020) \url{https://doi.org/10.1016/j.xplc.2020.100105}
+#' @source Source!
 #'
 #' @examples plot(E030, log="xy")
 #' plot(E030, xlim=c(0, 200), ylim=c(0, 10000000))
+#' tetmer(E030)
 "E030"
 
diff --git a/Tetmer/R/formulae.R b/Tetmer/R/formulae.R
index 73fb785..4d3d21b 100644
--- a/Tetmer/R/formulae.R
+++ b/Tetmer/R/formulae.R
@@ -64,8 +64,10 @@ factorTraaa <- expression(
 
 factorTraab <- expression(
   c(
-    tth^(-tdiverg*tth)*(-2-4*tth+exp(tdiverg*tth)*(2+7*tth+3*tth^2)) / (2+3*tth+tth^2),
+    exp(-tdiverg*tth)*(-2-4*tth+exp(tdiverg*tth)*(2+7*tth+3*tth^2)) / (2+3*tth+tth^2),
+
     (1 + ((2*exp(-tdiverg*tth) * (tth-1))/(2+tth)))/ (1 + tth),
+
     2*exp(-tdiverg*tth) / (2 + 3*tth + tth^2)
   )
 )
diff --git a/Tetmer/R/tetmer.R b/Tetmer/R/tetmer.R
index 68a3c12..0df6c01 100644
--- a/Tetmer/R/tetmer.R
+++ b/Tetmer/R/tetmer.R
@@ -43,7 +43,7 @@ tet.server <- function(input, output) {
     factors <<- getFactors(input)
 
     # read spectrum
-    spec <- prepare.spectrum(sp)
+    spec <- prepare.spectrum(.sp)
 
     if(input$fitmod=="man"){
 
@@ -51,7 +51,8 @@ tet.server <- function(input, output) {
       plotSpecApp(input, spec)
       addvertlines(input) # vertical lines, multiples of n
       pointsFit(input)
-      output$outText <- textOut(input)
+      output$outText <- textOut(input, 0, spec)
+      #gsFit(input, 0) # debug
       return()
 
     }
@@ -61,20 +62,24 @@ tet.server <- function(input, output) {
       plotSpecApp(input, spec)
       minFun <<- makeMinFun(input)
       startingVals <<- getStartingVals(input)
-      optimised <- doOptimisation(input)
+      optimised <- doOptimisation(input, spec)
       addvertlines(input, optimised)
       pointsFit(input, optimised)
-      output$outText <- textOut(input, optimised)
+      pointsExtrap(input, optimised)
+      pointsContam(input, optimised, spec)
+      output$outText <- textOut(input, optimised, spec)
+      #gsFit(input, optimised) # debug
 
     }# if auto
 
+
   })
 }
 
 
 # The user interface
-tet.ui <- fluidPage(titlePanel("Tetmer"),
-                    "Fitting paramters to k-mer spectra (by Hannes Becher)",
+tet.ui <- fluidPage(titlePanel("Tetmer v2.1.0"),
+                    "Fitting population paramters to k-mer spectra (by Hannes Becher)",
                     fluidRow(
                       column(8, plotOutput('plot')),
                       column(4,
@@ -109,11 +114,11 @@ tet.ui <- fluidPage(titlePanel("Tetmer"),
                       column(3,
                              conditionalPanel(condition = "input.fitmod == 'man'",
                                               wellPanel(
-                                                h4("3rd: Adjust parameters, make the fit match the data"),
-                                                numericInput('tkcov', 'k-mer  multiplicity', tkcov),
+                                                h4("3rd: Param ranges"),
+                                                numericInput('tkcov', 'Monoploid k-mer multiplicity', tkcov),
                                                 numericInput('tbias', 'Peak width', tbias),
-                                                numericInput('tth', 'theta', tth),
-                                                numericInput('tyadj', 'Haploid non-rep GS (Mbp)', tyadj)
+                                                numericInput('tth', 'θ', tth),
+                                                numericInput('tyadj', 'Monoploid non-rep GS (Mbp)', tyadj)
                                               ))),
                       column(3,
                              conditionalPanel(condition = "(input.fitmod == 'man') && (['tal', 'traab'].includes(input.mod))",
@@ -133,7 +138,7 @@ tet.ui <- fluidPage(titlePanel("Tetmer"),
                                               ))),
                       column(3,
                              conditionalPanel(condition = "input.fitmod == 'auto'",
-                                              wellPanel(h4("3rd: adjust parameter ranges for fitting"),
+                                              wellPanel(h4("3rd: Param ranges"),
 
                                                         sliderInput('akcov', 'k-mer  multiplicity',
                                                                     min=5, max = 200,
@@ -141,10 +146,10 @@ tet.ui <- fluidPage(titlePanel("Tetmer"),
                                                         sliderInput('abias', 'Peak width',
                                                                     min=0.01, max = 4,
                                                                     value=c(abiasl, abiash)),
-                                                        sliderInput('ath', "log10 of theta",
+                                                        sliderInput('ath', "log10 of θ",
                                                                     min=-4, max = 1, step = 0.05,
                                                                     value=c(athl, athh)),
-                                                        sliderInput('ayadj', 'Haploid non-rep GS (Mbp)',
+                                                        sliderInput('ayadj', 'Monoploid non-rep GS (Mbp)',
                                                                     min=1, max =2000,
                                                                     value=c(agsl, agsh))
                                               ))),
@@ -171,7 +176,7 @@ tet.ui <- fluidPage(titlePanel("Tetmer"),
 #' @examples \dontrun{tetmer(E028)}
 #' \dontrun{tetmer(E030)}
 tetmer <- function(sp){
-  sp <<- sp
+  .sp <<- sp
   shinyApp(ui = tet.ui, server = tet.server)
 }
 
@@ -180,9 +185,10 @@ tetmer <- function(sp){
 #' @slot name A string indicating the name of the spectrum (or sample)
 #'
 #' @slot data A `data.frame` with numeric cloumns mult and count.
+#' @slot k A `numeric` indiciating the k-mer length.
 #'
 #' @export
-setClass("spectrum", slots=list(name="character", data="data.frame"))
+setClass("spectrum", slots=list(name="character", data="data.frame", k="numeric"))
 
 
 
@@ -192,7 +198,8 @@ setClass("spectrum", slots=list(name="character", data="data.frame"))
 #'
 #'
 #' @param f A string indicating the path to th spectrum file
-#' @param nam (otional) A string indicaing the name of the spectrum
+#' @param nam (otional) A string indicating the name of the spectrum
+#' @param k (otional) A numeric indicating the k-mer length
 #' @param ... keyword arguments to be passed to `read.table()`
 #' @return A two-element list comprising a names string and a two-column data frame of the
 #' k-mer spectrum
@@ -200,20 +207,22 @@ setClass("spectrum", slots=list(name="character", data="data.frame"))
 #' @importFrom methods new
 #' @importFrom utils read.table
 #' @examples
-#' testdf <- data.frame(mult=1:100, count = 100:1)
-#' tf <- tempfile()
-#' write.table(testdf, tf)
-#' sp <- read.spectrum(tf, nam = "Test spectrum")
-#' unlink(tf)
-#' plot(sp)
+#' \dontrun{testdf <- data.frame(mult=1:100, count = 100:1)}
+#' \dontrun{tf <- tempfile()}
+#' \dontrun{write.table(testdf, tf)}
+#' \dontrun{sp <- read.spectrum(tf, nam = "Test spectrum")}
+#' \dontrun{unlink(tf)}
+#' \dontrun{plot(sp)}
 read.spectrum <- function(f,
-                          nam="spectrum",
+                          nam="MySpectrum",
+                          k=0,
                           ...){
   sp = read.table(f, ...)
   names(sp) = c("mult", "count")
   return(new("spectrum",
              name=nam,
-             data=sp
+             data=sp,
+             k=k
   )
   )
 }
@@ -321,6 +330,13 @@ plotSpecApp <- function(input, spec){
       plot(spec, xlim=c(0,input$txmax), ylim=c(0, input$tymax*1000),
            xlab="K-mer multiplicity (coverage)", ylab="K-mer count", type = 'n')
     }
+    legend("topright",
+           col=c(1, 2),
+           lwd=c(1, 2),
+           lty=c(0, 1),
+           pch=c(1, NA),
+           legend=c("Data", "Fit")
+           )
 
   }
   if(input$fitmod=="auto"){
@@ -331,9 +347,15 @@ plotSpecApp <- function(input, spec){
            type = 'n')
     }
     abline(v=input$axrange)
+    legend("topright",
+           col=c(1, 2, 2, 4),
+           lwd=c(1, 2, 2, 2),
+           lty=c(0, 1, 2, 1),
+           pch=c(1, NA, NA, NA),
+           legend=c("Data", "Fit", "Extrapolation", "Contaminations"))
 
   }
-  legend("topleft", col=c(1,2), lwd=2, lty=c(0,1), pch=c(1,NA), legend=c("Data","Fit"))
+
 
 
 }
@@ -369,11 +391,13 @@ pointsFit <- function(input, optimised=0){
   }
   if(input$fitmod=="auto"){
     if(input$mod %in% c("d", "tau", "traaa")){
+      # fit within range
     points(input$axrange[1]:input$axrange[2],
            colSums(eval(probs, envir = list(txmin=input$axrange[1], txmax=input$axrange[2], tkcov=optimised$par[1], tbias=optimised$par[2])) *
                      eval(factors, envir=list(tth=optimised$par[3])))*optimised$par[4]*1000000,
            col="red", type = 'l', lty=1, lwd=2
     )
+
     }
     if(input$mod %in% c("tal", "traab")){
       points(input$axrange[1]:input$axrange[2],
@@ -388,7 +412,62 @@ pointsFit <- function(input, optimised=0){
                                        tdiverg=optimised$par[5]))
                      )*optimised$par[4]*1000000,
              col="red", type = 'l', lty=1, lwd=2
+
+      )
+
+
+
+    }
+  }
+}
+
+# for debug only
+gsFit <- function(input, optimised=0){
+
+  if(input$fitmod=="man"){
+    if(input$mod %in% c("d", "tau", "traaa")){
+      print(sum(1:1000 *
+        colSums(eval(probs, envir = list(txmin= 1, txmax=1000, tkcov=input$tkcov, tbias=input$tbias)) *
+                  eval(factors, envir=list(tth=input$tth)))*input$tyadj*1000000
+      ) / input$tkcov /1000000
       )
+    }
+    if(input$mod %in% c("tal", "traab")){
+      print(sum(1:1000 *
+        colSums(eval(probs, envir = list(txmin= 1, txmax=1000, tkcov=input$tkcov, tbias=input$tbias)) *
+                  eval(factors, envir=list(tth=input$tth, tdiverg=input$tdiverg)))*input$tyadj*1000000
+      ) / input$tkcov /1000000
+      )
+    }
+  }
+  if(input$fitmod=="auto"){
+    if(input$mod %in% c("d", "tau", "traaa")){
+      # fit within range
+      print(sum(1:1000 *
+             colSums(eval(probs, envir = list(txmin=1,
+                                              txmax=1000,
+                                              tkcov=optimised$par[1],
+                                              tbias=optimised$par[2])) *
+                       eval(factors, envir=list(tth=optimised$par[3]))
+                     )*optimised$par[4]*1000000
+      )  / optimised$par[1] / 1000000
+      )
+
+    }
+    if(input$mod %in% c("tal", "traab")){
+      print(sum(1:1000 *
+             colSums(eval(probs,
+                          envir = list(txmin=1,
+                                       txmax=1000,
+                                       tkcov=optimised$par[1],
+                                       tbias=optimised$par[2]
+                          )) *
+                       eval(factors,
+                            envir=list(tth=optimised$par[3],
+                                       tdiverg=optimised$par[5]))
+             )*optimised$par[4]*1000000
+      )   / optimised$par[1] / 1000000)
+
 
 
     }
@@ -396,6 +475,70 @@ pointsFit <- function(input, optimised=0){
 }
 
 
+
+# only run in auto mode
+pointsExtrap <- function(input, optimised){
+  if(input$mod %in% c("d", "tau", "traaa")){
+    points(1:input$axrange[1],
+           colSums(eval(probs, envir = list(txmin=1, txmax=input$axrange[1], tkcov=optimised$par[1], tbias=optimised$par[2])) *
+                     eval(factors, envir=list(tth=optimised$par[3])))*optimised$par[4]*1000000,
+           col="red", type = 'l', lty=2, lwd=2
+    )
+  }
+  if(input$mod %in% c("tal", "traab")){
+    points(1:input$axrange[1],
+           colSums(eval(probs,
+                        envir = list(txmin=1,
+                                     txmax=input$axrange[1],
+                                     tkcov=optimised$par[1],
+                                     tbias=optimised$par[2]
+                        )) *
+                     eval(factors,
+                          envir=list(tth=optimised$par[3],
+                                     tdiverg=optimised$par[5]))
+           )*optimised$par[4]*1000000,
+           col="red", type = 'l', lty=2, lwd=2
+    )
+  }
+}
+
+pointsContam <- function(input, optimised, spect){
+  if(input$mod %in% c("d", "tau", "traaa")){
+    # contaminations (i.e., spectrum - fit)
+    points(1:input$axrange[1],
+
+           spect@data[1:input$axrange[1], 2] -
+           colSums(eval(probs,
+                        envir = list(txmin=1,
+                                     txmax=input$axrange[1],
+                                     tkcov=optimised$par[1],
+                                     tbias=optimised$par[2])
+                        ) *
+                     eval(factors, envir=list(tth=optimised$par[3]))
+                   ) * optimised$par[4]*1000000,
+           type='l', col=4, lwd=2
+    )
+
+  }
+  if(input$mod %in% c("tal", "traab")){
+
+    points(1:input$axrange[1],
+           spect@data[1:input$axrange[1], 2] -
+           colSums(eval(probs,
+                        envir = list(txmin=1,
+                                     txmax=input$axrange[1],
+                                     tkcov=optimised$par[1],
+                                     tbias=optimised$par[2]
+                        )) *
+                     eval(factors,
+                          envir=list(tth=optimised$par[3],
+                                     tdiverg=optimised$par[5]))
+                   ) * optimised$par[4] * 1000000,
+           type = 'l', col = 4, lwd=2
+    )
+}
+}
+
 #' Generate text for Tetmer window
 #'
 #' @param input Input from the GUI (contains plotting range, model, etc.)
@@ -404,14 +547,193 @@ pointsFit <- function(input, optimised=0){
 #' @return A string of the fitted parameters produced by \code{renderText()}
 #' to be displayed in the Tetmer window.
 #' @keywords internal
-textOut <- function(input, optimised=0){
+textOut <- function(input, optimised, spec){
+  if(spec@k > 0){
+    if(input$fitmod == "man"){
+      if(input$mod=="d"){
+        return(renderText(paste(
+          "DIPLOID MODEL, MANUAL FIT",
+          "\n         k-mer length:", spec@k,
+          "\n  monoploid k-mer cov:", input$tkcov,
+          "\n          θ per k-mer:", input$tth,
+          "\n     θ per nucleotide:", round(input$tth/spec@k, 4),
+          "\n     non-rep GS (Mbp):", input$tyadj,
+          "\n    bias (peak width):", input$tbias
+        ))
+        )
+      }
+      if(input$mod=="tau"){
+        return(
+          renderText(paste(
+            "AUTOTETRAPLOID MODEL, MANUAL FIT",
+            "\n         k-mer length:", spec@k,
+            "\n  monoploid k-mer cov:", input$tkcov,
+            "\n          θ per k-mer:", input$tth,
+            "\n     θ per nucleotide:", round(input$tth/spec@k, 4),
+            "\n     non-rep GS (Mbp):", input$tyadj,
+            "\n    bias (peak width):", input$tbias
+          ))
+        )
+      }
+      if(input$mod=="traaa"){
+        return(
+          renderText(paste(
+            "AUTOTRIPLOID MODEL, MANUAL FIT",
+            "\n         k-mer length:", spec@k,
+            "\n  monoploid k-mer cov:", input$tkcov,
+            "\n          θ per k-mer:", input$tth,
+            "\n     θ per nucleotide:", round(input$tth/spec@k, 4),
+            "\n     non-rep GS (Mbp):", input$tyadj,
+            "\n    bias (peak width):", input$tbias
+          ))
+        )
+      }
+      if(input$mod=="traab"){
+        return(
+          renderText(paste(
+            "ALLOTRIPLOID MODEL, MANUAL FIT",
+            "\n         k-mer length:", spec@k,
+            "\n  monoploid k-mer cov:", input$tkcov,
+            "\n          θ per k-mer:", input$tth,
+            "\n     θ per nucleotide:", round(input$tth/spec@k, 4),
+            "\n                    T:", input$tdiverg,
+            "\n     non-rep GS (Mbp):", input$tyadj,
+            "\n    bias (peak width):", input$tbias,
+            "\n     diverg per k-mer:", round(input$tth*input$tdiverg, 4),
+            "\ndiverg per nucleotide:", round(input$tth*input$tdiverg/spec@k, 4)
+          ))
+        )
+      }
+      if(input$mod=="tal"){
+        return(
+          renderText(paste(
+            "ALLOTETRAPLOID MODEL, MANUAL FIT",
+            "\n         k-mer length:", spec@k,
+            "\n  monoploid k-mer cov:", input$tkcov,
+            "\n          θ per k-mer:", input$tth,
+            "\n     θ per nucleotide:", round(input$tth/spec@k, 4),
+            "\n                    T:", input$tdiverg,
+            "\n     non-rep GS (Mbp):", input$tyadj,
+            "\n    bias (peak width):", input$tbias,
+            "\n     diverg per k-mer:", round(input$tth*input$tdiverg, 4),
+            "\ndiverg per nucleotide:", round(input$tth*input$tdiverg/spec@k, 4)
+          ))
+        )
+      }
+    } # if man
+    if(input$fitmod == "auto"){
+      if(input$mod=="tal"){
+        return(
+          renderText(paste(
+            "ALLOTETRAPLOID MODEL, AUTO FITTED",
+            "\n         k-mer length:", spec@k,
+            "\n  monoploid k-mer cov:", round(optimised$par[1],1),
+            "\n          θ per k-mer:", round(optimised$par[3],4),
+            "\n     θ per nucleotide:", round(optimised$par[3] / spec@k,5),
+            "\n                    T:", round(optimised$par[5],2),
+            "\n     non-rep GS (Mbp):", round(optimised$par[4],1),
+            "\n    bias (peak width):", round(optimised$par[2],1),
+            "\n diverg per k-mer:", round(optimised$par[3]*optimised$par[5], 4),
+            "\ndiverg per nucleotide:", round(optimised$par[3]*optimised$par[5]/spec@k, 4),
+            "\n\nSTARTING RANGES (MIN MAX)",
+            "\n  monoploid k-mer cov:", input$akcov[1], input$akcov[2],
+            "\n    log10 θ per k-mer:", input$ath[1], input$ath[2],
+            "\n                    T:", input$adiv[1], input$adiv[2],
+            "\n     non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
+            "\n    bias (peak width):", input$abias[1], input$abias[2],
+            "\n              x range:", input$axrange[1],input$axrange[2]
+          ))
+        )
+      } # if allotet
+      if(input$mod=="traab"){
+        return(
+          renderText(paste(
+            "ALLOTRIPLOID MODEL, AUTO FITTED",
+            "\n         k-mer length:", spec@k,
+            "\n  monoploid k-mer cov:", round(optimised$par[1],1),
+            "\n          θ per k-mer:", round(optimised$par[3],4),
+            "\n     θ per nucleotide:", round(optimised$par[3] / spec@k,5),
+            "\n                    T:", round(optimised$par[5],2),
+            "\n     non-rep GS (Mbp):", round(optimised$par[4],1),
+            "\n    bias (peak width):", round(optimised$par[2],1),
+            "\n diverg per k-mer:", round(optimised$par[3]*optimised$par[5], 4),
+            "\ndiverg per nucleotide:", round(optimised$par[3]*optimised$par[5]/spec@k, 4),
+            "\n\nSTARTING RANGES (MIN MAX)",
+            "\n  monoploid k-mer cov:", input$akcov[1], input$akcov[2],
+            "\n    log10 θ per k-mer:", input$ath[1], input$ath[2],
+            "\n                    T:", input$adiv[1], input$adiv[2],
+            "\n     non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
+            "\n    bias (peak width):", input$abias[1], input$abias[2],
+            "\n              x range:", input$axrange[1],input$axrange[2]
+          ))
+        )
+      } # if allotrip
+      if(input$mod=="tau"){
+        return(
+          renderText(paste(
+            "AUTOTETRAPLOID MODEL, AUTO FITTED",
+            "\n         k-mer length:", spec@k,
+            "\n  monoploid k-mer cov:", round(optimised$par[1],1),
+            "\n          θ per k-mer:", round(optimised$par[3],4),
+            "\n     θ per nucleotide:", round(optimised$par[3] / spec@k,5),
+            "\n     non-rep GS (Mbp):", round(optimised$par[4],1),
+            "\n    bias (peak width):", round(optimised$par[2],1),
+            "\n\nSTARTING RANGES (MIN MAX)",
+            "\n  monoploid k-mer cov:", input$akcov[1], input$akcov[2],
+            "\n    log10 θ per k-mer:", input$ath[1], input$ath[2],
+            "\n     non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
+            "\n    bias (peak width):", input$abias[1], input$abias[2],
+            "\n              x range:", input$axrange[1],input$axrange[2]
+          ))
+        )
+      } # if tau
+      if(input$mod=="traaa"){
+        return(
+          renderText(paste(
+            "AUTOTRIPLOID MODEL, AUTO FITTED",
+            "\n         k-mer length:", spec@k,
+            "\n monoploid k-mer cov:", round(optimised$par[1],1),
+            "\n          θ per k-mer:", round(optimised$par[3],4),
+            "\n     θ per nucleotide:", round(optimised$par[3] / spec@k,5),
+            "\n     non-rep GS (Mbp):", round(optimised$par[4],1),
+            "\n    bias (peak width):", round(optimised$par[2],1),
+            "\n\nSTARTING RANGES (MIN MAX)",
+            "\n  monoploid k-mer cov:", input$akcov[1], input$akcov[2],
+            "\n    log10 θ per k-mer:", input$ath[1], input$ath[2],
+            "\n     non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
+            "\n    bias (peak width):", input$abias[1], input$abias[2],
+            "\n              x range:", input$axrange[1],input$axrange[2]
+          ))
+        )
+      } # if traab
+      if(input$mod=="d"){
+        return(
+          renderText(paste(
+            "DIPLOID MODEL, AUTO FITTED",
+            "\n         k-mer length:", spec@k,
+            "\n  monoploid k-mer cov:", round(optimised$par[1],1),
+            "\n          θ per k-mer:", round(optimised$par[3],4),
+            "\n     θ per nucleotide:", round(optimised$par[3] / spec@k,5),
+            "\n     non-rep GS (Mbp):", round(optimised$par[4],1),
+            "\n    bias (peak width):", round(optimised$par[2],1),
+            "\n\nSTARTING RANGES (MIN MAX)",
+            "\n  monoploid k-mer cov:", input$akcov[1], input$akcov[2],
+            "\n    log10 θ per k-mer:", input$ath[1], input$ath[2],
+            "\n     non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
+            "\n    bias (peak width):", input$abias[1], input$abias[2],
+            "\n              x range:", input$axrange[1],input$axrange[2]
+          ))
+        )
+      } #if d
+    } # if auto fit
+  } # if k > 0
   if(input$fitmod == "man"){
     if(input$mod=="d"){
       return(renderText(paste(
         "DIPLOID MODEL, MANUAL FIT",
-        "\n    haploid k-mer cov:", input$tkcov,
-        "\n      per k-mer theta:", input$tth,
-        "\nhapl non-rep GS (Mbp):", input$tyadj,
+        "\n  monoploid k-mer cov:", input$tkcov,
+        "\n          θ per k-mer:", input$tth,
+        "\n     non-rep GS (Mbp):", input$tyadj,
         "\n    bias (peak width):", input$tbias
       ))
       )
@@ -420,9 +742,9 @@ textOut <- function(input, optimised=0){
       return(
         renderText(paste(
           "AUTOTETRAPLOID MODEL, MANUAL FIT",
-          "\n    haploid k-mer cov:", input$tkcov,
-          "\n      per k-mer theta:", input$tth,
-          "\nhapl non-rep GS (Mbp):", input$tyadj,
+          "\n  monoploid k-mer cov:", input$tkcov,
+          "\n          θ per k-mer:", input$tth,
+          "\n     non-rep GS (Mbp):", input$tyadj,
           "\n    bias (peak width):", input$tbias
         ))
       )
@@ -431,9 +753,9 @@ textOut <- function(input, optimised=0){
       return(
         renderText(paste(
           "AUTOTRIPLOID MODEL, MANUAL FIT",
-          "\n    haploid k-mer cov:", input$tkcov,
-          "\n      per k-mer theta:", input$tth,
-          "\nhapl non-rep GS (Mbp):", input$tyadj,
+          "\n  monoploid k-mer cov:", input$tkcov,
+          "\n          θ per k-mer:", input$tth,
+          "\n     non-rep GS (Mbp):", input$tyadj,
           "\n    bias (peak width):", input$tbias
         ))
       )
@@ -442,11 +764,12 @@ textOut <- function(input, optimised=0){
       return(
         renderText(paste(
           "ALLOTRIPLOID MODEL, MANUAL FIT",
-          "\n    haploid k-mer cov:", input$tkcov,
-          "\n      per k-mer theta:", input$tth,
+          "\n  monoploid k-mer cov:", input$tkcov,
+          "\n          θ per k-mer:", input$tth,
           "\n                    T:", input$tdiverg,
-          "\nhapl non-rep GS (Mbp):", input$tyadj,
-          "\n    bias (peak width):", input$tbias
+          "\n     non-rep GS (Mbp):", input$tyadj,
+          "\n    bias (peak width):", input$tbias,
+          "\n     diverg per k-mer:", round(input$tth*input$tdiverg, 4)
         ))
       )
     }
@@ -454,11 +777,12 @@ textOut <- function(input, optimised=0){
       return(
         renderText(paste(
           "ALLOTETRAPLOID MODEL, MANUAL FIT",
-          "\n    haploid k-mer cov:", input$tkcov,
-          "\n      per k-mer theta:", input$tth,
+          "\n  monoploid k-mer cov:", input$tkcov,
+          "\n          θ per k-mer:", input$tth,
           "\n                    T:", input$tdiverg,
-          "\nhapl non-rep GS (Mbp):", input$tyadj,
-          "\n    bias (peak width):", input$tbias
+          "\n     non-rep GS (Mbp):", input$tyadj,
+          "\n    bias (peak width):", input$tbias,
+          "\n     diverg per k-mer:", round(input$tth*input$tdiverg, 4)
         ))
       )
     }
@@ -468,17 +792,17 @@ textOut <- function(input, optimised=0){
       return(
         renderText(paste(
           "ALLOTETRAPLOID MODEL, AUTO FITTED",
-          "\n    haploid k-mer cov:", round(optimised$par[1],1),
-          "\n      per k-mer theta:", round(optimised$par[3],4),
+          "\n  monoploid k-mer cov:", round(optimised$par[1],1),
+          "\n          θ per k-mer:", round(optimised$par[3],4),
           "\n                    T:", round(optimised$par[5],2),
-          "\nhapl non-rep GS (Mbp):", round(optimised$par[4],1),
+          "\n     non-rep GS (Mbp):", round(optimised$par[4],1),
           "\n    bias (peak width):", round(optimised$par[2],1),
           "\n     per k-mer diverg:", round(optimised$par[3]*optimised$par[5], 3),
           "\n\nSTARTING RANGES (MIN MAX)",
-          "\n    haploid k-mer cov:", input$akcov[1], input$akcov[2],
-          "\nlog10 per k-mer theta:", input$ath[1], input$ath[2],
+          "\n  monoploid k-mer cov:", input$akcov[1], input$akcov[2],
+          "\n    log10 θ per k-mer:", input$ath[1], input$ath[2],
           "\n                    T:", input$adiv[1], input$adiv[2],
-          "\nhapl non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
+          "\n     non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
           "\n    bias (peak width):", input$abias[1], input$abias[2],
           "\n              x range:", input$axrange[1],input$axrange[2]
         ))
@@ -488,17 +812,17 @@ textOut <- function(input, optimised=0){
       return(
         renderText(paste(
           "ALLOTRIPLOID MODEL, AUTO FITTED",
-          "\n    haploid k-mer cov:", round(optimised$par[1],1),
-          "\n      per k-mer theta:", round(optimised$par[3],4),
+          "\n  monoploid k-mer cov:", round(optimised$par[1],1),
+          "\n          θ per k-mer:", round(optimised$par[3],4),
           "\n                    T:", round(optimised$par[5],2),
-          "\nhapl non-rep GS (Mbp):", round(optimised$par[4],1),
+          "\n     non-rep GS (Mbp):", round(optimised$par[4],1),
           "\n    bias (peak width):", round(optimised$par[2],1),
           "\n     per k-mer diverg:", round(optimised$par[3]*optimised$par[5], 3),
           "\n\nSTARTING RANGES (MIN MAX)",
-          "\n    haploid k-mer cov:", input$akcov[1], input$akcov[2],
-          "\nlog10 per k-mer theta:", input$ath[1], input$ath[2],
+          "\n  monoploid k-mer cov:", input$akcov[1], input$akcov[2],
+          "\n    log10 θ per k-mer:", input$ath[1], input$ath[2],
           "\n                    T:", input$adiv[1], input$adiv[2],
-          "\nhapl non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
+          "\n     non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
           "\n    bias (peak width):", input$abias[1], input$abias[2],
           "\n              x range:", input$axrange[1],input$axrange[2]
         ))
@@ -507,15 +831,15 @@ textOut <- function(input, optimised=0){
     if(input$mod=="tau"){
       return(
         renderText(paste(
-          "AUTOTETRAPLOID MODLE, AUTO FITTED",
-          "\n    haploid k-mer cov:", round(optimised$par[1],1),
-          "\n      per k-mer theta:", round(optimised$par[3],4),
-          "\nhapl non-rep GS (Mbp):", round(optimised$par[4],1),
+          "AUTOTETRAPLOID MODEL, AUTO FITTED",
+          "\n  monoploid k-mer cov:", round(optimised$par[1],1),
+          "\n          θ per k-mer:", round(optimised$par[3],4),
+          "\n     non-rep GS (Mbp):", round(optimised$par[4],1),
           "\n    bias (peak width):", round(optimised$par[2],1),
           "\n\nSTARTING RANGES (MIN MAX)",
-          "\n    haploid k-mer cov:", input$akcov[1], input$akcov[2],
-          "\nlog10 per k-mer theta:", input$ath[1], input$ath[2],
-          "\nhapl non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
+          "\n  monoploid k-mer cov:", input$akcov[1], input$akcov[2],
+          "\n    log10 θ per k-mer:", input$ath[1], input$ath[2],
+          "\n     non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
           "\n    bias (peak width):", input$abias[1], input$abias[2],
           "\n              x range:", input$axrange[1],input$axrange[2]
         ))
@@ -524,15 +848,15 @@ textOut <- function(input, optimised=0){
     if(input$mod=="traaa"){
       return(
         renderText(paste(
-          "AUTOTRIPLOID MODLE, AUTO FITTED",
-          "\n    haploid k-mer cov:", round(optimised$par[1],1),
-          "\n      per k-mer theta:", round(optimised$par[3],4),
-          "\nhapl non-rep GS (Mbp):", round(optimised$par[4],1),
+          "AUTOTRIPLOID MODEL, AUTO FITTED",
+          "\n  monoploid k-mer cov:", round(optimised$par[1],1),
+          "\n          θ per k-mer:", round(optimised$par[3],4),
+          "\n     non-rep GS (Mbp):", round(optimised$par[4],1),
           "\n    bias (peak width):", round(optimised$par[2],1),
           "\n\nSTARTING RANGES (MIN MAX)",
-          "\n    haploid k-mer cov:", input$akcov[1], input$akcov[2],
-          "\nlog10 per k-mer theta:", input$ath[1], input$ath[2],
-          "\nhapl non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
+          "\n  monoploid k-mer cov:", input$akcov[1], input$akcov[2],
+          "\n    log10 θ per k-mer:", input$ath[1], input$ath[2],
+          "\n     non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
           "\n    bias (peak width):", input$abias[1], input$abias[2],
           "\n              x range:", input$axrange[1],input$axrange[2]
         ))
@@ -542,14 +866,14 @@ textOut <- function(input, optimised=0){
       return(
         renderText(paste(
           "DIPLOID MODEL, AUTO FITTED",
-          "\n    haploid k-mer cov:", round(optimised$par[1],1),
-          "\n      per k-mer theta:", round(optimised$par[3],4),
-          "\nhapl non-rep GS (Mbp):", round(optimised$par[4],1),
+          "\n  monoploid k-mer cov:", round(optimised$par[1],1),
+          "\n          θ per k-mer:", round(optimised$par[3],4),
+          "\n     non-rep GS (Mbp):", round(optimised$par[4],1),
           "\n    bias (peak width):", round(optimised$par[2],1),
           "\n\nSTARTING RANGES (MIN MAX)",
-          "\n    haploid k-mer cov:", input$akcov[1], input$akcov[2],
-          "\nlog10 per k-mer theta:", input$ath[1], input$ath[2],
-          "\nhapl non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
+          "\n  monoploid k-mer cov:", input$akcov[1], input$akcov[2],
+          "\n    log10 θ per k-mer:", input$ath[1], input$ath[2],
+          "\n     non-rep GS (Mbp):", input$ayadj[1], input$ayadj[2],
           "\n    bias (peak width):", input$abias[1], input$abias[2],
           "\n              x range:", input$axrange[1],input$axrange[2]
         ))
@@ -579,7 +903,7 @@ prepare.spectrum <- function(spe){
   offs = sp@data[1,1] - 1
   if(offs > 0){
     # print("Padding with zeros to fill in lower multiplicities.")
-    prepDF = data.frame(mult=1:offs, count=rep(0, offs))
+    prepDF = data.frame(mult=1:offs, count=rep(NA_real_, offs))
     sp@data <- rbind(prepDF, sp@data)
   }
   return(sp)
@@ -670,7 +994,7 @@ getStartingVals <- function(input){
   }
 }
 
-doOptimisation <- function(input){
+doOptimisation <- function(input, sp){
   if(input$mod %in% c("tal", "traab")){
 
     return(
diff --git a/Tetmer/data-raw/save_spectra.R b/Tetmer/data-raw/save_spectra.R
index 9b9afcc..e2cc4d1 100644
--- a/Tetmer/data-raw/save_spectra.R
+++ b/Tetmer/data-raw/save_spectra.R
@@ -1,5 +1,5 @@
 library(usethis)
-E030 <- read.spectrum("~/Dropbox/manuscripts/2108_kmers_gs/data/analyseSpectra/spectra/E030full.hist.no0", "E030")
-E028 <- read.spectrum("~/git_repos/shiny-k-mers/data/A0", "E028")
+E030 <- read.spectrum("~/Dropbox/manuscripts/2108_kmers_gs/data/analyseSpectra/spectra/E030full.hist.no0", "E. anglica, E030", 21)
+E028 <- read.spectrum("~/git_repos/shiny-k-mers/data/A0", "E. arctica, E028", 27)
 use_data(E030, overwrite = T)
 use_data(E028, overwrite = T)
diff --git a/Tetmer/data/E028.rda b/Tetmer/data/E028.rda
index b8a0f90..920888f 100644
Binary files a/Tetmer/data/E028.rda and b/Tetmer/data/E028.rda differ
diff --git a/Tetmer/data/E030.rda b/Tetmer/data/E030.rda
index e1728d9..7aabd96 100644
Binary files a/Tetmer/data/E030.rda and b/Tetmer/data/E030.rda differ
diff --git a/Tetmer/man/E028.Rd b/Tetmer/man/E028.Rd
index 529b0b2..4793bcc 100644
--- a/Tetmer/man/E028.Rd
+++ b/Tetmer/man/E028.Rd
@@ -7,10 +7,11 @@
 \format{
 A \code{spectrum} object containing a name and data frame:
 \describe{
-\item{name}{The name E028}
-\item{data}{A dataframe with coulmns \describe{
-\item{mult}{K-mer multiplicity}
-\item{count}{The number of different k-mer at a given multiplicity}
+\item{\code{name}}{The name 'E. arctica, E028'}
+\item{\code{k}}{The k-mer length of 27}
+\item{\code{data}}{A dataframe with coulmns \describe{
+\item{\code{mult}}{K-mer multiplicity}
+\item{\code{count}}{The number of different k-mer at a given multiplicity}
 }
 }
 }
@@ -27,5 +28,6 @@ A k-mer spectrum if allotetraploid \emph{E. arctica}.
 \examples{
 plot(E028, log="xy")
 plot(E028, xlim=c(0, 200), ylim=c(0, 10000000))
+tetmer(E028)
 }
 \keyword{datasets}
diff --git a/Tetmer/man/E030.Rd b/Tetmer/man/E030.Rd
index 058bb0c..bda1d8a 100644
--- a/Tetmer/man/E030.Rd
+++ b/Tetmer/man/E030.Rd
@@ -7,16 +7,17 @@
 \format{
 A \code{spectrum} object containing a name and data frame:
 \describe{
-\item{name}{The name E028}
-\item{data}{A dataframe with coulmns \describe{
-\item{mult}{K-mer multiplicity}
-\item{count}{The number of different k-mer at a given multiplicity}
+\item{\code{name}}{The name 'E. anglica, E030'}
+\item{\code{k}}{The k-mer length of 21}
+\item{\code{data}}{A dataframe with coulmns \describe{
+\item{\code{mult}}{K-mer multiplicity}
+\item{\code{count}}{The number of different k-mer at a given multiplicity}
 }
 }
 }
 }
 \source{
-Becher et al. (2020) \url{https://doi.org/10.1016/j.xplc.2020.100105}
+Source!
 }
 \usage{
 E030
@@ -27,5 +28,6 @@ A k-mer spectrum of diploid \emph{E. anglica}.
 \examples{
 plot(E030, log="xy")
 plot(E030, xlim=c(0, 200), ylim=c(0, 10000000))
+tetmer(E030)
 }
 \keyword{datasets}
diff --git a/Tetmer/man/read.spectrum.Rd b/Tetmer/man/read.spectrum.Rd
index 4e84b4d..9fe04ce 100644
--- a/Tetmer/man/read.spectrum.Rd
+++ b/Tetmer/man/read.spectrum.Rd
@@ -4,12 +4,14 @@
 \alias{read.spectrum}
 \title{Read in a k-mer spectrum}
 \usage{
-read.spectrum(f, nam = "spectrum", ...)
+read.spectrum(f, nam = "MySpectrum", k = 0, ...)
 }
 \arguments{
 \item{f}{A string indicating the path to th spectrum file}
 
-\item{nam}{(otional) A string indicaing the name of the spectrum}
+\item{nam}{(otional) A string indicating the name of the spectrum}
+
+\item{k}{(otional) A numeric indicating the k-mer length}
 
 \item{...}{keyword arguments to be passed to \code{read.table()}}
 }
diff --git a/Tetmer/man/spectrum-class.Rd b/Tetmer/man/spectrum-class.Rd
index 6e912ff..eb59b67 100644
--- a/Tetmer/man/spectrum-class.Rd
+++ b/Tetmer/man/spectrum-class.Rd
@@ -13,5 +13,7 @@ A named k-mer spectrum class
 \item{\code{name}}{A string indicating the name of the spectrum (or sample)}
 
 \item{\code{data}}{A \code{data.frame} with numeric cloumns mult and count.}
+
+\item{\code{k}}{A \code{numeric} indiciating the k-mer length.}
 }}