setwd(open_dir)

##
what_you_specified = list.files(pattern = starthere)
for (x in what_you_specified[1])
{rawtable = read.csv(x, colClasses="character", header=T)}
rawtable  = data.frame(rawtable)

################################################################################
##LOOK HERE, list your choise of PBL correction
##Default is the 47 interpolation.

PBL_47_interpo = TRUE
PBL_49_scaling = FALSE
#PBL_47_interpo = FALSE

################################################################################
##Defines required constants for all the calculation, this session is largely###
##imported from the Matlab code presented in Huntington et al., 2009.
OrigPDBCO2_R13  = 0.0112372000
OrigPDBCO2_R18	= 0.0020790000
OrigPDBCO2_R17	= 0.0003799500
VSMOW_R18		= 0.0020052000		#from Gonfianti et al, 1993
VSMOW_R17		= 0.0003799000		#from Gonfianti et al, 1993
VPDBcalcite_R13	= 0.0112372000
VPDBcalcite_R18	= (30.9e-3+1)*VSMOW_R18
lambda		= 0.5164			#from Gonfianti et al, 1993
VPDBcalcite_R17	= VSMOW_R17*(VPDBcalcite_R18/VSMOW_R18)^lambda
VPDBCO2_R13		= VPDBcalcite_R13
VPDBCO2_R18		= VPDBcalcite_R18*1.01025	#isotopic fractionation of
#calcite at 25C (from Kim & Oneil
#1997, from Hoefs 1997)     
VPDBCO2_R17 = VSMOW_R17*(VPDBCO2_R18/VSMOW_R18)^lambda


################################################################################
##START OF THE CALCULATION######################################################
##Specify number of off-peak PBL cycles before and after the on-peak cycles#####

pre_PBL <- 4
post_PBL <- 5

##Here, "row" index is extracted from the "Row" column in isodat data export####
minrow	= min(as.numeric(rawtable$Row))
maxrow	= max(as.numeric(rawtable$Row))
nmbrow	= maxrow + 1 - minrow

##"cyclecount" is a numeric vector that archives the total number of half ######
##cycles number (before formatting) for each acquisition.
cyclecount	= vector(mode = "numeric", length = nmbrow)
for(n in minrow:maxrow)
{
  cyclecount[n - minrow + 1] = nrow(rawtable[which(rawtable$Row == n),])
}

##"cycleindex" is a numeric vector of the FULL CYCLES (sam side and ref side####
##as a full cycle) for each acquisition.
cycleindex	= (cyclecount+1)/2

##In the current data acquisition (Figure 2 in the paper), first 4.5 cycles#####
##are pre-PBL, last 5 cycles are the post-PBL.  In other words, user defined####
##on-peak cycle number is 14.
Urs_defined_cycle <- cycleindex - (pre_PBL+1 + post_PBL) - 1

##if the row number is not a continuous sequence, STOP.
if( length(which(cycleindex == 0)) != 0 )
  stop("Please make sure that Rows are continuous or contact the administrator.")

jstartp	= vector(mode = "numeric", length = nmbrow)

################################################################################
##PBL calculation create a baseline matrix to store all the baseline reading####
##as well as PREDICTED mass 44 intensity readings.
PBL_matrix <- matrix(nrow=(nmbrow), ncol=29)
colnames(PBL_matrix) <- c("row_ind",
                          "44_sa_start","44_sa_end",
                          "bl44_sa_start","bl44_sa_end","bl45_sa_start","bl45_sa_end","bl46_sa_start","bl46_sa_end",
                          "bl47_sa_start","bl47_sa_end","bl48_sa_start","bl48_sa_end","bl49_sa_start","bl49_sa_end",
                          "44_ref_start","44_ref_end",
                          "bl44_ref_start","bl44_ref_end","bl45_ref_start","bl45_ref_end","bl46_ref_start","bl46_ref_end",
                          "bl47_ref_start","bl47_ref_end","bl48_ref_start","bl48_ref_end","bl49_ref_start","bl49_ref_end"
)

##Create a scaler matrix that enables us to compute the PBL from 49 to other####
##masses.
bl49_Scaler_matrix <- matrix(nrow=( (pre_PBL-1)*nmbrow), ncol=25)
colnames(bl49_Scaler_matrix) <- c("row_ind",
                                  "44/49_sa_start","44/49_sa_end","45/49_sa_start","45/49_sa_end","46/49_sa_start","46/49_sa_end",
                                  "47/49_sa_start","47/49_sa_end","48/49_sa_start","48/49_sa_end","49/49_sa_start","49/49_sa_end",
                                  "44/49_ref_start","44/49_ref_end","45/49_ref_start","45/49_ref_end","46/49_ref_start","46/49_ref_end",
                                  "47/49_ref_start","47/49_ref_end","48/49_ref_start","48/49_ref_end","49/49_ref_start","49/49_ref_end"
)

bl_realtime_average <- matrix(nrow=( pre_PBL*nmbrow ), ncol=25)
colnames(bl_realtime_average) <- c("row_ind",
                                   "44bl_sa_start","44bl_sa_end","45bl_sa_start","45bl_sa_end","46bl_sa_start","46bl_sa_end",
                                   "47bl_sa_start","47bl_sa_end","48bl_sa_start","48bl_sa_end","49bl_sa_start","49bl_sa_end",
                                   "44bl_ref_start","44bl_ref_end","45bl_ref_start","45bl_ref_end","46bl_ref_start","46bl_ref_end",
                                   "47bl_ref_start","47bl_ref_end","48bl_ref_start","48bl_ref_end","49bl_ref_start","49bl_ref_end"
)

Sig49_44_realtime <- matrix(nrow=nmbrow, ncol=9)
colnames(Sig49_44_realtime) <- c("row_ind",
                                 "mass44_sa_start","mass_44_sa_end","sig49_sa_start","sig49_sa_end",
                                 "mass44_ref_start","mass_44_ref_end","sig49_ref_start","sig49_ref_end"
)

