##the transform can be inverted (use for sensor fusion?)
#  spec.inv <- Re(mvfft(spec.fft, inverse = T))


########################################################
#FAST FOURIER TRANSFORM
########################################################
SPEC.fwd <- function(
  tstamp=POST$UT,
  data,
  relmot=POST$veloTAmagn,
  height=POST$lengMCzLS,
  fre=freq,
  demean=F,
  detrend=T,
  taper=0.05)
{
#in:
  #tstamp: continuous timestamp in any float format; only for gapfilling
  #data: continuous matrix of variables to be transformed with same length as tstamp; NAs are interpolated
  #relmot: relative motion between observation platform and atmosphere with same length as tstamp; NAs are interpolated. E.g. |wind vector| for tower observation [m s-1]
  #height: measurement height, either single number or vector of same length as tstamp
  #fre: observation frequency: single integer [Hz]
  #demean: demean the data?
  #detrend: detrend the data?
  #taper: tapter the data for a fraction (<=0.5) at both ends, or F
#out:
  #transformation of variables in frequency space and calculation of various independent variables
  #sorted from lowest to highest frequency

#combine variables
  datadum <- cbind(data, height, relmot)

#gap filling
  for(i in 1:ncol(datadum)) {
    dummy <- approx(tstamp, datadum[,i], xout = tstamp)[[2]]	#linearly interpolate gaps
#    dummy <- dummy[2:(length(dummy)-1)]				#cut out potentially remaining NAs at the end points
    if (i == 1) {
      spec.in <- as.ts(dummy, frequency = fre)
    } else {
      spec.in <- ts.union(spec.in, as.ts(dummy, frequency = fre))
    }
  }

#cut out potentially remaining NAs at the end points
  spec.in <- ts.union(as.ts(tstamp, frequency = fre), spec.in)
  spec.in <- na.omit(spec.in)
  tstamp <- spec.in[,1]
  spec.in <- spec.in[,-1]
  
#separate variables
  height <- spec.in[,(ncol(spec.in)-1)]
  relmot <- spec.in[,ncol(spec.in)]
  spec.in <- as.matrix(spec.in[,1:(ncol(spec.in)-2)])
  dimnames(spec.in)[[2]] <- dimnames(data)[[2]]

#demeaning
if(demean == T & detrend == F) {
  for(i in 1:ncol(spec.in)) {
    spec.in[,i] <- spec.in[,i] - mean(spec.in[,i], na.rm=T)
  }
}

#detrending
#http://www.stat.rutgers.edu/home/rebecka/Stat565/lab42007.pdf: Before you estimate a power spectrum, trends should be removed. Rapid falloff in the spectrum makes the estimates very unreliable (bias is a problem). The spectrum() command removes linear trends by default if you omit detrend=F.
if(detrend == T) {
  for(i in 1:ncol(spec.in)) {
    spec.in[,i] <- spec.in[,i] - lm(spec.in[,i]~tstamp)$fitted.values
    #plot(spec.in[,1]~I(tstamp[2:(length(tstamp)-1)]), type="l")
    #abline(h=0, col=4, lwd=3)
    #lines(test$fitted.values~I(tstamp[2:(length(tstamp)-1)]), col=2, lwd=3)
  }
}

#create time series
spec.in <- as.ts(spec.in, frequency = fre)

#tapering (omit for EC)
#http://www.stat.rutgers.edu/home/rebecka/Stat565/lab42007.pdf: Tapering reduces bias and makes peaks look sharper. By tapering you avoid having far away frequency peaks b leakb  into other frequencies. The taper number has to be between 0 and 0.5. This is the fraction of data points that get down-weighted in the spectrum computation. Tapering reduces bias, increases variance, decreases the correlation between variables.
if(is.numeric(taper)) {
  spec.in <- spec.taper(x=spec.in, p=taper)
}

#independent variables                                                                                                                    
  #define relative frequency according to Stull
    fr <- (as.ts(1:attributes(spec.in)$dim[1], frequency = fre) - 1) / attributes(spec.in)$dim[1]
  #only 0 < frequency < Nyqvist
    fr_whr <- which(fr > 0 & fr <= 0.5)
  #corresponding 'observation' frequency
    fo <- fr * fre
  #corresponding wavelength according to taylor's hypothesis (mean relative motion, Foken, 2008 Eq. 2.106)
    la <- mean(relmot, na.rm=T) / fo
  #corresponding normalized frequency
    fn <- mean(height, na.rm=T) / la
  #corresponding wavenumber
    ka <- (2 * pi) / la

#dependent variables
  #fft returns unnormalized univariate Fourier transform of the sequence of values in z
  #therefore divide by length of time series
  #Note that for frequencies greater than .5 the Fourier transform is just the complex conjugate of the frequencies less than .5
    spec.fft <- mvfft(spec.in) / attributes(spec.in)$dim[1]

#unfolded spectral energy
  G <- Re(spec.fft * Conj(spec.fft))

#combine results
  export <- list(
    fr_whr,	#valid subset (0 < frequency < Nyqvist)
    fr,		#relative frequency (0...1) [-]
    fo,		#observation frequency (0 < frequency < fmax) [s-1]
    la,		#wavelength [m]
    fn,		#normalized frequency (fo * height / U) [-]
    ka,		#wavenumber (2 * pi * fo / U) [m-1]
    spec.in,	#continuous input data in time space (no NAs)
    spec.fft,	#variables in frequency space; complex numbers [input units]
    G		#unfolded spectral energy in frequency space[(input units)^2]
  )
  attributes(export)$names <- c("fr_whr", "fr_rel", "fr_obs", "wl_obs", "fr_nor", "wn_obs", "TScont", "FScomplex", "FSunfold")

#export results
  return(export)

########################################################
}
########################################################



