###########################################
# Collection of re-usable functions
# include into your scripts
# at the top (but after rm(list=ls()))
# source("myfunctions.R")
# #########################################



###########################################
# Basics & Examples How to code in R
###########################################

# the function that does nothing
empty<-function(){
  
}

# what does a well-trained parrot do?
parrot<-function(text) {
  return(text)  
}
# this requires library "animation"
fun<-function(){
  if(require(animation)){
    oops<-ani.options(interval=0.1,nmax=10)
     for ( i in 1:100) {
        y<-sin(i*pi/10)
        plot(c(1,100),c(-1,1),typ='n')
        points(i-1,y,pch='o',col="blue")
        points(i+1,y,pch='o',col="blue")
        points(i,y-0.1,pch='|',col="blue")
        points(i,y-0.2,pch='~',col="blue")
        
        ani.pause()
    }
    ani.options(oops)
  }
}


# loading a gray-shaded image of ASCII PPM format
loadppm<-function(file="data/labrador.ppm",n=128,m=105){
  res<-scan(file,what="",n=10)
  print(res)
  # NOTE ppm images are stored line by line
  # (columns filled first then rows)
  n<-as.numeric(res[9])
  m<-as.numeric(res[10])
  buffer<-read.table(file,skip=4)
  M<-matrix(0,n,m,byrow=FALSE)
  ir<-seq(1,3*n*m,3)
  ig<-seq(2,3*n*m,3)
  ib<-seq(3,3*n*m,3)
  R<-buffer[ir,]
  G<-buffer[ig,]
  B<-buffer[ib,]
  shade<-1/3*(R+G+B)
  M[]<-shade
  # use only the average of the 3three RGB
  return(t(M))
}


matrix_flip_v<-function(X){
  # using image, contour or filled.contour
  # R plots the data matrix up-side-down
  # i.e. the first row is at the bottom
  # here I flip the data matrix X vertically
  isize<-dim(X)
  nrow<-isize[1]
  ncol<-isize[2]
  # define the vertical flip matrix operation
  # M X 
  M<-matrix(0,nrow,nrow)
  for ( i in 1:nrow) {
    j<-nrow-i+1
    M[i,j]<-1
  }
  X<-M%*%X
  return(X)
}

###########################################################
# Plotting support functions
###########################################################


matrix_image<-function(z,x=NULL,y=NULL,xlab="column",ylab="row (-1)",main="matrix image",col=NULL){
  # plots a matrix z the correct way
  # columns represent columns, and rows are
  # rows in the right order as a print(X) would print
  # the data entries of the matrix X
  # NOTE: col expects a color palette
  nrow<-dim(z)[1]
  ncol<-dim(z)[2]
  if (is.null(x)) {x<-seq(1,ncol,1)}
  if (is.null(y)) {y<-(-1)*rev(seq(1,nrow,1))}
  buffer<-t(matrix_flip_v(z))
  if (is.null(col)) {
    image(x=x,y=y,z=buffer,xlab=xlab,ylab=ylab,main=main)
  } else {
    image(x=x,y=y,z=buffer,col=col,xlab=xlab,ylab=ylab,main=main)
  }    
  # return the buffer used for image plotting
  return(buffer)
}


plot_vector<-function(v,col=1,lwd=1,lty=1,scale=1,xoff=0,yoff=0){
  # by default start line from origin (0,0)
  # offset xoff and yoff are allowed
  # arrowsize per default 10% of length
  # at 45 degree angle
  x<-v[1]
  y<-v[2]
  alpha<-30/180*pi
  arrowsize<-0.1
  len<-sqrt(x^2+y^2)
  phi<-atan2(y,x)
  dlen<-arrowsize*len*scale
  dx<-dlen*cos(alpha)
  dy<-dlen*sin(alpha)
  lines(c(xoff,xoff+x),c(yoff,yoff+y),col=col,lty=lty,lwd=lwd)  
  # rotation matrix
  M<-matrix(c(cos(phi),sin(phi),-1*sin(phi),cos(phi)),2,2)
  x1<-len-dx
  x2<-x1
  y1<-dy
  y2<-(-1)*dy
  v1<-M%*%c(x1,y1)
  v2<-M%*%c(x2,y2)
  #print(paste(len,dlen,dx,dy,alpha/pi*180))

  lines(c(xoff+x,xoff+v1[1]),c(yoff+y,yoff+v1[2]),col=col,lwd=lwd,lty=lty)
  lines(c(xoff+x,xoff+v2[1]),c(yoff+y,yoff+v2[2]),col=col,lwd=lwd,lty=lty)
}

