diff --git a/R/g.analyse.R b/R/g.analyse.R
index e706b6aa9..bc2849490 100644
--- a/R/g.analyse.R
+++ b/R/g.analyse.R
@@ -140,6 +140,7 @@ g.analyse =  function(I,C,M,IMP,qlevels=c(),qwindow=c(0,24),quantiletype = 7,L5M
   HFENplusi = which(colnames(metashort) == "HFENplus")
   MADi = which(colnames(metashort) == "MAD")
   ENi = which(colnames(metashort) == "EN")
+  ENMOai = which(colnames(metashort) == "ENMOa")
   ANYANGLEi = which(colnames(M$metashort) %in% c("anglex","angley","anglez") ==  TRUE)
   if (length(ANYANGLEi) == 0) ANYANGLEi = -1
   if (length(ENMOi) == 0) ENMOi = -1
@@ -149,6 +150,7 @@ g.analyse =  function(I,C,M,IMP,qlevels=c(),qwindow=c(0,24),quantiletype = 7,L5M
   if (length(HFENplusi) == 0) HFENplusi = -1
   if (length(MADi) == 0) MADi = -1
   if (length(ENi) == 0) ENi = -1
+  if (length(ENMOai) == 0) ENMOai = -1
   #===============================================
   # Extract features from the imputed data
   qcheck = r5long
@@ -176,14 +178,13 @@ g.analyse =  function(I,C,M,IMP,qlevels=c(),qwindow=c(0,24),quantiletype = 7,L5M
   qlevels_names = output_avday$qlevels_names
   ML5AD=output_avday$ML5AD
   ML5AD_names = output_avday$ML5AD_names
-
   #--------------------------------------------------------------
   # Analysis per day
   if (doperday == TRUE) {
     output_perday = g.analyse.perday(selectdaysfile, ndays, firstmidnighti, time, nfeatures,
                                      window.summary.size, qwindow, midnightsi, metashort, averageday,
                                      ENMOi, LFENMOi, BFENi, ENi,
-                                     HFENi, HFENplusi, MADi, doiglevels, nfulldays, lastmidnight,
+                                     HFENi, HFENplusi, MADi,  ENMOai, doiglevels, nfulldays, lastmidnight,
                                      ws3, ws2, qcheck, fname, idloc, BodyLocation, wdayname,
                                      tooshort, includedaycrit, winhr,L5M5window, M5L5res,
                                      doquan, qlevels, quantiletype, doilevels, ilevels, iglevels, domvpa,
diff --git a/R/g.analyse.perday.R b/R/g.analyse.perday.R
index bbc4d7d1b..8177e6d04 100644
--- a/R/g.analyse.perday.R
+++ b/R/g.analyse.perday.R
@@ -1,6 +1,6 @@
 g.analyse.perday = function(selectdaysfile, ndays, firstmidnighti, time, nfeatures, 
                             window.summary.size, qwindow, midnightsi, metashort, averageday,
-                            ENMOi, LFENMOi, BFENi, ENi, HFENi, HFENplusi, MADi,
+                            ENMOi, LFENMOi, BFENi, ENi, HFENi, HFENplusi, MADi, ENMOai,
                             doiglevels, nfulldays,lastmidnight, ws3, ws2, qcheck,
                             fname, idloc, BodyLocation, wdayname, tooshort, includedaycrit,
                             winhr, L5M5window, M5L5res,
@@ -300,7 +300,7 @@ g.analyse.perday = function(selectdaysfile, ndays, firstmidnighti, time, nfeatur
               }
               # Starting filling output matrix daysummary with variables per day segment and full day.
               if (mi == ENMOi | mi == LFENMOi | mi == BFENi |
-                  mi == ENi | mi == HFENi | mi == HFENplusi | mi == MADi) {
+                  mi == ENi | mi == HFENi | mi == HFENplusi | mi == MADi | mi == ENMOai) {
                 if (length(varnum) > ((60/ws3)*60*5.5)) { # Calculate values
                   exfi = 0
                   for (winhr_value in winhr) {
diff --git a/R/g.part3.R b/R/g.part3.R
index 5d23dba44..abbf4fb40 100644
--- a/R/g.part3.R
+++ b/R/g.part3.R
@@ -14,7 +14,7 @@ g.part3 = function(metadatadir=c(),f0,f1,anglethreshold = 5,timethreshold = 5,
     dir.create(file.path(paste(metadatadir,"/meta",sep=""),"sleep.qc"))
   }
   #------------------------------------------------------
-  fnames  =dir(paste(metadatadir,"/meta/basic",sep=""))
+  fnames  =dir(paste(metadatadir,"/meta/ms2.out",sep=""))
   if (f1 > length(fnames) | f1 == 0) f1 = length(fnames)
   if (f0 > length(fnames) | f0 == 0) f0 = 1
   #========================================================================
@@ -67,14 +67,14 @@ g.part3 = function(metadatadir=c(),f0,f1,anglethreshold = 5,timethreshold = 5,
                                 .export=functions2passon, .errorhandling=errhand) %myinfix% { 
     tryCatchResult = tryCatch({
       # for (i in f0:f1) {
-      FI = file.info(paste(metadatadir,"/meta/basic/",fnames[i],sep=""))
+      FI = file.info(paste(metadatadir,"/meta/ms2.out/",fnames[i],sep=""))
       if (is.na(FI$size) == TRUE) FI$size = 0
       if (FI$size == 0 | is.na(FI$size) == TRUE | length(FI$size) == 0) {
         cat(paste("P3 file ",fnames[i],sep=""))
         cat("Filename not recognised")
       }
       fname = unlist(strsplit(fnames[i],".RData"))[1]
-      fname = unlist(strsplit(fname,"eta_"))[2]
+      # fname = unlist(strsplit(fname,"eta_"))[2]
       #=========================================================
       #check whether file has already been processed
       #by comparing filename to read with list of processed files
@@ -92,8 +92,8 @@ g.part3 = function(metadatadir=c(),f0,f1,anglethreshold = 5,timethreshold = 5,
         # Load previously stored meta-data from part1.R
         cat(paste(" ",i,sep=""))
         IMP = M = c()
-        load(paste(metadatadir,"/meta/basic/",fnames[i],sep=""))
-        load(paste(metadatadir,"/meta/ms2.out/",unlist(strsplit(fnames[i], "eta_"))[2],sep=""))
+        load(paste(metadatadir,"/meta/basic/meta_",fnames[i],sep=""))
+        load(paste(metadatadir,"/meta/ms2.out/",fnames[i],sep=""))
         if (M$filecorrupt == FALSE & M$filetooshort == FALSE) {
           # IMP = g.impute(M,I,strategy=1,hrs.del.start=0,hrs.del.end=0,maxdur=0)
           SLE = g.sib.det(M,IMP,I,twd=c(-12,12),timethreshold=timethreshold,anglethreshold=anglethreshold,
diff --git a/R/g.part4.R b/R/g.part4.R
index 37059416f..1e28e2a72 100644
--- a/R/g.part4.R
+++ b/R/g.part4.R
@@ -421,6 +421,7 @@ g.part4 = function(datadir=c(),metadatadir=c(),f0=f0,f1=f1,idloc=1,loglocation =
                   calendardate[j] = DD$calendardate
                 }
                 spo = DD$spo
+                reversetime2 = reversetime3 = c()
                 if (daysleeper[j] == TRUE) {
                   if (loaddaysi == 1) {
                     w1 = which(spo[,3] >= 18) #only use periods ending after 6pm
@@ -461,14 +462,12 @@ g.part4 = function(datadir=c(),metadatadir=c(),f0=f0,f1=f1,idloc=1,loglocation =
                     name2 = paste("spo_day2",k,sep="")
                     tmpCmd = paste("spo = rbind(",name1,",",name2,")",sep="")
                     eval(parse(text=tmpCmd))
-                    #reverse back the timestamps to remember that these timestamps were coming from different days
-                    spo[which(spo[,3] >= 36),3] =  spo[which(spo[,3] >= 36),3] - 24
-                    spo[which(spo[,2] >= 36),2] =  spo[which(spo[,2] >= 36),2] - 24
                   }
                 }
                 # spo is now a matrix of onset and wake for each sleep period (episode)
                 for (evi in 1:nrow(spo)) { #Now classify as being part of the SPT window or not
                   if (spo[evi,2] < SptWake & spo[evi,3] > SptOnset) { # = acconset < logwake  & accwake > logonset
+                    
                     spo[evi,4] = 1 #nocturnal = all acc periods that end after diary onset and start before diary wake
                     # REDEFINITION OF ONSET/WAKE OF THIS PERIOD OVERLAPS
                     if (relyonsleeplog == TRUE) { #if TRUE then sleeplog value is assigned to accelerometer-based value for onset and wake up
@@ -481,6 +480,14 @@ g.part4 = function(datadir=c(),metadatadir=c(),f0=f0,f1=f1,idloc=1,loglocation =
                     }
                   }
                 }
+                if (daysleeper[j] == TRUE) {
+                  # for the labelling above it was needed to have times > 36, but for the plotting 
+                  # time in the second day needs to be returned to a normal 24 hour scale. 
+                  reversetime2 = which(spo[,2] >= 36)
+                  reversetime3 = which(spo[,3] >= 36)
+                  if (length(reversetime2) > 0) spo[reversetime2,2] = spo[reversetime2,2] - 24
+                  if (length(reversetime3) > 0) spo[reversetime3,3] =  spo[reversetime3,3] - 24
+                }
                 #------------------------------------------------------------------------
                 # Variable 'spo' contains all the sleep periods FOR ONE SLEEP DEFINITION
                 # Variable 'spocum' contains all the sleep periods FOR MULTIPLE SLEEP DEFINITIONS
@@ -556,6 +563,8 @@ g.part4 = function(datadir=c(),metadatadir=c(),f0=f0,f1=f1,idloc=1,loglocation =
                   
                   #------------------------------------
                   # ACCELEROMETER
+                  if (is.matrix(spocum.t) == FALSE) spocum.t = as.matrix(spocum.t) # seems needed in rare occasions
+                  if (ncol(spocum.t) < 4 & nrow(spocum.t) > 3) spocum.t = t(spocum.t) # seems needed in rare occasions
                   if (length(which(as.numeric(spocum.t[,4]) == 1)) > 0) {
                     rtl = which(spocum.t[,4] == 1)
                     nightsummary[sumi,3] =spocum.t[rtl[1],2]
diff --git a/R/g.part5.R b/R/g.part5.R
index 7bf2b3c39..c690777c6 100644
--- a/R/g.part5.R
+++ b/R/g.part5.R
@@ -39,9 +39,15 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
   fnames.ms3 = sort(dir(paste(metadatadir,"/meta/ms3.out",sep="")))
   fnames.ms4 = sort(dir(paste(metadatadir,"/meta/ms4.out",sep="")))
   fnames.ms5 = sort(dir(paste(metadatadir,"/meta/ms5.out",sep="")))
-  # fnames.ms5raw = sort(dir(paste(metadatadir,"/meta/ms5.outraw",sep="")))
-  # results
+  # path results folder
   results = paste(metadatadir,"/results",sep="")
+  # path to sleeplog milestonedata, if it exists:
+  sleeplogRDA = paste(metadatadir,"/meta/sleeplog.RData",sep="")
+  if (file.exists(sleeplogRDA) == TRUE){
+    load(sleeplogRDA)
+  } else {
+    sleeplog = c()
+  }
   #------------------------------------------------
   # specify parameters
   ffdone = fnames.ms5 #ffdone is now a list of files that have already been processed by g.part5
@@ -52,7 +58,6 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
   di = 1
   cnt = 1
   fnames.ms3 = sort(fnames.ms3)
-  # fnames_original= sort(fnames_original)
   if (f1 == 0) length(fnames.ms4)
   if (f1 > length(fnames.ms4)) f1 = length(fnames.ms4)
   boutdur.mvpa = sort(boutdur.mvpa,decreasing = TRUE)
@@ -89,7 +94,7 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
   if (do.parallel == TRUE) {
     cat(paste0('\n Busy processing ... see ',metadatadir,'/ms5', ' for progress\n'))
   }
-  # check whether we are indevelopment mode:
+  # check whether we are in development mode:
   GGIRinstalled = is.element('GGIR', installed.packages()[,1])
   packages2passon = functions2passon = NULL
   GGIRloaded = "GGIR" %in% .packages()
@@ -107,8 +112,6 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
   output_list =foreach::foreach(i=f0:f1,  .packages = packages2passon, 
                                 .export=functions2passon, .errorhandling=errhand) %myinfix% { # the process can take easily 1 minute per file, so probably there is a time gain by doing it parallel
     tryCatchResult = tryCatch({
-    
-    # for (i in f0:f1) {
     if (length(ffdone) > 0) {
       if (length(which(ffdone == fnames.ms3[i])) > 0) { 
         skip = 1 #skip this file because it was analysed before")
@@ -121,7 +124,11 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
     if (overwrite == TRUE) skip = 0
     # skip files from ms3 if there is no equivalent in ms4
     selp = which(fnames.ms4 == fnames.ms3[i])
-    if (file.exists(paste(metadatadir,"/meta/ms4.out/",fnames.ms4[selp],sep="")) == FALSE) {
+    if (length(selp) > 0) {
+      if (file.exists(paste(metadatadir,"/meta/ms4.out/",fnames.ms4[selp],sep="")) == FALSE) {
+        skip = 1
+      }
+    } else {
       skip = 1
     }
     if (skip == 0) {
@@ -137,8 +144,6 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
       # load output g.part2
       selp = which(fnames.ms2 == fnames.ms3[i]) # so, fnames.ms3[i] is the reference point for filenames
       load(file=paste(metadatadir,"/meta/ms2.out/",fnames.ms2[selp],sep=""))
-      # daysummary = SUM$daysummary # commented out because not used
-      # summary = SUM$summary # commented out  because not used
       # load output g.part4
       selp = which(fnames.ms4 == fnames.ms3[i])
       load(file=paste(metadatadir,"/meta/ms4.out/",fnames.ms4[selp],sep=""))
@@ -154,27 +159,27 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
                             *length(unique(timewindow))),nfeatures)
       di = 1
       fi = 1
+      sptwindow_HDCZA_end = c() # if it is not loaded from part3 milestone data then this will be the default
       if (length(idindex) > 0 & nrow(summarysleep) > 1) { #only attempt to load file if it was processed for sleep
         summarysleep_tmp = summarysleep
         #======================================================================
         # load output g.part1
         selp = which(fnames.ms1 == paste("meta_",fnames.ms3[i],sep=""))
         if (length(selp) != 1) {
-          cat("Warning: File not processed with part1")
+          cat("Warning: Milestone data part 1 could not be retrieved")
         }
         load(paste0(metadatadir,"/meta/basic/",fnames.ms1[selp]))
         # load output g.part3
         load(paste0(metadatadir,"/meta/ms3.out/",fnames.ms3[i]))
         # extract key variables from the mile-stone data: time, acceleration and elevation angle
-        # IMP = g.impute(M,I,strategy=strategy,hrs.del.start=hrs.del.start,
-        #                hrs.del.end=hrs.del.end,maxdur=maxdur)
         time = IMP$metashort[,1]
-        ACC = IMP$metashort[,acc.metric] * 1000 #note that this is imputed ACCELERATION because we use this for describing behaviour
-        
+        #note that this is imputed ACCELERATION because we use this for describing behaviour:
+        ACC = IMP$metashort[,acc.metric] * 1000 
         if (length(which(names(IMP$metashort) == "anglez")) == 0) {
           cat("Warning: anglez not extracted. Please check that do.anglez == TRUE")
         }
-        angle = as.numeric(as.matrix(M$metashort[,which(names(IMP$metashort) == "anglez")])) #note that this is the non-imputed angle because we use this here only for the visualisation
+        #note that this is the non-imputed angle because we use this here only for the visualisation:
+        angle = as.numeric(as.matrix(M$metashort[,which(names(IMP$metashort) == "anglez")])) 
         nonwear = IMP$rout[,5]
         nonwear = rep(nonwear,each=(IMP$windowsizes[2]/IMP$windowsizes[1]))
         if (length(nonwear) > length(time)) {
@@ -193,18 +198,23 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
         for (j in def) { # loop through sleep definitions (defined by angle and time threshold in g.part3)        
           #========================================================
           # SUSTAINED INACTIVITY BOUTS
-          Sbackup = S
+          # These are stored in part 3 milestone data as start- and end-times
+          # in the following code we convert the them into indices of the recording sequence
           S2 = S[which(S$definition==j),] # simplify to one definition
-          sibdetection = rep(0,length(time))
-          tmp = S2 #[which(S2$night==h),]
+          sibdetection = rep(0,length(time)) # initialize output vector that will hold the sibs
           s0s1 = c()
+          # pr0 and pr1 define the indices relative the start of the recording
+          # and specifying a 6 day time window 
+          # always relative to the most recently processed sustained inactivity bout
+          # s0 and s1 are the indices within that time window 
+          # that match the start and end of the next sustained inactivity bout 
           pr0 = 1
           pr1 = pr0 + ((60/ws3)*1440*6)
           pr2 = length(time)
           if (nrow(S2) > 0) {
             gik.ons = as.character(S2$sib.onset.time)
             gik.end = as.character(S2$sib.end.time)
-            for (g in 1:nrow(S2)) { # sleep periods
+            for (g in 1:nrow(S2)) { # sustained inactivity bouts
               lastpr0 = pr0
               pr1 = pr0 + ((60/ws3)*1440*6)
               if (pr1 > pr2) pr1 = pr2
@@ -216,15 +226,14 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
               }
               s0 = which(timebb == gik.ons[g])[1]
               s1 = which(timebb == gik.end[g])[1]
-              #Change 1 - 03/10/2017
               if ( timebb[1] != as.character(timebb[1])){ #not s0 because s0 does not exist yet if classes differ
                 timebb = as.character(timebb)
                 s0 = which(timebb == gik.ons[g])[1]
                 s1 = which(timebb == gik.end[g])[1]
               }
-              #End of change 1
               if (is.na(s0) == TRUE) s0 = which(timebb == paste(gik.ons[g]," 00:00:00",sep=""))[1]
               if (is.na(s1) == TRUE) s1 = which(timebb == paste(gik.end[g]," 00:00:00",sep=""))[1]
+              # add pr0 to make s0 and s1 be relative to the start of the recording
               s0 = s0 + pr0 - 1
               s1 = s1 + pr0 - 1
               pr0 = s1
@@ -254,14 +263,17 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
           # following code was move to here, because otherwise it would repeated remove the last night in the loop          
           #           ignore last night because the following day is potentially not complete e.g. by reduce protocol compliance
           #=========
-          # SUSTAINED INACTIVITY BOUTS are normally not assessed for the time between the first midnight and the first noon
-          # because it is not used for sleep reports (first night is ignored)
+          # If excludefirstlast was set to TRUE in part 4 then 
+          # SUSTAINED INACTIVITY BOUTS are not assessed for the time between
+          # the first midnight and the first noon
+          # because it is not used for sleep reports (first night is ignored in this case).
           # Therefore, expand the SIB sibdetection to include this time segment.
+          #----------------------------------------------------------------------
           # Step 1: Identify time window that needs to be processed
-          redo1 = nightsi[1] - ((60/ws3)*60)
+          redo1 = nightsi[1] - ((60/ws3)*60) # 1 hour before first midnight
           if (redo1 < 1) redo1 = 1
-          redo2 = nightsi[1] + (14*(60/ws3)*60)
-          # specify defintion of sustained inactivity bout    
+          redo2 = nightsi[1] + (14*(60/ws3)*60) # 14 hours after first midngiht
+          # Specify defintion of sustained inactivity bout
           anglethreshold = as.numeric(unlist(strsplit(j,"A"))[2])
           tempi = unlist(strsplit(unlist(strsplit(j,"A"))[1],"T"))
           timethreshold = as.numeric(tempi[length(tempi)])
@@ -298,7 +310,7 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
           }
           sibdetection[redo1:redo2] = sdl1
           #========================================================
-          # DIURNAL BINARY CLASSIFICATION INTO SLEEP OR WAKE PERIOD 
+          # DIURNAL BINARY CLASSIFICATION INTO NIGHT (onset-wake) OR DAY (wake-onset) PERIOD 
           # Note that the sleep date timestamp corresponds to day before night
           w0 = w1 = rep(0,length(summarysleep_tmp2$calendardate))
           diur = rep(0,length(time)) #diur is a variable to indicate the diurnal rhythm: 0 if wake/daytime, 1 if sleep/nighttime
@@ -307,11 +319,12 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
             summarysleep_tmp2 = summarysleep_tmp2[-pko,]
           }
           if (nrow(summarysleep_tmp2) > 0) {
-            for (k in 1:length(summarysleep_tmp2$calendardate)){
+            for (k in 1:length(summarysleep_tmp2$calendardate)){ # loop through nights from part 4
+              # combine date (the tt object) and time (inside the _ts columns of summarysleep_tmp2)
               tt = unlist(strsplit(as.character(summarysleep_tmp2$calendardate[k]),"/"))
               w0[k] = paste(tt[3],"-",tt[2],"-",tt[1]," ",as.character(summarysleep_tmp2$acc_onset_ts[k]),sep="")
               w1[k] = paste(tt[3],"-",tt[2],"-",tt[1]," ",as.character(summarysleep_tmp2$acc_wake_ts[k]),sep="")
-              if (summarysleep_tmp2$acc_onset[k] >= 24) {
+              if (summarysleep_tmp2$acc_onset[k] >= 24) { # if time is beyond 24 then change the date
                 # I am making the hard assumption here that a day has 24 hours...future improvement needed
                 w0[k] = as.character(as.POSIXlt(w0[k],tz=desiredtz) + (24*3600))
               }
@@ -323,7 +336,6 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
               w1c = as.character(as.POSIXlt(w1[k],tz=desiredtz))
               s0 = which(as.character(time) == w0c)[1]
               s1 = which(as.character(time) == w1c)[1]
-              
               if (length(s0) == 0) {
                 w0c = paste0(w0c," 00:00:00")
                 s0 = which(as.character(time) == w0c)[1]
@@ -363,24 +375,45 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
                 diur[s0:s1] = 1
               }
             }
-            
-            if (diur[1] != 0 & (60/ws3)*60 < nightsi[1]) { # the statement with unique that was here results is ambiguous because it can have multiple values
-              # waketi = which(diff(detection) == -1)[1]
-              # onseti = which(diff(detection) == 1)[1] # added 20/4/2018 by VvH
-              waketi = which(diff(diur) == -1)[1]
-              onseti = which(diff(diur) == 1)[1] # added 20/4/2018 by VvH
-            } else {
-              waketi = 0
-              onseti = 0 # added 20/4/2018 by VvH
-            }
-            if (waketi > 0 & onseti == 0) {
-              diur[1:waketi] = 1
-            } else if (waketi > 0 & onseti > 0) { # added 20/4/2018 by VvH
-              diur[onseti:waketi] = 1
+            # Note related to if first and last night were ignored in part 4:
+            # - diur lack sthe first and last night at this point in the code.
+            # - nightsi has all the midnights, so it is possible to check here
+            # whether a wakeup time is missing on the first full day.
+            # - if it is missing, then we will impute it in order for part5 to
+            # the wake-to-wake analys on the second recording day.
+            # Previously we accounted only for this later on in the code, whcih
+            # did not benefit the experoted timeseries.
+            firstwake = which(diff(diur) == -1)[1]
+            if (firstwake > nightsi[2]) { # test whether wake for second day is missing
+              Nepochsinhour = (60/ws3) * 60
+              if (length(sleeplog) > 0) {
+                # use sleeplog for waking up after first night
+                wake_night1 = sleeplog$sleepwake[which(sleeplog$id == id & sleeplog$night == 1)]
+                wake_night1_index =c()
+                if (length(wake_night1) != 0) {
+                  wake_night1_num = as.numeric(unlist(strsplit(wake_night1,":"))) / c(1,60,3600)
+                  wake_night1_hour = sum(wake_night1_num)
+                  wake_night1_index = nightsi[1] + round(wake_night1_hour * Nepochsinhour)
+                } else { # use HDCZA algorithm as plan B
+                  wake_night1_index = round(sptwindow_HDCZA_end[1] * Nepochsinhour)
+                }
+              }
+              if (length(sptwindow_HDCZA_end) > 0 & length(sleeplog) == 0) {
+                # use HDCZA algortihm for waking up after first night
+                # if there was no sleep log
+                wake_night1_index = round(sptwindow_HDCZA_end[1] * Nepochsinhour)
+              }
+              if (length(wake_night1_index) == 0) {
+                # use waking up from next day and subtract 24 hours,
+                # the final option if neither of the above routes works
+                wake_night1_index = (firstwake - (24* ((60/ws3)*60))) + 1
+              }
+              diur[nightsi[1]:wake_night1_index] = 1
             }
             for (TRLi in threshold.lig) {
               for (TRMi in threshold.mod) {
                 for (TRVi in threshold.vig) {
+                  # derive behavioral levels (class), e.g. MVPA, inactivity bouts, etc.
                   levels = identify_levels(time,diur,sibdetection,ACC,
                                            TRLi,TRMi,TRVi,
                                            boutdur.mvpa,boutcriter.mvpa,
@@ -425,7 +458,8 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
                         nightsi = nightsi[which(nightsi > FM[1] & nightsi < FM[length(FM)])]
                       }
                     } else {
-                      startend_sleep = which(abs(diff(diur))==1)  # newly added on 31-3-2019, because if first night is missing then nights needs to allign with diur
+                      # newly added on 31-3-2019, because if first night is missing then nights needs to allign with diur
+                      startend_sleep = which(abs(diff(diur))==1) 
                       Nepochsin12Hours =  (60/ws3)*60*12
                       nightsi = nightsi[which(nightsi >= (startend_sleep[1] - Nepochsin12Hours) &
                                               nightsi <= (startend_sleep[length(startend_sleep)] + Nepochsin12Hours))]  # newly added on 25-11-2019
@@ -433,12 +467,11 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
                       plusrow = 1
                     }
                     for (wi in 1:(nrow(summarysleep_tmp2)+plusrow)) { #loop through 7 windows (+1 to include the data after last awakening)
-                      #check that this is a meaningful day
+                      # Check that this is a meaningful day (that is all the qqq variable is used for),
+                      # before storing it.
                       qqq = rep(0,2)
-                      # check that it is possible to find both windows in the data for this day
-                      # qqq definitions changed so we get acc_onset and acc_wake regarding to the same day of measurement in the same row. 
-                      # Also, it allows for the analysis of the first day for those studies in which the accelerometer is 
-                      # started during the morning and the first day is of interest.
+                      # Check that it is possible to find both windows (WW and MM)
+                      # in the data for this day.
                       if (timewindowi == "MM") {
                         if (wi==1) {
                           qqq[1] = 1
@@ -446,29 +479,27 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
                         } else if (wi<=nrow(summarysleep_tmp2)) {
                           qqq[1] = nightsi[wi-1] + 1
                           qqq[2] = nightsi[wi]
+                          qqq_backup = qqq
                         } else if (wi>nrow(summarysleep_tmp2)) {
-                          qqq[1] = qqq[2] + 1
-                          #qqq[2] = length(diur)
+                          qqq[1] = qqq_backup[2] + 1
                           qqq[2] = nightsi[wi]
                           if(is.na(qqq[2])==TRUE | length(diur) < qqq[2]) qqq[2] = length(diur)
                         }
                       } else if(timewindowi == "WW") {
-                        if (wi==1){
-                          if (waketi > 0) {
-                            qqq[1] = waketi +1
-                          } else {
-                            qqq[1] = 1
-                          }
+                        if (wi==1){ # start measurement till first waking up
+                          qqq[1] = 1
                           qqq[2]=which(diff(diur) == -1)[wi] + 1
-                        } else if (wi<=nrow(summarysleep_tmp2)) {
+                        } else if (wi<=nrow(summarysleep_tmp2)) { # all full wake to wake days
                           qqq[1] = which(diff(diur) == -1)[wi-1] + 1
                           qqq[2] = which(diff(diur) == -1)[wi]
                         } else if (wi>nrow(summarysleep_tmp2)) {
-                          qqq[1] = qqq[2] + 1
-                          qqq[2] = length(diur)
+                          # time after last reliable waking up (this can be more than 24 hours)
+                          # ignore this day, because if the night was ignored for sleep analysis
+                          # then the description of the day in part 5 including that night is
+                          # not informative.
+                          qqq = c(NA, NA)
                         }
                       }
-                      
                       if (length(which(is.na(qqq)==TRUE)) == 0) { #if it is a meaningful day then none of the values in qqq should be NA
                         fi = 1
                         # START STORING BASIC INFORMATION
@@ -496,14 +527,17 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
                                                    "daysleeper","cleaningcode","sleeplog_used","acc_available"); fi = fi + 12
                           
                         } else {
-                          if(timewindowi == "MM" & wi == nrow(summarysleep_tmp2)+plusb) { # for the last day a ot of information is missing, so fill in defaults 
-                            dsummary[di,fi:(fi+1)] = c(weekdays(as.Date(summarysleep_tmp2$calendardate[wi-1], format="%e/%m/%Y")), #+1
-                                                       as.character(as.Date(summarysleep_tmp2$calendardate[wi-1], format="%e/%m/%Y"))) #+1
+                          if(timewindowi == "MM" & wi == nrow(summarysleep_tmp2)+plusb) { 
+                            # for the last day a ot of information is missing, so fill in defaults 
+                            options(encoding = "UTF-8")
+                            Sys.setlocale("LC_TIME", "C") # set language to English
+                            recDate = as.Date(summarysleep_tmp2$calendardate[wi-1], format="%e/%m/%Y")
+                            dsummary[di,fi:(fi+1)] = c(weekdays(recDate), as.character(recDate))
                             addone = 0
                             if (dsummary[di,fi] == remember_previous_weekday) {
                               addone = 1
-                              dsummary[di,fi:(fi+1)] = c(weekdays(as.Date(summarysleep_tmp2$calendardate[wi-1], format="%e/%m/%Y")+1), #+1
-                                                         as.character(as.Date(summarysleep_tmp2$calendardate[wi-1], format="%e/%m/%Y")+1)) #+1
+                              recDate = as.Date(summarysleep_tmp2$calendardate[wi-1], format="%e/%m/%Y") + 1
+                              dsummary[di,fi:(fi+1)] = c(weekdays(recDate), as.character(recDate))
                             }
                             ds_names[fi:(fi+1)] = c("weekday","calendardate");  fi = fi + 2
                             dsummary[di,fi] = j
@@ -570,44 +604,21 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
                           ds_names[fi] = "acc_available";      fi = fi + 1
                         }
                         # 
-                        # define time windows
-                        # qqq definitions changed so we get acc_onset and acc_wake regarding to the same day of measurement in the same row.
+                        # define time windows:
+                        # We will get acc_onset and acc_wake 
+                        # regarding to the same day of measurement in the same row.
+                        # which differs between MM and WW
                         # Also, it allows for the analysis of the first day for those studies in which the accelerometer is started during the morning and the first day is of interest.
+                        # qqq1 is the start of the day
+                        # qqq2 is the end of the day
+                        qqq1 = qqq[1] # added 26 Feb 2020
+                        qqq2 = qqq[2] # added 26 Feb 2020
+                       
                         if (timewindowi == "MM") {
-                          if (wi==1) {
-                            # if (waketi > 0) {
-                            #   qqq1 = waketi + 1
-                            # } else {
-                            qqq1 = 1
-                            # }
-                            qqq2 = nightsi[wi]
-                          } else if (wi<=nrow(summarysleep_tmp2)) {
-                            qqq1 = nightsi[wi-1] + 1
-                            qqq2 = nightsi[wi]
-                          } else if (wi>nrow(summarysleep_tmp2)) {
-                            qqq1 = qqq2 + 1
-                            #qqq2 = length(diur)
-                            qqq2 = nightsi[wi]
-                            if(is.na(qqq2)==TRUE | length(diur) < qqq2) qqq2 = length(diur)
-                          }
                           dsummary[di, fi] = "MM"
                         } else if (timewindowi == "WW") {
-                          if (wi==1){
-                            if (waketi > 0) qqq1 = waketi + 1
-                            else {qqq1 = 1}
-                            qqq2 = which(diff(diur) == -1)[wi] + 1
-                          } else if (wi<=nrow(summarysleep_tmp2)) {
-                            qqq1 = which(diff(diur) == -1)[wi-1] + 1
-                            qqq2 = which(diff(diur) == -1)[wi]
-                          } else if (wi>nrow(summarysleep_tmp2)) {
-                            qqq1 = qqq2 + 1
-                            # qqq2 = length(diur)
-                            qqq2 = nightsi[wi]
-                            if(is.na(qqq2)==TRUE | length(diur) < qqq2) qqq2 = length(diur)
-                          }
                           dsummary[di, fi] = "WW"
                         }
-                        
                         ds_names[fi] = "window";      fi = fi + 1    
                         # keep track of threshold value
                         dsummary[di,fi] = TRLi
@@ -618,13 +629,13 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
                         ds_names[fi] = "TRVi";      fi = fi + 1
                         wlih = ((qqq2-qqq1)+1)/((60/ws3)*60)
                         if (qqq1 > length(LEVELS)) qqq1 = length(LEVELS)
-                        if (wlih > 30 & length(summarysleep_tmp2$night) > 1) { # scenario when day is missing and code reaches out to two days before this day
-                          # if (summarysleep_tmp2$night[wi] - summarysleep_tmp2$night[wi-1] != 1) {
-                          qqq1 = (qqq2 - (24* ((60/ws3)*60))) + 1 # code now uses only 24hours before waking up
-                          if (qqq1 < 1) qqq1 = 1
-                          wlih = ((qqq2-qqq1)+1)/((60/ws3)*60)
-                          # }
-                        }
+                        # # This part should be redundant now, so commented out (26-Feb 2020):
+                        # if (wlih > 30 & length(summarysleep_tmp2$night) > 1) { 
+                        #   # scenario when day is missing and code reaches out to two days before this day
+                        #   qqq1 = (qqq2 - (24* ((60/ws3)*60))) + 1 # code now uses only 24hours before waking up
+                        #   if (qqq1 < 1) qqq1 = 1
+                        #   wlih = ((qqq2-qqq1)+1)/((60/ws3)*60)
+                        # }
                         dsummary[di,fi] = wlih
                         ds_names[fi] = "window_length_in_hours";      fi = fi + 1
                         dsummary[di,fi] = (length(which(nonwear[qqq1:qqq2] == 1)) * ws3) / 3600
@@ -887,12 +898,13 @@ g.part5 = function(datadir=c(),metadatadir=c(),f0=c(),f1=c(),strategy=1,maxdur=7
           if (length(emptycols) > 0) emptycols = emptycols[which(emptycols > lastcolumn)]
           if (length(emptycols) > 0) output = output[-emptycols]
         }
-        save(output,file=paste(metadatadir,ms5.out,"/",fnames.ms3[i],sep=""))
+        if (length(output) > 0) {
+          save(output,file=paste(metadatadir,ms5.out,"/",fnames.ms3[i],sep=""))
+        }
         rm(output,dsummary)
       }
     }
     }) # END tryCatch
-    
     return(tryCatchResult) 
   }
   if (do.parallel == TRUE) {
diff --git a/R/g.plot.R b/R/g.plot.R
index 69f083854..f895bc1a9 100644
--- a/R/g.plot.R
+++ b/R/g.plot.R
@@ -32,7 +32,8 @@ g.plot = function(IMP,M,I,durplot) {
   if (length(fname2) > 1) fname = paste(fname2[1],fname2[2],sep="")
   plot.new()	
   par(fig=c(0,1,0,1),new=T)
-  plot(seq(0,durplot),seq(0,durplot),col="white",type="l",axes=F,xlab="",ylab="",main=paste("device: ",monn," | filename: ",fname,sep=""),cex.main=0.6)#dummy plot
+  plot(seq(0,durplot),seq(0,durplot),col="white",type="l",axes=F, 
+       xlab="",ylab="",main=paste("device: ",monn," | filename: ",fname,sep=""),cex.main=0.6)#dummy plot
   lim = par("usr")
   Y0 = c(-1)
   Y1 = c(durplot+1)
@@ -60,7 +61,9 @@ g.plot = function(IMP,M,I,durplot) {
   if (length(which(colnames(M$metashort) == "ENMO")) > 0) {
     accel = as.numeric(as.matrix(M$metashort[,which(colnames(M$metashort) == "ENMO")]))
   } else {
-    accel = as.numeric(as.matrix(M$metashort[,which(colnames(M$metashort) == "BFEN")]))
+    # Index of first acceleration metric in metashort other than ENMO, timestamp, or angle-related
+    IndOtEnmo = which(colnames(M$metashort) %in% c("ENMO","timestamp","anglex","angley","anglez") == FALSE)[1]
+    accel = as.numeric(as.matrix(M$metashort[,IndOtEnmo]))
   }
   accel2 =cumsum(c(0,accel))
   select = seq(1,length(accel2),by=ws2/ws3)
diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd
index 5d0b40a6d..60731138f 100755
--- a/inst/NEWS.Rd
+++ b/inst/NEWS.Rd
@@ -3,6 +3,12 @@
 \newcommand{\cpkg}{\href{http://CRAN.R-project.org/package=#1}{\pkg{#1}}}
 \section{Changes in version 1.11-1 (GitHub only, release date:28-02-2020)}{
 \itemize{
+  \item Metric ENMOa now facilitated for MVPA calculation
+  \item Bug fixed with part4 daysleeper handling
+  \item Part5 WW window calculation improved, first day now uses sleeplog or
+  HDCZA algorithm estimate, and last day is ignored if no sleep estimates are
+  available. This also affects csv exports by argument save_ms5rawlevels.
+  \item Added explanation to vignette on how to use published cut-points.
   \item Axivity AX6 (acc + gyro) in cwa format now supported for file reading,
   actually using the gyro data for feature calculation is future work.
   In the mean time, gyro signal will be ignored by the rest of GGIR.
diff --git a/man/g.analyse.perday.Rd b/man/g.analyse.perday.Rd
index bb80295ff..274ce66bd 100644
--- a/man/g.analyse.perday.Rd
+++ b/man/g.analyse.perday.Rd
@@ -11,7 +11,7 @@ output matrix, \link{g.analyse}.
 g.analyse.perday(selectdaysfile, ndays, firstmidnighti, 
   time, nfeatures, window.summary.size, qwindow, midnightsi,
   metashort, averageday, ENMOi, LFENMOi, BFENi, ENi,
-  HFENi, HFENplusi, MADi, doiglevels, nfulldays, lastmidnight, 
+  HFENi, HFENplusi, MADi, ENMOai, doiglevels, nfulldays, lastmidnight, 
   ws3, ws2, qcheck, fname, idloc, BodyLocation, wdayname,
   tooshort, includedaycrit, winhr,L5M5window, M5L5res, 
   doquan, qlevels, quantiletype, doilevels, 
@@ -40,6 +40,7 @@ need to be stored in the output matrix}
 \item{HFENi}{column index of metahosrt where metric is stored}
 \item{HFENplusi}{column index of metahosrt where metric is stored}
 \item{MADi}{column index of metahosrt where metric is stored}
+\item{ENMOai}{column index of metahosrt where metric is stored}
 \item{doiglevels}{Boolean to indicate whether iglevels should be calculated}
 \item{nfulldays}{Number of days between the first and last midnight in the recording}
 \item{lastmidnight}{see \link{g.detecmidnight}}
diff --git a/tests/testthat/test_chainof5parts.R b/tests/testthat/test_chainof5parts.R
index c0a63465b..aba49cf1b 100644
--- a/tests/testthat/test_chainof5parts.R
+++ b/tests/testthat/test_chainof5parts.R
@@ -118,7 +118,7 @@ test_that("chainof5parts", {
   load(rn[1])
   expect_true(dir.exists(dirname))
   expect_true(file.exists(rn[1]))
-  expect_that(nrow(output),equals(6)) # changed because part5 now gives also first and last day
+  expect_that(nrow(output),equals(5)) # changed because part5 now gives also first and last day
   expect_that(ncol(output),equals(134))
   expect_that(round(as.numeric(output$acc_wake[2]),digits=4),equals(35.9958))
 
@@ -196,7 +196,7 @@ test_that("chainof5parts", {
   expect_true(dir.exists(dirname))
   expect_true(file.exists(vis_sleep_file))
   expect_that(round(nightsummary$acc_dur_sibd[1],digits=4),equals(0))
-  expect_that(round(nightsummary$acc_SptDuration[1],digits=2),equals(24))
+  expect_that(round(nightsummary$acc_SptDuration[1],digits=2),equals(7))
   expect_true(as.logical(nightsummary$acc_available[1]))
   expect_false(as.logical(nightsummary$sleeplog_used[1]))
   #----------------------------------------------------------------------
diff --git a/vignettes/GGIR.Rmd b/vignettes/GGIR.Rmd
index 86dd1f703..ff3db97fa 100644
--- a/vignettes/GGIR.Rmd
+++ b/vignettes/GGIR.Rmd
@@ -1,7 +1,7 @@
 ---
 title: "Accelerometer data processing with GGIR"
 author: "Vincent van Hees"
-date: "September 6 2019"
+date: "February 28 2020"
 output:
   rmarkdown::html_vignette:
     toc : true
@@ -39,7 +39,7 @@ R package GGIR would not have been possible without the support of the following
 
 - Dr. Vincent van Hees as main developer & maintainer.
 - Jairo Migueles has helped to improve sections of the code.
-- Dr. Evgeny Mirkes, Dr. Jing Hua Zhao, and Dr. Dan Jackson provided code to enable binary accelerometer data to be read in R.
+- Dr. Evgeny Mirkes, Dr. Dan Jackson, and Dr. Jing Hua Zhao contributed to the binary data reading functionality.
 - Dr. Zhou Fang contributed to a part of the code for the auto-calibration functionality.
 - Joe Heywood contributed code to correctly extract specific time segments with correct timestamps from an GENEActiv recording.
 - Joan Capdevila Pujol contributed with a couple of bug fixes.
@@ -49,6 +49,7 @@ R package GGIR would not have been possible without the support of the following
 
 
 
+
 **Institutions:**
 
 - The MRC Epidemiology Unit at the University of Cambridge, Cambridge, UK
@@ -56,6 +57,7 @@ R package GGIR would not have been possible without the support of the following
 - Netherlands eScience Center, Amsterdam, NL
 - University of Exeter, Exeter, UK
 - Centre for Longitudinal Studies, London, UK
+- Vincent van Hees as self-employed consultant
 
 **Funders:**
 
@@ -175,6 +177,61 @@ knitr::include_graphics("sleeplogexample.jpg")
 - `timewindow` to specify whether days should be defined from midnight to midnight, from waking-up to waking-up, or both (all output is generated twice).
 - Configure durations of bouts: `boutdur.mvpa`, `boutdur.in`, and `boutdur.lig`. Note that this can be a vector of multiple values indicating the minimum and maximum duration of subsequent bout types, e.g. 1-5 minutes MVPA, 5-10 minutes MVPA, and longer than 10 minutes MVPA.
 
+### Using published cut-points
+
+Cut-points to estimate time spent in acceleration levels that are roughly liked to levels of energy metabolism have been proposed by:
+
+- Esliger et al 2011: wrist and waist in adults. doi: 10.1249/MSS.0b013e31820513be.
+- Schaefer et al 2014: wrist in 6-11 year old children. doi: 10.1249/MSS.0000000000000150.
+- Roscoe et al 2017: wrist in 4-5 year old pre-school children. doi: 10.1007/s00431-017-2948-2.
+- Phillips et al 2013: wrist and hip in 8-14 year olds. doi: 10.1016/j.jsams.2012.05.013.
+- Vaha-Ypya et al 2015: hip in adults. 10.1371/journal.pone.0134813
+- Hildebrand et al 2014 and 2016: wrist and hip in 7-11 and 21-61 years old. in doi: 10.1111/sms.12795, 10.1249/MSS.0000000000000289.
+- Please let me know if you think I missed a publication.
+
+The first four publications make use of acceleration metrics that sum their values per epoch rather than average them per epoch like GGIR does. However, no need to worry: It is still possible to use those cut-points with a simple cut-point scaling trick and by selecting the right acceleration metric as detailed below.
+
+**Esliger 2011, Phillips 2013:**
+
+-	In GGIR use metric ENMOa instead of ENMO with arguments do.enmoa = TRUE, do.enmo = FALSE, and acc.metric=”ENMOa”.
+- `threshold.lig = (LightCutPointFromPaper/80) * 1000`
+- `threshold.mod = (ModerateCutPointFromPaper/80) * 1000`
+- `threshold.vig = (VigorousCutPointFromPaper/80) * 1000`
+- `mvpathreshold = (ModerateCutPointFromPaper/80) * 1000`
+-	In the part2 results you will need the MVPA estimates that are related to ENMOa, not ENMO.
+-	In the part 5 results everything will be based on the new cut-points.
+
+**Roscoe 2017**
+
+-	In GGIR use metric ENMOa instead of ENMO with arguments do.enmoa = TRUE, do.enmo = FALSE, and acc.metric=”ENMOa”.
+- `threshold.lig = (LightCutPointFromPaper/85.7) * 1000`
+- `threshold.mod = (ModerateCutPointFromPaper/85.7) * 1000`
+- `threshold.vig = (VigorousCutPointFromPaper/85.7) * 1000`
+- `mvpathreshold = (ModerateCutPointFromPaper/85.7) * 1000`
+-	In the part2 results you will need the MVPA estimates that are related to ENMOa, not ENMO.
+-	In the part 5 results everything will be based on the new cut-points.
+
+**Schaeffer 2014:**
+
+-	In GGIR use metric EN instead of ENMO with arguments do.en = TRUE, do.enmo = FALSE, and acc.metric=”EN”.
+-	Specify Schaeffer cut-points as:
+- `threshold.lig = (LightCutPointFromPaper/75) * 1000`
+- `threshold.mod = (ModerateCutPointFromPaper/75) * 1000`
+- `threshold.vig = (VigorousCutPointFromPaper/75) * 1000`
+- `mvpathreshold = (ModerateCutPointFromPaper/75) * 1000`-	In the part2 results you will need the MVPA estimates that are related to EN, not ENMO.
+-	In the part 5 results everything will be based on the new cut-points.
+
+**Vaha-Ypya et al 2015:**
+
+-	Use default setting do.mad= TRUE, acc.metric=”MAD”
+-	Use the cut-points as provided by Vaha-Ypya directly. No need for scaling.
+
+**Hildebrand 2014 and Hildebrand 2016:**
+
+-	Use default setting do.enmo= TRUE, acc.metric=”ENMO”
+-	Use the cut-points as provided by Hildebrand directly. No need for scaling.
+
+
 ### Example call
 
 So, if you consider all the arguments above you me may end up with a call to g.shell.GGIR that could look as follows.