########################################################
#PLOT POWER SPECTRA OF UP TO THREE VARIABLES
########################################################
SPEC.plot <- function(name, IDE, DEP, legtext, labx=NA, laby=NA, colvec=NULL, limx=NULL, limy=NULL, setoff=NULL, powerc=NULL, inisub=c(0.05,Inf)) {
#name: (no default) output location and filename for plot
#IDE: (no default) one independent variable, e.g. frequencey, wavenumber...
#DEP: (no default) spectra of up to three dependent variables, same length as IDE
#legtext: (no default) one discription (character or expression) for each variable for the legend
#labx: (NA) description for abscissa
#laby: (NA) description for ordinate
#colvec: (NULL) colors for plotting, as many colors as dependent variables
#limx: (NULL) limits for abscissa
#limy: (NULL) limits for ordinate
#setoff: (NULL) offsets different variables from each other to improve legibility (only of powerc is set)
#powerc: (NULL) coefficient of the power law c("-5/3", "-2/3") for unweighted / weighted spectra
#inisub: (c(0.05,Inf)) frequency range of initial subrange, used to fit the power law

#defaults for color vector
  coln <- ncol(as.matrix(DEP))
  if(is.null(colvec)) colvec <- c(2:(coln+1))
  if(is.null(setoff)) setoff <- rep(1,coln)

#call graphics device
  png(filename=name, width=830*2, height=1170,pointsize=2000/50,bg = "white")
    cexvar=3; par(mfrow=c(1,3),las=0,cex.axis=cexvar*0.4,cex.lab=cexvar*0.4,font.lab=2, mar=c(4,4,2,2),mgp=c(2.8,0.8,0), family="times",lwd=cexvar*1,cex.main=cexvar*0.6)

  #plotting
    #pick variables
      ide <- IDE
      if(is.null(powerc)) dep <- as.matrix(DEP)[,1] else dep <- setoff[1] * as.matrix(DEP)[,1]

    #estimate slope for reference line
      if(!is.null(powerc)) {
	powercn <- as.integer(strsplit(powerc,"/")[[1]])
	powercn <- powercn[1] / powercn[2]
	whr <- which(ide >= inisub[1] & ide < inisub[2])
	lm_slope <- lmrob(I(dep[whr]^(1/powercn)) ~ 0 + ide[whr])
	slope <- lm_slope$coefficients
	powerl <- (slope * ide)^powercn
      }

    #actual plotting
      plot(dep ~ ide, type="l", log="xy", xaxt="n", yaxt="n", xlim=limx, ylim=limy, xlab=labx, ylab=laby, col=colvec[1])
      if(!is.null(powerc)) lines(powerl ~ ide, col=1, lty=2) #, lwd=5

    #max. two additional variables
    if(coln > 1 & coln <= 3) {
      for (i in c(2:coln)) {
	#pick next dependend variable
	  if(is.null(powerc)) dep <- as.matrix(DEP)[,i] else dep <- setoff[i] * as.matrix(DEP)[,i]

	#estimate slope for reference line
	if(!is.null(powerc)) {
	  lm_slope <- lmrob(I(dep[whr]^(1/powercn)) ~ 0 + ide[whr])
	  slope <- lm_slope$coefficients
	  powerl <- (slope * ide)^powercn
	}

	#actual plotting
	  lines(dep ~ ide, col=colvec[i])
	  if(!is.null(powerc)) lines(powerl ~ ide, col=1, lty=2) #, lwd=5
      }
    }
    eaxis(side=1,labels=NA,las=0)
    eaxis(side=2,labels=NA,las=0)
    if(is.null(powerc)) {
      legend(x="bottomleft", lty=rep(1,coln), bty="n", col=colvec, cex=cexvar*0.4, xjust = 0, yjust = 0, lwd=rep(par()$lwd,coln), legend=legtext)
    } else {
      legend(x="bottomleft", lty=c(rep(1,coln),2), bty="n", col=c(colvec, 1), cex=cexvar*0.4, xjust = 0, yjust = 0, lwd=c(rep(par()$lwd,coln),par()$lwd), legend=c(legtext, substitute(paste("f"^{powerc}, " law", sep=""), list(powerc=powerc))))
    }
    box()

#close graphics device
  dev.off()
########################################################
}
########################################################