myplot<-function(xrange=c(0,1),yrange=c(0,1),frame=FALSE,xlab="x",ylab="y",main=""){
  cex<-1
  scale<-sqrt(xrange[2]^2+yrange[2]^2)
  scale<-0.2*scale
  plot(xrange,yrange,typ='n',frame=frame)
  plot_vector(c(0,yrange[2]),scale=scale)
  text(xrange[2],-0.75,xlab)
  plot_vector(c(xrange[2],0),scale=scale)
  text(-0.35,yrange[2],ylab) 
  lines(xrange,c(0,0))
  lines(c(0,0),yrange)
  
}


#####################################################
# Section: Physics, Meteorology, Climatology  
#####################################################


# Mangus formula for water vapor saturation pressure
# as function of Temperature
magnus<-function(t){
  # t: temperature in C
  # returns pressure in hPa
  z<-(17.625*t)/(t+243.04)
  res<-6.1094*exp(z)
  return(res)
}

c2f<-function(x){
  # Fahrenheit to Celsius conversion function
  # x is local variable in the function
  # that gets the object from the function call
  # (e.g. a number, or a vector of numbers)
  res<-x*9/5+32
  return(res)
}

f2c<-function(x){
  # Celsius to Fahrenheit conversion function
  res<-(x-32)*5/9
  return(res)
}

############################################
# Statistical Calculations
############################################





center<-function(x){
  # center a data sample
  res<-x-mean(x,na.rm=TRUE)
  return(res)
}


standardize<-function(x){
  # center and then scale to variance = 1
  res<-center(x)
  xsd<-sd(x,na.rm=TRUE)
  res<-res/xsd
  return(res)
}


anomaly<-function(x,clim){
  ##########################################################
  # function to remove the climatological seasonal cycle
  # x: input vector with monthly mean data
  # clim: input vector with the seasonal cycle (12 months)
  #########################################################
  inum<-length(x)
    icyc<-length(clim)
    # counter for the index position in clim
    m<-0
    print(paste(inum,icyc))
  
    for (i in 1:inum) {
      m<-m+1
      if (m>12) {m<-1}
      x[i]<-x[i]-clim[m]
    }  
    return(x)
}

# here I repeat the climatology
# and use then vector subtraction
anomaly2<-function(x,clim){
  inum<-length(x)
  icyc<-length(clim)
  irep<- inum/icyc
  buffer<-rep(clim,irep)
  res<-x-buffer
  return(res)
}


  

#2014-02-24 Gaussian distribution

fgauss<-function(x,center=0,s=1){
  # if nothing is provided for
  # center and standard deviation
  # it is used 0 and 1 respectively
  res<-1/(s*sqrt(2*pi)) * exp( -1*(x-center)^2/(2*s^2))
  return(res)
}



# test functions

t.test.cor<-function(r,n,p){
  # assume gaussian data 
  # independent samples, or n effective sample size
  if (!is.na(r)) {
    t<-r/sqrt(1-r^2)*sqrt(n-2)
    # two-sided test
    p.value<-pt(t,df=n-2)
    if (p.value<p/2. | p.value > (1-p/2.)) {res<-1} else {res<-0}
  } else {
    res<-NA
  }
  return(res)
}