################################################################################
##We now loop to fill the matrices created from above.
for (blRowIndex in minrow:maxrow)
{
  bl <- blRowIndex - minrow + 1
  
  PBL_matrix[bl,1] <- blRowIndex
  
  selectRow  = rawtable[which(rawtable$Row == blRowIndex),]
  totline	= nrow(selectRow)
  totcycle	= (totline + 1)/2
  RowRearran	= matrix(nrow = totcycle, ncol=17)
  
  RowRearran[1, 1]	= selectRow[1, 1]
  RowRearran[1, 2]	= selectRow[1, 2]
  RowRearran[1, 3]	= selectRow[1, 3]
  RowRearran[1, 4:11]	= 0
  RowRearran[1, 12:17]= as.numeric( selectRow[1, 4:9] )
  for( n in 2:totcycle )
  {
    RowRearran[n, 1]	= selectRow[2*(n-1), 1]
    RowRearran[n, 2]	= selectRow[2*(n-1), 2]
    RowRearran[n, 3]	= selectRow[2*(n-1), 3]
    RowRearran[n, 4:11]	= as.numeric( selectRow[2*(n-1), 4:11] )
    RowRearran[n, 12:17]= as.numeric( selectRow[2*n-1, 4:9] )
  }
  colnames(RowRearran) = c("row", "identifier1", "method", "sa_intensity_44", "sa_intensity_45", "sa_intensity_46", "sa_intensity_47", "sa_intensity_48", "sa_intensity_49", "d_13C", "d_18O", "ref_intensity_44", "ref_intensity_45", "ref_intensity_46", "ref_intensity_47", "ref_intensity_48", "ref_intensity_49")
  
  meas_index_sa <- seq(from=(pre_PBL+1 + 1), by=1, to=(totcycle-post_PBL))
  meas_index_ref <- seq(from=(pre_PBL+1 + 1), by=1, to=(totcycle-post_PBL))
  
  ##This is the reverse of a local 2-pt averaging.
  PBL_matrix[bl,2] = 2*as.numeric(RowRearran[min(meas_index_sa),4]) - as.numeric(RowRearran[(min(meas_index_sa)+1),4])
  PBL_matrix[bl,3] = 2*as.numeric(RowRearran[max(meas_index_sa),4]) - as.numeric(RowRearran[(max(meas_index_sa)-1),4])
  
  PBL_matrix[bl,16] = 2*as.numeric(RowRearran[min(meas_index_ref),12]) - as.numeric(RowRearran[(min(meas_index_ref)+1),12])
  PBL_matrix[bl,17] = 2*as.numeric(RowRearran[max(meas_index_ref),12]) - as.numeric(RowRearran[(max(meas_index_ref)-1),12])
  
  ##The following can be looped through to simplify.
  PBL_matrix[bl,4] = as.numeric(RowRearran[(min(meas_index_sa)-1),4])
  PBL_matrix[bl,5] = as.numeric(RowRearran[(max(meas_index_sa)+1),4])
  PBL_matrix[bl,6] = as.numeric(RowRearran[(min(meas_index_sa)-1),5])
  PBL_matrix[bl,7] = as.numeric(RowRearran[(max(meas_index_sa)+1),5])
  PBL_matrix[bl,8] = as.numeric(RowRearran[(min(meas_index_sa)-1),6])
  PBL_matrix[bl,9] = as.numeric(RowRearran[(max(meas_index_sa)+1),6])
  PBL_matrix[bl,10] = as.numeric(RowRearran[(min(meas_index_sa)-1),7])
  PBL_matrix[bl,11] = as.numeric(RowRearran[(max(meas_index_sa)+1),7])
  PBL_matrix[bl,12] = as.numeric(RowRearran[(min(meas_index_sa)-1),8])
  PBL_matrix[bl,13] = as.numeric(RowRearran[(max(meas_index_sa)+1),8])
  PBL_matrix[bl,14] = as.numeric(RowRearran[(min(meas_index_sa)-1),9])
  PBL_matrix[bl,15] = as.numeric(RowRearran[(max(meas_index_sa)+1),9])
  
  PBL_matrix[bl,18] = as.numeric(RowRearran[(min(meas_index_ref)-1),12])
  PBL_matrix[bl,19] = as.numeric(RowRearran[(max(meas_index_ref)+1),12])
  PBL_matrix[bl,20] = as.numeric(RowRearran[(min(meas_index_ref)-1),13])
  PBL_matrix[bl,21] = as.numeric(RowRearran[(max(meas_index_ref)+1),13])
  PBL_matrix[bl,22] = as.numeric(RowRearran[(min(meas_index_ref)-1),14])
  PBL_matrix[bl,23] = as.numeric(RowRearran[(max(meas_index_ref)+1),14])
  PBL_matrix[bl,24] = as.numeric(RowRearran[(min(meas_index_ref)-1),15])
  PBL_matrix[bl,25] = as.numeric(RowRearran[(max(meas_index_ref)+1),15])
  PBL_matrix[bl,26] = as.numeric(RowRearran[(min(meas_index_ref)-1),16])
  PBL_matrix[bl,27] = as.numeric(RowRearran[(max(meas_index_ref)+1),16])
  PBL_matrix[bl,28] = as.numeric(RowRearran[(min(meas_index_ref)-1),17])
  PBL_matrix[bl,29] = as.numeric(RowRearran[(max(meas_index_ref)+1),17])
  ##Now we have a baseline data matrix.
  
  Scaler_start = (bl-1)*(pre_PBL-1) + 1
  Scaler_end = bl*(pre_PBL-1)
  
  bl49_Scaler_matrix[Scaler_start:Scaler_end,1] <- blRowIndex
  ## sample side
  bl49_Scaler_matrix[Scaler_start:Scaler_end,2] <- as.numeric(RowRearran[((2:pre_PBL)+1),4])/as.numeric(RowRearran[((2:pre_PBL)+1),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,3] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),4])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,4] <- as.numeric(RowRearran[((2:pre_PBL)+1),5])/as.numeric(RowRearran[((2:pre_PBL)+1),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,5] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),5])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,6] <- as.numeric(RowRearran[((2:pre_PBL)+1),6])/as.numeric(RowRearran[((2:pre_PBL)+1),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,7] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),6])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,8] <- as.numeric(RowRearran[((2:pre_PBL)+1),7])/as.numeric(RowRearran[((2:pre_PBL)+1),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,9] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),7])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,10] <- as.numeric(RowRearran[((2:pre_PBL)+1),8])/as.numeric(RowRearran[((2:pre_PBL)+1),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,11] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),8])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,12] <- as.numeric(RowRearran[((2:pre_PBL)+1),9])/as.numeric(RowRearran[((2:pre_PBL)+1),9])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,13] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),9])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),9])
  ## reference side
  bl49_Scaler_matrix[Scaler_start:Scaler_end,14] <- as.numeric(RowRearran[((2:pre_PBL)+1),12])/as.numeric(RowRearran[((2:pre_PBL)+1),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,15] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),12])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,16] <- as.numeric(RowRearran[((2:pre_PBL)+1),13])/as.numeric(RowRearran[((2:pre_PBL)+1),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,17] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),13])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,18] <- as.numeric(RowRearran[((2:pre_PBL)+1),14])/as.numeric(RowRearran[((2:pre_PBL)+1),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,19] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),14])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,20] <- as.numeric(RowRearran[((2:pre_PBL)+1),15])/as.numeric(RowRearran[((2:pre_PBL)+1),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,21] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),15])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,22] <- as.numeric(RowRearran[((2:pre_PBL)+1),16])/as.numeric(RowRearran[((2:pre_PBL)+1),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,23] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),16])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,24] <- as.numeric(RowRearran[((2:pre_PBL)+1),17])/as.numeric(RowRearran[((2:pre_PBL)+1),17])
  bl49_Scaler_matrix[Scaler_start:Scaler_end,25] <- as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),17])/as.numeric(RowRearran[((totcycle-post_PBL+3):totcycle),17])
  
  
  blIndex_start <- (bl-1)*(pre_PBL) + 1
  blIndex_end <- bl*(pre_PBL)
  ##Now, average baseline.
  bl_realtime_average[blIndex_start:blIndex_end,1] <- blRowIndex
  bl_realtime_average[blIndex_start:blIndex_end,2] <- as.numeric(RowRearran[((1:pre_PBL)+1),4])
  bl_realtime_average[blIndex_start:blIndex_end,3] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),4])
  bl_realtime_average[blIndex_start:blIndex_end,4] <- as.numeric(RowRearran[((1:pre_PBL)+1),5])
  bl_realtime_average[blIndex_start:blIndex_end,5] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),5])
  bl_realtime_average[blIndex_start:blIndex_end,6] <- as.numeric(RowRearran[((1:pre_PBL)+1),6])
  bl_realtime_average[blIndex_start:blIndex_end,7] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),6])
  bl_realtime_average[blIndex_start:blIndex_end,8] <- as.numeric(RowRearran[((1:pre_PBL)+1),7])
  bl_realtime_average[blIndex_start:blIndex_end,9] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),7])
  bl_realtime_average[blIndex_start:blIndex_end,10] <- as.numeric(RowRearran[((1:pre_PBL)+1),8])
  bl_realtime_average[blIndex_start:blIndex_end,11] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),8])
  bl_realtime_average[blIndex_start:blIndex_end,12] <- as.numeric(RowRearran[((1:pre_PBL)+1),9])
  bl_realtime_average[blIndex_start:blIndex_end,13] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),9])
  bl_realtime_average[blIndex_start:blIndex_end,14] <- as.numeric(RowRearran[((1:pre_PBL)+1),12])
  bl_realtime_average[blIndex_start:blIndex_end,15] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),12])
  bl_realtime_average[blIndex_start:blIndex_end,16] <- as.numeric(RowRearran[((1:pre_PBL)+1),13])
  bl_realtime_average[blIndex_start:blIndex_end,17] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),13])
  bl_realtime_average[blIndex_start:blIndex_end,18] <- as.numeric(RowRearran[((1:pre_PBL)+1),14])
  bl_realtime_average[blIndex_start:blIndex_end,19] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),14])
  bl_realtime_average[blIndex_start:blIndex_end,20] <- as.numeric(RowRearran[((1:pre_PBL)+1),15])
  bl_realtime_average[blIndex_start:blIndex_end,21] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),15])
  bl_realtime_average[blIndex_start:blIndex_end,22] <- as.numeric(RowRearran[((1:pre_PBL)+1),16])
  bl_realtime_average[blIndex_start:blIndex_end,23] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),16])
  bl_realtime_average[blIndex_start:blIndex_end,24] <- as.numeric(RowRearran[((1:pre_PBL)+1),17])
  bl_realtime_average[blIndex_start:blIndex_end,25] <- as.numeric(RowRearran[((totcycle-post_PBL+2):(totcycle)),17])
  
  ##### mass49 signal for realtime regression.
  Sig49_44_realtime[bl, 1] <- blRowIndex
  Sig49_44_realtime[bl, 2] <- as.numeric(RowRearran[(min(meas_index_sa)),4])
  Sig49_44_realtime[bl, 3] <- as.numeric(RowRearran[(max(meas_index_sa)),4])
  Sig49_44_realtime[bl, 4] <- ( as.numeric(RowRearran[(min(meas_index_sa)),9]) - (2*as.numeric(RowRearran[(min(meas_index_sa)-1),9])-as.numeric(RowRearran[(min(meas_index_sa)-2),9])) )
  Sig49_44_realtime[bl, 5] <- ( as.numeric(RowRearran[(max(meas_index_sa)),9]) - (3*as.numeric(RowRearran[(max(meas_index_sa)+2),9])-2*as.numeric(RowRearran[(max(meas_index_sa)+3),9])) )
  
  Sig49_44_realtime[bl, 6] <- as.numeric(RowRearran[(min(meas_index_sa)),12])
  Sig49_44_realtime[bl, 7] <- as.numeric(RowRearran[(max(meas_index_sa)),12])
  Sig49_44_realtime[bl, 8] <- ( as.numeric(RowRearran[(min(meas_index_sa)),17]) - (2*as.numeric(RowRearran[(min(meas_index_sa)-1),17])-as.numeric(RowRearran[(min(meas_index_sa)-2),17])) )
  Sig49_44_realtime[bl, 9] <- ( as.numeric(RowRearran[(max(meas_index_sa)),17]) - (3*as.numeric(RowRearran[(max(meas_index_sa)+2),17])-2*as.numeric(RowRearran[(max(meas_index_sa)+3),17])) )
}