########################################################
#GENERATE COSPECTRA
########################################################
COSP.fwd <- function(FS, TS, cparho=1) {
#FS: frequency space data, raw spectra as complex numbers
#TS: time space data, identical variables and length as FS
#cparho: (1) single factor (specific heat at constant pressure * air density) for conversion from kinematic to energetic units; only used for comparison with Reynolds-Covariance, not for output

#calculate cospectra
  for(i in 2:ncol(FS)) {
    if(i == 2) {
      COSPdata <- Re(FS[,1]) * Re(FS[,i]) + Im(FS[,1]) * Im(FS[,i])
    } else {
      COSPdata <- cbind(COSPdata, Re(FS[,1]) * Re(FS[,i]) + Im(FS[,1]) * Im(FS[,i]))
    }
  }
  COSPdata <- as.matrix(COSPdata)
  if(ncol(COSPdata) == 1) {
    dimnames(COSPdata)[[2]] <- list(paste("Co(w'", dimnames(FS)[[2]][2:ncol(FS)], "')", sep=""))
  } else {
    dimnames(COSPdata)[[2]] <- paste("Co(w'", dimnames(FS)[[2]][2:ncol(FS)], "')", sep="")
  }

#the sum should be the same as the covariance:
  COSPsum <- sapply(1:ncol(COSPdata), function(x) sum(COSPdata[2:nrow(COSPdata),x], na.rm=T)) * cparho
  REYsum <- stats:::cov(TS)[1,c(2:ncol(TS))] * cparho
    names(COSPsum) <- names(REYsum) <- dimnames(COSPdata)[[2]]
  #relative difference
    DIFFperc <- (COSPsum - REYsum) / REYsum * 100

#generate output
  OUTsum <- rbind(COSPsum, REYsum, DIFFperc)
  export <- list(COSPdata, OUTsum)
    attributes(export)$names <- c("FScosp", "sum")

#return results
  return(export)
########################################################
}
########################################################



########################################################
#PLOT COSPECTRA
########################################################
COSP.plot <- function(
  name,
  IDE,
  DEP,
  title,
  labx=NA,
  laby=NA,
  colvec=c(1,2,1),
  limx=NULL,
  limy=NULL,
  weight=T,
  wl=F,
  norm=NULL,
  part_cov=1,
  modspec=T,
  zeta,
  FX=NULL,
  devran=NULL,
  plotting=T,
  peak_shift=NULL,
  devran_plot=NULL,
  OUT_plotting
  ) {
#name: (no default) output location and filename for plot
#IDE: (no default) one independent variable, e.g. frequencey, wavenumber...
#DEP: (no default) spectra of up to three dependent variables, same length as IDE
#title: (no default) one discription (character or expression) for each variable for the title
#labx: (NA) description for abscissa
#laby: (NA) description for ordinate
#colvec: (c(1,2,1)) colors for plotting of measurement, model and evaluation range
#limx: (NULL) limits for abscissa
#limy: (NULL) limits for ordinate
#weight (T) weight the spectra for independent variable (wavelength / frequency)?
#wl: (F) is wavelength used as independent variable?
#norm: (NULL) normalizes sum(cospectrum) to a given quantity, e.g. unity or the total flux per variable, else NULL
#part_cov: (1) allows to plot measured cospectrum with scaling factor that is different from norm
#modspec: (T) shall model cospetra be calculated and plotted?
#zeta: (no default) atmospheric stability to choose the reference cospectrum
#FX: (NULL) frequency at which fCO(f) reaches its maximum value; will be determined from measured data if NULL
#devran: (NULL) range to assess sum((cospectrum - model)[whr]) / sum(model) [%]; only valid for equidistant data
#peak_shift: peak peak_shift to make observation peak and peak in theory match to each other
#devran_plot: whether we would like to draw the limit line
#OUT_plotting: determine which kind of result to export

#prepare variables
  ide <- IDE
  fs <- as.matrix(DEP)

#weighting of the cospectrum: divide for wavelenght, multiply for frequencies
  if(weight == T) {
    if(wl==T) fs <- as.matrix(fs / ide) else fs <- as.matrix(fs * ide)
  }


#call graphics device
  if(plotting == T) {
    
    png(filename=name, width = 1000, height = 1000, units = "px", pointsize = 20, bg = "white")
    cexvar=3; par(mfrow=c(2,2), las=0, cex.axis=cexvar*0.4, cex.lab=cexvar*0.4, font.main=1, font.lab=1,
      mar=c(4,4,2,2), mgp=c(2.6,0.8,0), family="times", lwd=cexvar, cex.main=cexvar*0.5)
    
  }
    
#normalize each cospectrum
if(!is.null(norm)) {
  if(length(norm) == 1) norm <- rep(norm, ncol(fs))
  FSsum <- sapply(1:ncol(fs), function(x) sum(fs[,x], na.rm=TRUE))
  fs <- sapply(1:length(FSsum), function(x) as.matrix(fs[,x] / FSsum[x] * norm[x]))
}

##loop around variables
##
i <- 1
  #Model cospectrum after Massman, 2005 (in Lee, 2005); continuous approximation of the Kaimal (1972) cospectra
    if(modspec == T) {
      #(inertial subrange) slope parameter; 3/4 for -7/3 (cospectra), and 3/2 for -5/3 (spectra) power law
      	m=3/4
      #broadness parameter; 1/2 for unstable, 7/6 for stable stratification 
      	if(zeta <= 0) mue=1/2  else mue=7/6
      #frequency at which fCO(f) reaches its maximum value in theory
      	if(!is.null(FX)) fx=FX else fx=ide[which(fs[,i] == max(fs[,i],na.rm=TRUE))]
     
      #frequency at which fCO(f) reaches its maximum value for observation data
      if (is.null(peak_shift)){
        
        fx_obs <- ide[which(fs[,i] == max(fs[,i],na.rm=TRUE))]
        
        peak_shift <- fx_obs / fx
      }
      
#       if(!is.null(devran)){
#         if(peak_shift <= devran[1] / freq){
#           #if the part we shift is smaller than the range we calculate, then we set peak_shift = 1, flag = 1
#           peak_shift <- 1
#           flag <- 1
#         } else {
#           flag <- 0
#         }
#       }
      
      ide_obs <- ide
      ide_mod <- ide_obs

#       where <- which(ide_mod < 1/3600 | ide_mod > freq / 2)
#       if (length(where) > 0){
#         ide_mod <- ide_mod[-where]
#         fs <- fs[-where,i]
#       } else {
#         ide_obs <- ide
#         fs <- fs[,i]
#       }
# 
#      ide_mod <- ide_obs         

      #normalize each cospectrum
      if(!is.null(norm)) {
        #if(length(norm) == 1) norm <- rep(norm, ncol(fs))
        FSsum <- sum(fs, na.rm=TRUE)
        fs <- fs / FSsum * norm
      }
      
      
      #calculate non-scaled, frequency-weighted model Cospectrum (fCo or nCo)
      	COmM <- ((ide_mod / peak_shift) / fx) / (
      	  ( 1 + m * ((ide_mod / peak_shift) / fx)^(2 * mue) )^( (1/(2*mue)) * ((m+1)/m) )
      	)
		if(weight == F) COmM <- COmM / ide_mod 
      #individual normalisation parameter for each variable
      	A0 <- sum(fs, na.rm=T) / sum(COmM, na.rm=T)
      	COmM <- A0 * COmM
    }

  #Assessment of deviation(cospectrum - model) in given frequency range
    #only valid for equidistantly binned data; else multiply with bin width
    if(modspec == T & !is.null(devran)) {
      #range to assess    
        whr <- which(ide_obs >= devran[1] & ide_obs <= devran[2])
      #weighted or unweighted?
        if(weight == T) {
          
          if(wl==T) {
            fac <- 1 - ( sum(((fs - COmM) * ide_obs)[whr], na.rm=T) / abs( sum(fs * ide_obs, na.rm=T) ) )
          } else {
            fac <- 1 - ( sum(((fs - COmM) / ide_obs)[whr], na.rm=T) / abs( sum(fs / ide_obs, na.rm=T) ) )
          }
        } else {

          crit <- 1
          crit_run <- 0
           fac_0 <- 1
          
          while(crit > 1e-2) {
            crit_run <- crit_run + 1
            
            #fac_dum <- 1 - ( sum((fs[,i] - COmM * fac_0)[whr], na.rm=T) / sum(fs[,i], na.rm=T) )
			fac_dum <- 1 - ( sum((fs / fac_0 - COmM)[whr], na.rm=T) / sum(fs, na.rm=T) )
                      
            crit <- fac_dum - fac_0
            fac_0 <- fac_dum
                        
          }
            fac <- fac_dum
          
        }
      #store      
        if(i == 1) corfac <- fac else corfac <- c(corfac, fac)
    }

  #plot cospectrum  
  
    if(plotting == T) {
      
      plot(fs * part_cov ~ ide_obs, type="l", log="x", xaxt="n", yaxt="n", xlim=limx, ylim=range(COmM,(fs * part_cov),na.rm=TRUE), main=title[i], xlab=labx,
        ylab=laby, col=colvec[1])
      
        abline(h=0, col="grey")
      
        if(modspec == T) lines(COmM ~ ide_mod, lty=2, col=colvec[2], lwd=5)
        if(modspec == T & !is.null(devran_plot)) abline(v=c(devran_plot), lty=2, col=colvec[3])
      
        eaxis(side=1,labels=NA,las=0)
        eaxis(side=2,labels=NA,las=0)
      
        legend(x="topleft", lty=c(1,2,2,1,1), bty="n", col=c(colvec,NA,NA), xjust = 0, yjust = 0,
          lwd=c(par()$lwd,5,par()$lwd,1,1), legend = c(
          	"measured",
          	"Massman (2005)",
          	"evaluated range",
          	paste("stability", " = ", round(zeta,1), sep=""),
          	paste("N = ", length(ide), sep="")
        ))
      
        box()
    }

##end loop around variables

##


#close graphics device
  if(plotting == T) {
    dev.off()
  }

#store return variables


#export results
if (OUT_plotting == FALSE){
  export_FALSE <- list(
    corfac
  )
  
  rm(peak_shift)
  
  return(export_FALSE)
}

if(OUT_plotting == TRUE){
  export_TRUE <- list(
    fs,
    COmM,
    peak_shift)
  
  rm(peak_shift)
  
  return(export_TRUE)  
}

########################################################
}
########################################################