################################################################################
################################################################################
##DEFINE output matrix!

calc_output	= matrix(ncol = 30, nrow = sum(Urs_defined_cycle))
colnames(calc_output)  = c("row_index_number", "method_info", "number_cycle", "sa_d13_VPDB_mean", "sa_d18_VSMOW_mean", "d45", "d46", "d47", "d48", "d49",
                           "D45", "D46", "D47", "D48", "D49", "FullD47", "FullD48", "FullD49", "FullD47_mean", "FullD47_sd", "FullD48_mean", "FullD48_sd", "FullD49_mean", "FullD49_sd",
                           "sa_d13_mean", "sa_d13_sd", "sa_d18_mean", "sa_d18_sd", "cum_ave_FullD47", "cum_sd_Full47")

for (i in minrow:maxrow)
{
  
  subset	= rawtable[which(rawtable$Row == i),]
  nmbline	= nrow(subset)
  nmbcycle	= (nmbline + 1)/2
  subtable	= matrix(nrow = nmbcycle, ncol=19)
  ##NOTE, we had inserted one more sequence line in the acquisition.
  MeasurementLoop_cycle_realtime <- nmbcycle - pre_PBL-1 - post_PBL
  MeasurementLoop_cycle_usr <- MeasurementLoop_cycle_realtime - 1  
  MeasurementLoop_inDex <- seq(from=(pre_PBL+2), by=1, length=(MeasurementLoop_cycle_usr+1) )
  
  SwitchIndex <- which(PBL_matrix[,1] == i)
  
  subtable[1, 1]	= subset[1, 1]
  subtable[1, 2]	= subset[1, 2]
  subtable[1, 3]	= subset[1, 3]
  subtable[1, 4:11]	= 0
  subtable[1, 12:19]= as.numeric( subset[1, 4:11] )
  
  for( n in 2:nmbcycle)
  {
    subtable[n, 1]	= subset[2*(n-1), 1]
    subtable[n, 2]	= subset[2*(n-1), 2]
    subtable[n, 3]	= subset[2*(n-1), 3]
    subtable[n, 4:11]	= as.numeric( subset[2*(n-1), 4:11] )
    subtable[n, 12:19]= as.numeric( subset[2*n-1, 4:11] )
  }
  colnames(subtable) = c("row", "identifier1", "method", "sa_intensity_44", "sa_intensity_45", "sa_intensity_46", "sa_intensity_47", "sa_intensity_48", "sa_intensity_49", "d_13C", "d_18O", "ref_intensity_44", "ref_intensity_45", "ref_intensity_46", "ref_intensity_47", "ref_intensity_48", "ref_intensity_49", "d13_ref","d18_ref")
  
  
  ##At sample side, we ignore the first cycle.
  sa_44_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex[-1], 4])
  sa_45_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex[-1], 5])
  sa_46_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex[-1], 6])
  sa_47_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex[-1], 7])
  sa_48_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex[-1], 8])
  sa_49_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex[-1], 9])
  
  ##At ref side, we have one more line.
  ref_44_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex, 12])
  ref_45_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex, 13])
  ref_46_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex, 14])
  ref_47_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex, 15])
  ref_48_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex, 16])
  ref_49_raw_raw	= as.numeric(subtable[MeasurementLoop_inDex, 17])
  
  ################################################################################
  ################################################################################
  blScaler_interpoInd <- c( ((1:pre_PBL)+1), ((nmbcycle-post_PBL+2):nmbcycle) )[-c(1,5)]  
  bl_Scaler_subtable <- bl49_Scaler_matrix[which(bl49_Scaler_matrix[,1] == i),]
  
  blScaler44_2cluster_sa = c(bl_Scaler_subtable[,2], bl_Scaler_subtable[,3])
  blScaler45_2cluster_sa = c(bl_Scaler_subtable[,4], bl_Scaler_subtable[,5])
  blScaler46_2cluster_sa = c(bl_Scaler_subtable[,6], bl_Scaler_subtable[,7])
  blScaler47_2cluster_sa = c(bl_Scaler_subtable[,8], bl_Scaler_subtable[,9])
  blScaler48_2cluster_sa = c(bl_Scaler_subtable[,10], bl_Scaler_subtable[,11])
  blScaler49_2cluster_sa = c(bl_Scaler_subtable[,12], bl_Scaler_subtable[,13])
  
  blScaler44_2cluster_ref = c(bl_Scaler_subtable[,14], bl_Scaler_subtable[,15])
  blScaler45_2cluster_ref = c(bl_Scaler_subtable[,16], bl_Scaler_subtable[,17])
  blScaler46_2cluster_ref = c(bl_Scaler_subtable[,18], bl_Scaler_subtable[,19])
  blScaler47_2cluster_ref = c(bl_Scaler_subtable[,20], bl_Scaler_subtable[,21])
  blScaler48_2cluster_ref = c(bl_Scaler_subtable[,22], bl_Scaler_subtable[,23])
  blScaler49_2cluster_ref = c(bl_Scaler_subtable[,24], bl_Scaler_subtable[,25])
  
  blScaler44_interpo_sa <- approx(blScaler_interpoInd, blScaler44_2cluster_sa, method=c("linear"), n=19)
  blScaler45_interpo_sa <- approx(blScaler_interpoInd, blScaler45_2cluster_sa, method=c("linear"), n=19)
  blScaler46_interpo_sa <- approx(blScaler_interpoInd, blScaler46_2cluster_sa, method=c("linear"), n=19)
  blScaler47_interpo_sa <- approx(blScaler_interpoInd, blScaler47_2cluster_sa, method=c("linear"), n=19)
  blScaler48_interpo_sa <- approx(blScaler_interpoInd, blScaler48_2cluster_sa, method=c("linear"), n=19)
  blScaler49_interpo_sa <- approx(blScaler_interpoInd, blScaler49_2cluster_sa, method=c("linear"), n=19)
  
  blScaler44_interpo_ref <- approx(blScaler_interpoInd, blScaler44_2cluster_ref, method=c("linear"), n=19)
  blScaler45_interpo_ref <- approx(blScaler_interpoInd, blScaler45_2cluster_ref, method=c("linear"), n=19)
  blScaler46_interpo_ref <- approx(blScaler_interpoInd, blScaler46_2cluster_ref, method=c("linear"), n=19)
  blScaler47_interpo_ref <- approx(blScaler_interpoInd, blScaler47_2cluster_ref, method=c("linear"), n=19)
  blScaler48_interpo_ref <- approx(blScaler_interpoInd, blScaler48_2cluster_ref, method=c("linear"), n=19)
  blScaler49_interpo_ref <- approx(blScaler_interpoInd, blScaler49_2cluster_ref, method=c("linear"), n=19)
  
  blScaler44_ref <- blScaler44_interpo_ref$y[4:14]
  blScaler45_ref <- blScaler45_interpo_ref$y[4:14]
  blScaler46_ref <- blScaler46_interpo_ref$y[4:14]
  blScaler47_ref <- blScaler47_interpo_ref$y[4:14]
  blScaler48_ref <- blScaler48_interpo_ref$y[4:14]
  blScaler49_ref <- blScaler49_interpo_ref$y[4:14]
  
  blScaler44_sa <- blScaler44_interpo_sa$y[5:14]
  blScaler45_sa <- blScaler45_interpo_sa$y[5:14]
  blScaler46_sa <- blScaler46_interpo_sa$y[5:14]
  blScaler47_sa <- blScaler47_interpo_sa$y[5:14]
  blScaler48_sa <- blScaler48_interpo_sa$y[5:14]
  blScaler49_sa <- blScaler49_interpo_sa$y[5:14]
  
  ##Calculate the real time 49 signal response.
  Sig49_44_realtime_reg <- lm(c(Sig49_44_realtime[,4:5], Sig49_44_realtime[,8:9])~c(Sig49_44_realtime[,2:3], Sig49_44_realtime[,6:7]))
  Sig49_44_coef <- Sig49_44_realtime_reg$coefficients[2]
  Sig49_44_intp <- Sig49_44_realtime_reg$coefficients[1]
  
  Sig49_44_realtime_reg_sa <- lm(c(Sig49_44_realtime[,4:5])~c(Sig49_44_realtime[,2:3]))
  Sig49_44_coef_sa <- Sig49_44_realtime_reg_sa$coefficients[2]
  Sig49_44_intp_sa <- Sig49_44_realtime_reg_sa$coefficients[1]
  
  Sig49_44_realtime_reg_ref <- lm(c(Sig49_44_realtime[,8:9])~c(Sig49_44_realtime[,6:7]))
  Sig49_44_coef_ref <- Sig49_44_realtime_reg_ref$coefficients[2]
  Sig49_44_intp_ref <- Sig49_44_realtime_reg_ref$coefficients[1]
  
  #####
  subTableIndex <- which(bl_realtime_average[,1] == i)
  bl_44_endPoints_sa <- c( (bl_realtime_average[subTableIndex,2]), (bl_realtime_average[subTableIndex,3]) )[-c(1,2,5,6)]
  bl_45_endPoints_sa <- c( (bl_realtime_average[subTableIndex,4]), (bl_realtime_average[subTableIndex,5]) )[-c(1,2,5,6)]
  bl_46_endPoints_sa <- c( (bl_realtime_average[subTableIndex,6]), (bl_realtime_average[subTableIndex,7]) )[-c(1,2,5,6)]
  bl_47_endPoints_sa <- c( (bl_realtime_average[subTableIndex,8]), (bl_realtime_average[subTableIndex,9]) )[-c(1,2,5,6)]
  bl_48_endPoints_sa <- c( (bl_realtime_average[subTableIndex,10]), (bl_realtime_average[subTableIndex,11]) )[-c(1,2,5,6)]
  bl_49_endPoints_sa <- c( (bl_realtime_average[subTableIndex,12]), (bl_realtime_average[subTableIndex,13]) )[-c(1,2,5,6)]
  bl_44_endPoints_ref <- c( (bl_realtime_average[subTableIndex,14]), (bl_realtime_average[subTableIndex,15]) )[-c(1,2,5,6)]
  bl_45_endPoints_ref <- c( (bl_realtime_average[subTableIndex,16]), (bl_realtime_average[subTableIndex,17]) )[-c(1,2,5,6)]
  bl_46_endPoints_ref <- c( (bl_realtime_average[subTableIndex,18]), (bl_realtime_average[subTableIndex,19]) )[-c(1,2,5,6)]
  bl_47_endPoints_ref <- c( (bl_realtime_average[subTableIndex,20]), (bl_realtime_average[subTableIndex,21]) )[-c(1,2,5,6)]
  bl_48_endPoints_ref <- c( (bl_realtime_average[subTableIndex,22]), (bl_realtime_average[subTableIndex,23]) )[-c(1,2,5,6)]
  bl_49_endPoints_ref <- c( (bl_realtime_average[subTableIndex,24]), (bl_realtime_average[subTableIndex,25]) )[-c(1,2,5,6)]
  
  linear_interpo_temp <- function(data, bl_interpo_Index)
  {
    temp_mdl <- lm(as.numeric(data)~as.numeric(bl_interpo_Index))
    pred_index <- c(min(bl_interpo_Index):max(bl_interpo_Index))
    pred_value <- temp_mdl$coefficients[1] + temp_mdl$coefficients[2] * pred_index
    x <- pred_index
    y <- pred_value
    return(list(x,y))
  }
  
  bl_interpo_Index <- c( ((1:pre_PBL)+1), ((nmbcycle-post_PBL+1):nmbcycle) )[-c(1,2,5,6,7)]
  
  bl44_sa_linear_spline <- linear_interpo_temp(bl_44_endPoints_sa, bl_interpo_Index)
  bl45_sa_linear_spline <- linear_interpo_temp(bl_45_endPoints_sa, bl_interpo_Index)
  bl46_sa_linear_spline <- linear_interpo_temp(bl_46_endPoints_sa, bl_interpo_Index)
  bl47_sa_linear_spline <- linear_interpo_temp(bl_47_endPoints_sa, bl_interpo_Index)
  bl48_sa_linear_spline <- linear_interpo_temp(bl_48_endPoints_sa, bl_interpo_Index)
  bl49_sa_linear_spline <- linear_interpo_temp(bl_49_endPoints_sa, bl_interpo_Index)
  
  bl44_ref_linear_spline <- linear_interpo_temp(bl_44_endPoints_ref, bl_interpo_Index)
  bl45_ref_linear_spline <- linear_interpo_temp(bl_45_endPoints_ref, bl_interpo_Index)
  bl46_ref_linear_spline <- linear_interpo_temp(bl_46_endPoints_ref, bl_interpo_Index)
  bl47_ref_linear_spline <- linear_interpo_temp(bl_47_endPoints_ref, bl_interpo_Index)
  bl48_ref_linear_spline <- linear_interpo_temp(bl_48_endPoints_ref, bl_interpo_Index)
  bl49_ref_linear_spline <- linear_interpo_temp(bl_49_endPoints_ref, bl_interpo_Index)
  
  ###########options for PBL correction HERE, with ifelse.
  if(PBL_49_scaling == TRUE)
  {
    ## subtract off the signal component first!!
    sa_49_raw_bl = sa_49_raw_raw - (Sig49_44_coef_sa * sa_44_raw_raw + Sig49_44_intp_sa)
    ref_49_raw_bl = ref_49_raw_raw - (Sig49_44_coef_ref * ref_44_raw_raw + Sig49_44_intp_ref)
    
    ref_44_raw <- ref_44_raw_raw - ref_49_raw_bl*blScaler44_ref
    ref_45_raw <- ref_45_raw_raw - ref_49_raw_bl*blScaler45_ref
    ref_46_raw <- ref_46_raw_raw - ref_49_raw_bl*blScaler46_ref
    ref_47_raw <- ref_47_raw_raw - ref_49_raw_bl*blScaler47_ref
    ref_48_raw <- ref_48_raw_raw - ref_49_raw_bl*blScaler48_ref
    ref_49_raw <- ref_49_raw_raw - ref_49_raw_bl*blScaler49_ref
    
    sa_44_raw <- sa_44_raw_raw - sa_49_raw_bl*blScaler44_sa
    sa_45_raw <- sa_45_raw_raw - sa_49_raw_bl*blScaler45_sa
    sa_46_raw <- sa_46_raw_raw - sa_49_raw_bl*blScaler46_sa
    sa_47_raw <- sa_47_raw_raw - sa_49_raw_bl*blScaler47_sa
    sa_48_raw <- sa_48_raw_raw - sa_49_raw_bl*blScaler48_sa
    sa_49_raw <- sa_49_raw_raw - sa_49_raw_bl*blScaler49_sa
  }else
    if(PBL_47_interpo == TRUE)
    {
      bl_pick_min <- which(bl44_sa_linear_spline[[1]] == min(MeasurementLoop_inDex))
      bl_pick_max <- which(bl44_sa_linear_spline[[1]] == max(MeasurementLoop_inDex))
      
      sa_44_raw <- sa_44_raw_raw - bl44_sa_linear_spline[[2]][(bl_pick_min+1):bl_pick_max]
      sa_45_raw <- sa_45_raw_raw - bl45_sa_linear_spline[[2]][(bl_pick_min+1):bl_pick_max]
      sa_46_raw <- sa_46_raw_raw - bl46_sa_linear_spline[[2]][(bl_pick_min+1):bl_pick_max]
      sa_47_raw <- sa_47_raw_raw - bl47_sa_linear_spline[[2]][(bl_pick_min+1):bl_pick_max]
      sa_48_raw <- sa_48_raw_raw - bl48_sa_linear_spline[[2]][(bl_pick_min+1):bl_pick_max]
      sa_49_raw <- sa_49_raw_raw - bl49_sa_linear_spline[[2]][(bl_pick_min+1):bl_pick_max]
      
      ref_44_raw <- ref_44_raw_raw - bl44_ref_linear_spline[[2]][bl_pick_min:bl_pick_max]
      ref_45_raw <- ref_45_raw_raw - bl45_ref_linear_spline[[2]][bl_pick_min:bl_pick_max]
      ref_46_raw <- ref_46_raw_raw - bl46_ref_linear_spline[[2]][bl_pick_min:bl_pick_max]
      ref_47_raw <- ref_47_raw_raw - bl47_ref_linear_spline[[2]][bl_pick_min:bl_pick_max]
      ref_48_raw <- ref_48_raw_raw - bl48_ref_linear_spline[[2]][bl_pick_min:bl_pick_max]
      ref_49_raw <- ref_49_raw_raw - bl49_ref_linear_spline[[2]][bl_pick_min:bl_pick_max]
    }else{
      sa_44_raw <- sa_44_raw_raw
      sa_45_raw <- sa_45_raw_raw
      sa_46_raw <- sa_46_raw_raw
      sa_47_raw <- sa_47_raw_raw
      sa_48_raw <- sa_48_raw_raw
      sa_49_raw <- sa_49_raw_raw
      ref_44_raw <- ref_44_raw_raw
      ref_45_raw <- ref_45_raw_raw
      ref_46_raw <- ref_46_raw_raw
      ref_47_raw <- ref_47_raw_raw
      ref_48_raw <- ref_48_raw_raw
      ref_49_raw <- ref_49_raw_raw
    }
  
  ncycles <- MeasurementLoop_cycle_usr + 1
  ## or, ncycles <- length(MeasurementLoop_inDex)
  diff_44	= ref_44_raw[2:ncycles] - sa_44_raw
  diff_47_44	= (ref_47_raw/ref_44_raw)[2:ncycles] - (sa_47_raw/sa_44_raw)
  
  # list calcated sa_d13C_VPDB and sa_d18O_VSMOW from isodat 3.0
  sa_d13C	= as.numeric(subtable[MeasurementLoop_inDex[-1], 10])
  sa_d18O	= as.numeric(subtable[MeasurementLoop_inDex[-1], 11])
  
  RefGasd13 = as.numeric(subtable[1,18])
  RefGasd18 = as.numeric(subtable[1,19])
  
  n		= length(sa_d13C)
  f_ref_49	= ref_47_raw/1000
  f_sa_49 <- (f_ref_49[2:ncycles] + f_ref_49[1:(ncycles-1)])/2+sa_49_raw-(ref_49_raw[2:ncycles]+ref_49_raw[1:(ncycles-1)])/2
  
  ##Define and calculate R_i=[mass_i]/[mass_44], where i=[45,46,47,48,49] of sample & standard  
  
  Rsa45 	= sa_45_raw/sa_44_raw
  Rsa46 	= sa_46_raw/sa_44_raw
  Rsa47 	= sa_47_raw/sa_44_raw
  Rsa48		= sa_48_raw/sa_44_raw
  Rsa49		= f_sa_49/sa_44_raw
  
  Rref45	= ref_45_raw/ref_44_raw
  Rref46	= ref_46_raw/ref_44_raw
  Rref47	= ref_47_raw/ref_44_raw
  Rref48	= ref_48_raw/ref_44_raw
  Rref49	= f_ref_49/ref_44_raw
  
  ##d_i=(Rsa_i/Rst_i)-1)*1000, where i = [45,46,47,48,49] of sample & standard. Rst_i is R_i value for hypothetical CO2 with d13C=0, d18O=0, D47=0
  
  d45		= (Rsa45/( (Rref45[2:ncycles] + Rref45[1:(ncycles-1)])/2)-1)*1000
  d46		= (Rsa46/( (Rref46[2:ncycles] + Rref46[1:(ncycles-1)])/2)-1)*1000
  d47		= (Rsa47/( (Rref47[2:ncycles] + Rref47[1:(ncycles-1)])/2)-1)*1000
  d48		= (Rsa48/( (Rref48[2:ncycles] + Rref48[1:(ncycles-1)])/2)-1)*1000
  d49		= (Rsa49/( (Rref49[2:ncycles] + Rref49[1:(ncycles-1)])/2)-1)*1000
  
  ################################################################################
  ##Scramble the referene gas	side################################################
  ################################################################################
  
  # REF GAS: calculate R13, R17, R18 
  # minimum value of (R13+2R18-R45)^2+(2R18+2R13R17+(R17)^2-R46)^2
  
  RefGasR13	= (1+RefGasd13/1000)*VPDBcalcite_R13;
  RefGasR18	= (1+RefGasd18/1000)*VSMOW_R18;
  RefGasR17	= VSMOW_R17*(RefGasR18/VSMOW_R18)^lambda;
  
  # REF GAS: calculate abundance of isotopes, assuming stochastic
  # distribution of isotopoes among isotopologues
  
  RefGasC12 = 1/(1+RefGasR13)			#[12C]=1/(1+R13)
  RefGasC13 = RefGasR13 * RefGasC12		#[13C]=R13/(1+R13)
  RefGasC16 = 1/(1+RefGasR17+RefGasR18)	#[16O]=1/(1+R17+R18)
  RefGasC17 = RefGasR17 * RefGasC16		#[17O]=R17/(1+R17+R18)
  RefGasC18 = RefGasR18 * RefGasC16		#[18O]=R18/(1+R17+R18)
  
  # REF GAS: calculate randomly distributed abundance of CO2 
  # [44]* = [12][16][16]
  # [45]* = [13][16][16]+2[12][16][17]
  # [46]* = 2[12][16][18]+[12][17][17]+2[13][16][17]
  # [47]* = 2[13][16][18]+[13][17][17]+2[12][17][18]
  
  RefGasC121616	= RefGasC12 * RefGasC16^2
  RefGasC131616	= RefGasC13 * RefGasC16^2
  RefGasC121716	= RefGasC12 * RefGasC17 * RefGasC16 * 2
  RefGasC131716	= RefGasC13 * RefGasC17 * RefGasC16 * 2
  RefGasC121816	= RefGasC12 * RefGasC18 * RefGasC16 * 2
  RefGasC121717	= RefGasC12 * RefGasC17^2
  RefGasC131816	= RefGasC13 * RefGasC18 * RefGasC16 * 2
  RefGasC131717	= RefGasC13 * RefGasC17^2
  RefGasC121817	= RefGasC12 * RefGasC18 * RefGasC17 * 2
  RefGasC131817	= RefGasC13 * RefGasC18 * RefGasC17 * 2
  RefGasC121818	= RefGasC12 * RefGasC18^2
  RefGasC131818	= RefGasC13 * RefGasC18^2
  
  RefGasC44		= RefGasC121616
  RefGasC45		= RefGasC131616 + RefGasC121716
  RefGasC46		= RefGasC131716 + RefGasC121816 + RefGasC121717
  RefGasC47		= RefGasC131816 + RefGasC121817 + RefGasC131717
  RefGasC48		= RefGasC131817 + RefGasC121818
  RefGasC49		= RefGasC131818
  
  # REF GAS: calculate ratios by dividing abundance by [44]*
  
  RefGasR45		= RefGasC45/RefGasC44
  RefGasR46		= RefGasC46/RefGasC44
  RefGasR47		= RefGasC47/RefGasC44
  RefGasR48		= RefGasC48/RefGasC44
  RefGasR49		= RefGasC49/RefGasC44
  
  # SAMPLE: calculate R45 through R49, inclusive
  
  R45			= (1+d45/1000)*RefGasR45;
  R46			= (1+d46/1000)*RefGasR46;
  R47			= (1+d47/1000)*RefGasR47;
  R48			= (1+d48/1000)*RefGasR48;
  R49			= (1+d49/1000)*RefGasR49;     
  
  ################################################################################
  ##Scramble the sample gas	side##################################################
  ################################################################################
  
  # SAMPLE: calculate R13, R17, R18 
  
  ScrSampleR13	= (1 + sa_d13C/1000) * VPDBcalcite_R13
  ScrSampleR18	= (1 + sa_d18O/1000) * VSMOW_R18
  ScrSampleR17	= VSMOW_R17*(ScrSampleR18/VSMOW_R18)^lambda
  
  # SAMPLE: calculate abundance of isotopes, assuming stochastic
  # distribution of isotopoes among isotopologues
  
  ScrSampleC12	= 1/(1 + ScrSampleR13)				#[12C]=1/(1+R13)
  ScrSampleC13	= ScrSampleR13 * ScrSampleC12			#[13C]=R13/(1+R13)
  ScrSampleC16	= 1/(1 + ScrSampleR17 + ScrSampleR18)	#[16O]=1/(1+R17+R18)
  ScrSampleC17	= ScrSampleR17 * ScrSampleC16			#[17O]=R17/(1+R17+R18)
  ScrSampleC18	= ScrSampleR18 * ScrSampleC16			#[18O]=R18/(1+R17+R18)
  
  # SAMPLE: calculate randomly distributed abundance of CO2 
  # [44] = [12][16][16]
  # [45] = [13][16][16]+2[12][16][17]
  # [46] = 2[12][16][18]+[12][17][17]+2[13][16][17]
  # [47] = 2[13][16][18]+[13][17][17]+2[12][17][18]
  
  ScrSampleC121616	= ScrSampleC12 * ScrSampleC16^2
  ScrSampleC131616	= ScrSampleC13 * ScrSampleC16^2
  ScrSampleC121716	= ScrSampleC12 * ScrSampleC17 * ScrSampleC16 * 2
  ScrSampleC131716	= ScrSampleC13 * ScrSampleC17 * ScrSampleC16 * 2
  ScrSampleC121816	= ScrSampleC12 * ScrSampleC18 * ScrSampleC16 * 2
  ScrSampleC121717	= ScrSampleC12 * ScrSampleC17^2
  ScrSampleC131816	= ScrSampleC13 * ScrSampleC18 * ScrSampleC16 * 2
  ScrSampleC131717	= ScrSampleC13 * ScrSampleC17^2
  ScrSampleC121817	= ScrSampleC12 * ScrSampleC18 * ScrSampleC17 * 2
  ScrSampleC131817	= ScrSampleC13 * ScrSampleC18 * ScrSampleC17 * 2
  ScrSampleC121818	= ScrSampleC12 * ScrSampleC18^2
  ScrSampleC131818	= ScrSampleC13 * ScrSampleC18^2
  
  ScrSampleC44	= ScrSampleC121616
  ScrSampleC45	= ScrSampleC131616 + ScrSampleC121716
  ScrSampleC46	= ScrSampleC131716 + ScrSampleC121816 + ScrSampleC121717
  ScrSampleC47	= ScrSampleC131816 + ScrSampleC121817 + ScrSampleC131717
  ScrSampleC48	= ScrSampleC131817 + ScrSampleC121818
  ScrSampleC49	= ScrSampleC131818
  
  # SAMPLE: calculate ratios by dividing abundance by [44]
  
  ScrSampleR45 = ScrSampleC45/ScrSampleC44
  ScrSampleR46 = ScrSampleC46/ScrSampleC44
  ScrSampleR47 = ScrSampleC47/ScrSampleC44
  ScrSampleR48 = ScrSampleC48/ScrSampleC44
  ScrSampleR49 = ScrSampleC49/ScrSampleC44
  
  
  ################################################################################
  #####		Compare measured sample R values to calculated scrambled R values ######
  ################################################################################
  
  # D_i values for the sample
  
  D45			=(R45/ScrSampleR45 -1)*1000
  D46			=(R46/ScrSampleR46 -1)*1000
  D47			=(R47/ScrSampleR47 -1)*1000
  D48			=(R48/ScrSampleR48 -1)*1000
  D49			=(R49/ScrSampleR49 -1)*1000
  
  
  # mass 47, 48, 49 anomalies for the sample, showing isotopologue
  # abunndance with respect to the stochastic distribution.
  
  FullD47 = D47 - D46 - D45	
  FullD48 = D48 - 2*D46			
  FullD49 = D49 - 2*D46 - D45	
  
  ################################################################################
  #####		Please refer to the test data and compare results with the following####
  ################################################################################
  # The results calculated with PBL correction:
  # D45 = 0.05809299
  # D46 = -0.002678686
  # D47 = -0.7574661
  # D48 = -0.2864781
  # D49 = 92.49767
  # FullD47 = -0.8128804
  # FullD48 = -0.2811207
  # FullD49 = 92.44493
  
  # The results calculated w/o PBL correction:
  # D45 = 0.05864693
  # D46 = -0.001676277
  # D47 = -0.9770799
  # D48 = -0.7226561
  # D49 = -797.6971
  # FullD47 = -1.034051
  # FullD48 = -0.7193035
  # FullD49 = -797.7524
  
  j = (i - minrow + 1)
  
  if (j == 1)
  {
    jstart = 1
    jend = jstart +Urs_defined_cycle[j] - 1
  }else{
    jstart = sum(Urs_defined_cycle[1:(j-1)]) + 1
    jend = jstart +Urs_defined_cycle[j] - 1
  }
  
  jstartp[j] = jstart
  
  calc_output[jstart:jend, 1]		= subtable[1,1]
  calc_output[jstart:jend, 2]		= subtable[1,3]
  calc_output[jstart:jend, 3]		= ncycles
  calc_output[jstart:jend, 4]		= sa_d13C
  calc_output[jstart:jend, 5]		= sa_d18O
  calc_output[jstart:jend, 6]		= d45
  calc_output[jstart:jend, 7]		= d46
  calc_output[jstart:jend, 8]		= d47
  calc_output[jstart:jend, 9]		= d48
  calc_output[jstart:jend, 10]		= d49
  calc_output[jstart:jend, 11]		= D45
  calc_output[jstart:jend, 12]		= D46
  calc_output[jstart:jend, 13]		= D47
  calc_output[jstart:jend, 14]		= D48
  calc_output[jstart:jend, 15]		= D49
  calc_output[jstart:jend, 16]		= FullD47
  calc_output[jstart:jend, 17]		= FullD48
  calc_output[jstart:jend, 18]		= FullD49
  calc_output[jstart:jend, 19]		= mean(FullD47)
  calc_output[jstart:jend, 20]		= sd(FullD47)
  calc_output[jstart:jend, 21]		= mean(FullD48)
  calc_output[jstart:jend, 22]		= sd(FullD49)
  calc_output[jstart:jend, 23]		= mean(FullD49)
  calc_output[jstart:jend, 24]		= sd(FullD49)
  calc_output[jstart:jend, 25]		= mean(sa_d13C)
  calc_output[jstart:jend, 26]		= sd(sa_d13C)
  calc_output[jstart:jend, 27]		= mean(sa_d18O)
  calc_output[jstart:jend, 28]		= sd(sa_d18O)
  calc_output[jstart:jend, 29]		= 0
  calc_output[jstart:jend, 30]		= 0
  
}

################################################################################
##Summarize all the data and export as a seperate file##########################
setwd(save_dir)
write.csv(calc_output, paste("data_summary_", x), quote=F)

################################################################################
plot(calc_output[,16])
abline(h=mean(as.numeric(calc_output[,16])), col="grey", lwd=1.5)
mean(as.numeric(calc_output[,16]))
################################################################################
#####                 END             	########################################
################################################################################