########################################################
#PLOT OGIVES
##cumulation only valid with raw data or equidistant binned data!
########################################################
OGIV.plot <- function(name, IDR, RAW, IDS, SMO, thresh, Fperi, title=NULL, labx=NULL, laby=NULL, colvec=c("grey",1,1), LIMX=NULL, LIMY=NULL, weight=F, wl=T, norm=1) {
#name: (no default) output location and filename for plot
#IDR: (no default) one independent variable for the raw cospectra, e.g. frequencey, wavenumber...
#RAW: (no default) raw cospectra of up to four dependent variables, same length as IDR
#IDS: (no default) one independent variable for the smoothed cospectra, e.g. frequencey, wavenumber...
#SMO: (no default) smoothed cospectra of up to four dependent variables, same length as IDR
#thresh: (no default) flux contribution level to evaluate e.g. corresponding wavelenght
#Fperi: (no default) flux averaging period to evaluate corresponding contribution level
#title: (no default) one discription (character or expression) for each variable for the title
#labx: (NA) description for abscissa
#laby: (NA) description for ordinate
#colvec: (c(1,2,1)) colors for plotting of measurement, model and evaluation range
#LIMX: (NULL) limits for abscissa
#LIMY: (NULL) limits for ordinate
#weight (F) weight the cospectra for independent variable (wavelength / frequency)?
#wl: (T) is wavelength used as independent variable?
#norm: (1) normalizes max(cospectrum) to a given quantity, e.g. total flux per variable, else unity

#prepare local variables
  idr <- IDR
  raw <- as.matrix(RAW)
  ids <- IDS
  smo <- as.matrix(SMO)

#sort for increasing wavelenght (a) or decreasing frequency (b) -> always from short- to longer wavelength contributions
  if(wl == T) {
    raw <- as.matrix(raw[sort.list(idr, decreasing = F),])	#first dependent variable
    idr <- idr[sort.list(idr, decreasing = F)]			#then independent variable
    smo <- as.matrix(smo[sort.list(ids, decreasing = F),])	#first dependent variable
    ids <- ids[sort.list(ids, decreasing = F)]			#then independent variable
  } else {
    raw <- as.matrix(raw[sort.list(idr, decreasing = T),])
    idr <- idr[sort.list(idr, decreasing = T)]
    smo <- as.matrix(smo[sort.list(ids, decreasing = T),])
    ids <- ids[sort.list(ids, decreasing = T)]
  }

#weighting of the cospectrum: divide for wavelength, multiply for frequencies
#only valid for cospectrum, not for Ogive!
  if(weight == T) {
    if(wl==T) {
#      raw <- as.matrix(raw / idr)
      smo <- as.matrix(smo / ids)
    } else {
#      raw <- as.matrix(raw * idr)
      smo <- as.matrix(smo * ids)
  }}

#make sure that all fluxes are cumulating towards positive
    signum <- sign(colSums(raw))	#check the sign(tendency) of the different Cospectra
    smo <- t(signum * t(smo))
    raw <- t(signum * t(raw))

#normalize the smoothed cospectrum
  if(!is.null(norm)) {
    if(length(norm) == 1) norm <- rep(norm, ncol(raw))
    smo <- sapply(1:length(norm), function(x) as.matrix(smo[,x] / max(smo[,x]) * norm[x]))
  }

#call graphics device
  png(filename=name, width = 1000, height = 1000, units = "px", pointsize = 20, bg = "white")
  cexvar=3; par(mfrow=c(2,2), las=0, cex.axis=cexvar*0.4, cex.lab=cexvar*0.4, font.lab=2, mar=c(4,4,2,2), mgp=c(2.6,0.8,0), family="times", lwd=cexvar, cex.main=cexvar*0.6)


##loop around variables
for(i in 1:ncol(raw)) {
##

  #Ogive: cumulated flux contribution
    OGIVE <- cumsum(raw[,i])
    ogive <- OGIVE / max(OGIVE) * norm[i]

  #determine wavelength for given flux contribution level
    #find first value exceeding given contribution level
      zcross <- zeroCross(ogive - thresh, slope="positive") + 1
    #flag whether zero crossing is singular (s) or multiple (m)
      if(length(zcross == 1)) Og.flag <- "s" else Og.flag <- "m"
    #wavelength corresponding to given contribution level
      wathr <- idr[(zcross[1]-1)]

  #determine flux contribution level for given wavelength
    Fnear <- nearest(idr, Fperi)
    Fleve <- ogive[Fnear]

  #store ogive results
    OGdum <- c(wathr, Fleve)
    if(i == 1) OGout <- OGdum else OGout <- cbind(OGout, OGdum)

  #plotting range
    if(!is.null(LIMX)) limx <- LIMX else limx <- range(idr)
    if(!is.null(LIMY)) limy <- LIMY else limy <- range(c(smo[,i], ogive))

  #plotting
    plot(smo[,i] ~ ids, col=colvec[1], type="l", xlim=limx, ylim=limy, main=title[i], xlab=labx, ylab=laby, las=1, log="x", xaxt = "n")
      lines(ogive ~ idr, col=colvec[2])
      abline(h=0)
      abline(h=thresh, lty=2, col=colvec[3])
      abline(v=wathr, lty=2, col=colvec[3])
      legend(x="topleft", lty=c(1,1,1,1,2), bty="n", col=c(NA,NA,col=colvec,NA,NA,NA,NA), cex=cexvar*0.4, pt.cex=cexvar/2, pt.lwd=cexvar*2/3, xjust = 0, yjust = 0, legend = c(
	"","",
	"Cospectrum",
	"Ogive (Og)",
	paste(thresh, " max(Og)", sep=""),
	paste("== ", round(wathr,2)," X", sep="")
      ))
      legend(x="topright", bty="n", cex=cexvar*0.4, pt.cex=cexvar/2, xjust = 0, yjust = 0, legend = c(
	"","","","",
	paste(Fperi, " X ==", sep=""),
	paste(round(Fleve,2), " max(Og)", sep="")
      ))
      eaxis(side=1,labels=NA,las=0)
      box()

##end loop around variables
}
##


#close graphics device
  dev.off()

#export Ogive results
  dimnames(OGout) <- list(c(paste(thresh, "(Og)==", sep=""), paste(Fperi, " X==", sep="")), dimnames(raw)[[2]])
  return(OGout)

########################################################
}
########################################################