########################### Library_Rcran-Switzerland_EPFL.R ###########################
# This library provides useful R-cran functions (see 'http://cran.r-project.org/' for additional info on the R-cran statistical software) associated with the 'Switzerland_EPFL' data set hosted by the NASA. The user is refered to the associated documentation ('EPFL-LTE_Network_disdro-NASA.pdf') for additional information on the data set and the different function presented here.

# Environmental Remote Sensing Lab., EPFL
# March 2012.


######################## Network Station Coordinates ########################
get.network_coordinates <- function(id_station,all=FALSE){

    ## Returns the coordinates (lat,long,alt,estY,nordX) of EPFL-LTE network of 16 disdrometers that have been deployed in Lausanne, Switzerland.
    
    ## Input: 
    ## id_station = a list with the desired stations (caution: id = 10,11,12,13,20,21,...)
    ## all = logical, if TRUE the coordinates of all stations available for this campaign are returned
    
    ## Output:
    ## CoordM = a matrix of station coordinates (id;lat;long;alt;Est;Nord), 1 row per station
    
    ## Source: M.Schleiss, October 2009
    ## Modified by J. Jaffrain, Oct 12th 2010.
    
    ## Remarks:
    ## All coordinates were measured with Garmin GPS Dakota 20. 
    ## EstY and NordX coordinates were obtained using the online NAVREF
    ## projection tool provided by the swiss topographic institute 
    ## http://www.swisstopo.admin.ch/internet/swisstopo/fr/home/apps/calc/navref.html
    ## All altitudes are assumed constant at 400m.

    TotalCoordM <- matrix(NA,nrow=16,ncol=6)
    TotalCoordM[1,]  <- c(10,46.520500,6.565200,400,532977,152507)
    TotalCoordM[2,]  <- c(11,46.520433,6.562833,400,532795,152502)
    TotalCoordM[3,]  <- c(12,46.521900,6.565183,400,532977,152663)
    TotalCoordM[4,]  <- c(13,46.521267,6.566767,400,533098,152591)

    TotalCoordM[5,]  <- c(20,46.519800,6.570500,400,533383,152425)
    TotalCoordM[6,]  <- c(21,46.519583,6.572317,400,533522,152399)
    TotalCoordM[7,]  <- c(22,46.521200,6.572583,400,533544,152579)
    TotalCoordM[8,]  <- c(23,46.520533,6.571100,400,533429,152506)

    TotalCoordM[9,]  <- c(30,46.518333,6.563933,400,532877,152267)
    TotalCoordM[10,] <- c(31,46.519650,6.563900,400,532876,152414)
    TotalCoordM[11,] <- c(32,46.518700,6.562733,400,532785,152309)
    TotalCoordM[12,] <- c(33,46.517633,6.564583,400,532926,152189)

    TotalCoordM[13,] <- c(40,46.521017,6.569733,400,533325,152561)
    TotalCoordM[14,] <- c(41,46.519500,6.567883,400,533181,152394)
    TotalCoordM[15,] <- c(42,46.520600,6.567850,400,533180,152516)
    TotalCoordM[16,] <- c(43,46.521400,6.567867,400,533182,152605)


    if(all==TRUE){CoordM <- TotalCoordM}
    if(all==FALSE){
        Nstations <- length(id_station)
	if(Nstations==0){stop("no station has been specified")}
	CoordM <- data.frame(matrix(NA,nrow=Nstations,ncol=6))
	for(i in 1:Nstations){
	    id <- id_station[i]
	    j  <- which(TotalCoordM[,1]==id)
	    Nj <- length(j)
	    if(Nj!=1){print(sprintf("warning: station %i not found",id))}
	    if(Nj==1){CoordM[i,] <- TotalCoordM[j,]}
	}
    }
    for(j in 2:6){CoordM[,j] <- as.numeric(CoordM[,j])}
    return(CoordM)
}

##################### Distance Between Network Stations #####################
get.network_station_distance <- function(id_station1,id_station2){

    ## Returns the distance (in meters) between two LTE network stations
    ## If the two stations are identical, the returned distance is zero
    ## If one of the stations does not exist, the function returns NA
    
    ## Input:
    ## id_station1 = number of the first station
    ## id_station2 = number of the second station
    
    ## Output:
    ## dist = the distance (in meters) between the two stations

    rt <- 6378100

    if(is.na(id_station1) || is.na(id_station2)){
	print("Error: NA in station number")
	stop()
    }
    if(id_station1==id_station2){dist <- 0}
    else{
	coords <- get.network_coordinates(c(id_station1,id_station2))
	lat1   <- coords[1,2]*pi/180
	long1  <- coords[1,3]*pi/180
	lat2   <- coords[2,2]*pi/180
	long2  <- coords[2,3]*pi/180
	dist   <- rt*acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(long2-long1))
    }
    return(dist)  
}

########################## Parsivel_classes ##############################
# This function provides the middle and the spread of Parsivel's 32 velocity and diameter classes.
Parsivel_classes <- function(test=0){
	## Equivolumetric drop diameter classes in mm ###
	class_size <- vector(mode = "numeric", length = 32)
	class_size[1] = 0.062
	class_size[2] = 0.187
	class_size[3] = 0.312
	class_size[4] = 0.437
	class_size[5] = 0.562
	class_size[6] = 0.687
	class_size[7] = 0.812
	class_size[8] = 0.937
	class_size[9] = 1.062
	class_size[10] = 1.187
	class_size[11] = 1.375
	class_size[12] = 1.625
	class_size[13] = 1.875
	class_size[14] = 2.125
	class_size[15] = 2.375
	class_size[16] = 2.750
	class_size[17] = 3.250
	class_size[18] = 3.750
	class_size[19] = 4.250
	class_size[20] = 4.750
	class_size[21] = 5.500
	class_size[22] = 6.500
	class_size[23] = 7.500
	class_size[24] = 8.500
	class_size[25] = 9.500
	class_size[26] = 11.000
	class_size[27] = 13.000
	class_size[28] = 15.000
	class_size[29] = 17.000
	class_size[30] = 19.000
	class_size[31] = 21.500
	class_size[32] = 24.500
	
	## Classes Spread in mm
	class_size_spread <- c(rep(0.125,10),rep(0.250,5),rep(0.500,5),rep(1.000,5),rep(2.000,5),rep(3.000,2))

	## Fall velocity classes in m.s^-1 ###
	class_speed <- vector(mode = "numeric", length = 32)
	class_speed[1] = 0.050
	class_speed[2] = 0.150
	class_speed[3] = 0.250
	class_speed[4] = 0.350
	class_speed[5] = 0.450
	class_speed[6] = 0.550
	class_speed[7] = 0.650
	class_speed[8] = 0.750
	class_speed[9] = 0.850
	class_speed[10] = 0.950
	class_speed[11] = 1.100
	class_speed[12] = 1.300
	class_speed[13] = 1.500
	class_speed[14] = 1.700
	class_speed[15] = 1.900
	class_speed[16] = 2.200
	class_speed[17] = 2.600
	class_speed[18] = 3.000
	class_speed[19] = 3.400
	class_speed[20] = 3.800
	class_speed[21] = 4.400
	class_speed[22] = 5.200
	class_speed[23] = 6.000
	class_speed[24] = 6.800
	class_speed[25] = 7.600
	class_speed[26] = 8.800
	class_speed[27] = 10.400
	class_speed[28] = 12.000
	class_speed[29] = 13.600
	class_speed[30] = 15.200
	class_speed[31] = 17.600
	class_speed[32] = 20.800
    
    list <- list(class_size,class_speed,class_size_spread)
    names(list) <- c("class_size","class_speed","class_size_spread")
    return(list)
}

######################### Uncertainty associated with Parsivel measurements ################################
uncertainty_Parsiv <- function(vect,var_type,dt,freq=9.4,fill.na=FALSE){
	## This function was written in order to give the sampling uncertainty values associated with Parsivel measurements quantified by Jaffrain and Berne [JHM 2011].
	##
	## Inputs:
	##    'vect': the vector of value to calculate the uncertainty (ex: rain rates, reflectivities,...)
	##    'var_type': a character vector with the name of the variable of interest (choice between: 'Nt', 'D0','Dmw', 'R', 'Z' or 'Zdr')
	##    'dt': temporal resolution [in sec]
	## Optional arguments:
	##    'freq': numeric value with the radar frequency to consider (9.4, 5.6 or 2.8).
	##    'fill.na': logical specifying if NA values of uncertainty should be set to the last uncertainty values recorded at this temporal resolution (only available for Nt, R, Zh and kdp).
	## Outputs:
	##    'uncer': a matrix with 2 columns: 1st= relative uncertainty values (%), 2nd=rel* input vect (absolute uncertainty)

	# Uncertainty values are calculated using Parsivel 01 and 02 during the 'Roof_2008' campaign (15 months of collocated emasurements).
	# Moments were calculated from raw DSD measurements after velocity-size relationship filtering ('Kruger_2002' filtering with + or - 60% tolerance).
	# Moreover, uncertainty is quantified as std calculated from interquantiles ((Q90-Q10)/2*1.28).
	# Uncertainty is quantified from at least 30 values.
	# For additional info, see 'Jaffrain_JHM_2011'.

	# Written by Joel Jaffrain, May 2010
	# Last update: November, 15th 2010

  if(is.null(dim(vect)) == FALSE){stop("The input should be a vector !")}
  if(sum(round(freq,digits=1) == c(9.4,5.6,2.8))==0){stop("Provide a useful frequency band: 9.4 GHz(X), 5.6 GHz(C) or 2.8 GHz(S)")}
  
  list_tres <- c(20,60,120,180,240,300,600,900,1800,2700,3600)
  n_tres <- length(list_tres)
  pos_dt <- which(dt == list_tres)
  uncer <- matrix(NA,ncol=2,nrow=length(vect))
  n_meas <- length(vect)

  if(var_type == "Nt"){
    #### Nt ####
    class_inf <- c(0,50,100,200,400,600,800,1000,1500,2000)
    mat <- matrix(NA,nrow=9,ncol=n_tres)
    mat[1,] <- c(0.25,0.17,0.13,0.12,0.11,0.09,0.07,0.07,0.06,0.05,0.04)
    mat[2,] <- c(0.22,0.14,0.11,0.09,0.08,0.08,0.07,0.06,0.05,0.06,0.05)
    mat[3,] <- c(0.17,0.11,0.09,0.07,0.07,0.07,0.06,0.06,0.05,0.05,0.05)
    mat[4,] <- c(0.13,0.08,0.07,0.06,0.06,0.06,0.05,0.04,0.04,0.04,0.04)
    mat[5,] <- c(0.10,0.07,0.05,0.05,0.05,0.04,0.04,0.04,0.03,0.04,NA)
    mat[6,] <- c(0.09,0.06,0.05,0.05,0.04,0.04,0.04,NA,NA,NA,NA)
    mat[7,] <- c(0.08,0.06,0.05,0.05,0.04,NA,NA,NA,NA,NA,NA)
    mat[8,] <- c(0.07,0.05,0.05,0.04,NA,NA,NA,NA,NA,NA,NA)
    mat[9,] <- c(0.07,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
  }else if(var_type == "D0"){
    #### D0 ####
    class_inf <- c(0.6,0.7,0.8,0.9,1,1.25,1.5,1.75,2,2.25,2.5)
    mat <- matrix(NA,nrow=10,ncol=n_tres)
    mat[1,] <- c(0.05,0.03,0.03,0.02,0.02,0.02,0.02,0.02,0.01,NA,NA)
    mat[2,] <- c(0.06,0.04,0.03,0.03,0.03,0.02,0.02,0.02,0.02,0.02,0.02)
    mat[3,] <- c(0.06,0.04,0.03,0.03,0.03,0.03,0.02,0.02,0.02,0.02,0.02)
    mat[4,] <- c(0.08,0.05,0.04,0.04,0.04,0.03,0.03,0.03,0.03,0.03,0.03)
    mat[5,] <- c(0.10,0.06,0.05,0.05,0.04,0.04,0.04,0.04,0.04,0.04,0.03)
    mat[6,] <- c(0.11,0.07,0.06,0.05,0.05,0.04,0.04,0.04,0.04,0.03,0.03)
    mat[7,] <- c(0.15,0.10,0.08,0.06,0.07,0.06,0.05,0.06,0.03,0.03,NA)
    mat[8,] <- c(0.19,0.10,0.09,0.08,0.07,0.06,0.06,NA,NA,NA,NA)
    mat[9,] <- c(0.22,0.14,0.10,0.12,0.07,NA,NA,NA,NA,NA,NA)
    mat[10,] <- c(0.27,0.15,NA,NA,NA,NA,NA,NA,NA,NA,NA)
  }else if(var_type == "Dmw"){
    #### Dmw ####
    class_inf <- c(0.6,0.7,0.8,0.9,1,1.25,1.5,1.75,2,2.25,2.5,3)
    mat <- matrix(NA,nrow=11,ncol=n_tres)
    mat[1,] <- c(0.04,0.03,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
    mat[2,] <- c(0.05,0.03,0.03,0.02,0.02,0.02,0.02,0.02,0.01,0.01,  NA)
    mat[3,] <- c(0.06,0.04,0.03,0.03,0.03,0.03,0.02,0.02,0.02,0.02,0.02)
    mat[4,] <- c(0.07,0.05,0.04,0.04,0.03,0.03,0.03,0.03,0.03,0.02,0.02)
    mat[5,] <- c(0.09,0.06,0.05,0.05,0.04,0.04,0.04,0.03,0.03,0.03,0.03)
    mat[6,] <- c(0.12,0.08,0.06,0.06,0.05,0.05,0.04,0.04,0.04,0.04,0.03)
    mat[7,] <- c(0.14,0.09,0.07,0.06,0.06,0.06,0.05,0.05,0.04,0.05,0.04)
    mat[8,] <- c(0.17,0.11,0.08,0.07,0.07,0.06,0.05,0.05,0.04,  NA,  NA)
    mat[9,] <- c(0.19,0.13,0.12,0.11,0.11,0.08,0.09,0.07,  NA,  NA,  NA)
    mat[10,] <- c(0.23,0.17,0.15,0.14,0.10,  NA,  NA,  NA,  NA,  NA,  NA)
    mat[11,] <- c(0.31,0.19,0.14,0.13,0.10,  NA,  NA,  NA,  NA,  NA,  NA)
  }else if(var_type == "R"){
    #### Rdsd ####
    class_inf <- c(0.1,2,4,6,8,10,15,20,25,30)
    mat <- matrix(NA,nrow=9,ncol=n_tres)
    mat[1,] <- c(0.25,0.18,0.15,0.13,0.12,0.12,0.11,0.10,0.10,0.10,0.09)
    mat[2,] <- c(0.21,0.15,0.12,0.11,0.11,0.10,0.10,0.10,0.08,0.09,0.09)
    mat[3,] <- c(0.20,0.14,0.12,0.11,0.09,0.09,0.08,0.07,0.07,NA,NA)
    mat[4,] <- c(0.18,0.13,0.11,0.10,0.09,0.09,0.09,0.07,NA,NA,NA)
    mat[5,] <- c(0.16,0.11,0.09,0.10,0.06,0.08,NA,NA,NA,NA,NA)
    mat[6,] <- c(0.16,0.12,0.11,0.08,0.08,NA,NA,NA,NA,NA,NA)
    mat[7,] <- c(0.16,0.10,NA,NA,NA,NA,NA,NA,NA,NA,NA)
    mat[8,] <- c(0.13,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
    mat[9,] <- c(0.13,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
  }else if(var_type == "W"){
    stop("Sampling uncertainty values are not implemented for W so far.")
			class_inf <- c(seq(0.01,0.04,0.01),seq(0.06,0.1,0.02),0.2)
  }else if(var_type == "Z"){
    #### Zh ####
    class_inf <- c(seq(10,50,5))
    mat <- matrix(NA,nrow=8,ncol=n_tres)
    if(freq == 9.4){
	mat[1,] <- c(0.12,0.09,0.08,0.07,0.07,0.06,0.06,0.05,0.06,0.06,0.05)
	mat[2,] <- c(0.11,0.08,0.07,0.06,0.06,0.06,0.05,0.05,0.05,0.05,0.05)
	mat[3,] <- c(0.09,0.06,0.05,0.05,0.05,0.05,0.04,0.04,0.04,0.04,0.04)
	mat[4,] <- c(0.08,0.06,0.05,0.05,0.04,0.04,0.04,0.04,0.05,0.05,0.04)
	mat[5,] <- c(0.09,0.07,0.06,0.06,0.06,0.06,0.06,0.06,0.05,0.05,0.05)
	mat[6,] <- c(0.14,0.11,0.09,0.07,0.07,0.07,0.06,0.05,0.04,0.03,0.03)
	mat[7,] <- c(0.17,0.10,0.09,0.08,0.07,0.06,0.05,0.04,  NA,  NA,  NA)
	mat[8,] <- c(0.15,0.13,0.06,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
    }else if(freq == 5.6){
	mat[1,] <- c(0.12,0.09,0.08,0.07,0.07,0.06,0.06,0.05,0.06,0.06,0.05)
	mat[2,] <- c(0.11,0.08,0.07,0.06,0.06,0.06,0.05,0.05,0.05,0.05,0.05)
	mat[3,] <- c(0.09,0.07,0.06,0.05,0.05,0.05,0.04,0.04,0.04,0.03,0.03)
	mat[4,] <- c(0.09,0.06,0.05,0.04,0.04,0.04,0.04,0.03,0.03,0.03,0.03)
	mat[5,] <- c(0.09,0.06,0.04,0.04,0.04,0.04,0.03,0.03,0.04,0.03,0.03)
	mat[6,] <- c(0.08,0.05,0.04,0.04,0.04,0.04,0.04,0.05,0.07,0.06,0.06)
	mat[7,] <- c(0.06,0.11,0.06,0.13,0.13,0.13,0.09,  NA,  NA,  NA,  NA)
	mat[8,] <- c(0.17,0.19,0.20,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
    }else if(freq == 2.8){
	mat[1,] <- c(0.12,0.09,0.08,0.07,0.07,0.06,0.06,0.05,0.06,0.06,0.06)
	mat[2,] <- c(0.11,0.08,0.07,0.06,0.06,0.06,0.05,0.05,0.05,0.05,0.05)
	mat[3,] <- c(0.09,0.07,0.06,0.05,0.05,0.05,0.04,0.04,0.04,0.04,0.04)
	mat[4,] <- c(0.09,0.06,0.05,0.04,0.04,0.04,0.04,0.04,0.03,0.03,0.03)
	mat[5,] <- c(0.09,0.06,0.05,0.04,0.04,0.04,0.03,0.04,0.03,0.03,0.03)
	mat[6,] <- c(0.09,0.06,0.05,0.04,0.04,0.04,0.03,0.04,0.03,0.02,0.02)
	mat[7,] <- c(0.09,0.08,0.07,0.06,0.05,0.04,  NA,  NA,  NA,  NA,  NA)
	mat[8,] <- c(0.11,0.04,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
    }else{stop("{Fct 'uncertainty_Parsiv'}: Radar frequency not available.")}
  }else if(var_type == "Zdr"){
    #### Zdr ####
    class_inf <- c(seq(0.1,0.5,0.1),0.75,1,1.5,2,3)
    mat <- matrix(NA,nrow=9,ncol=n_tres)
    if(freq == 9.4){
	mat[1,] <- c(0.16,0.12,0.11,0.10,0.10,0.10,0.08,0.08,0.08,0.08,0.08)
	mat[2,] <- c(0.22,0.18,0.16,0.15,0.14,0.13,0.12,0.11,0.11,0.11,0.10)
	mat[3,] <- c(0.28,0.22,0.20,0.19,0.18,0.17,0.17,0.16,0.14,0.13,0.13)
	mat[4,] <- c(0.31,0.26,0.23,0.23,0.21,0.21,0.17,0.17,0.15,0.14,0.14)
	mat[5,] <- c(0.36,0.31,0.27,0.25,0.25,0.24,0.21,0.21,0.22,0.19,0.21)
	mat[6,] <- c(0.41,0.37,0.35,0.32,0.31,0.30,0.30,0.31,0.31,0.32,0.30)
	mat[7,] <- c(0.47,0.43,0.40,0.39,0.40,0.41,0.38,0.35,0.32,0.35,0.32)
	mat[8,] <- c(0.50,0.46,0.44,0.43,0.41,0.39,0.36,0.34,0.30,0.28,0.24)
	mat[9,] <- c(0.49,0.44,0.41,0.39,0.38,0.32,0.31,0.35,0.23,0.23,0.16)
    }else if(freq == 5.6){
	mat[1,] <- c(0.16,0.12,0.11,0.10,0.10,0.10,0.08,0.08,0.08,0.08,0.08)
	mat[2,] <- c(0.22,0.18,0.17,0.15,0.15,0.14,0.12,0.11,0.11,0.11,0.10)
	mat[3,] <- c(0.28,0.22,0.20,0.19,0.18,0.17,0.17,0.15,0.14,0.13,0.14)
	mat[4,] <- c(0.30,0.26,0.23,0.23,0.21,0.20,0.17,0.16,0.16,0.13,0.12)
	mat[5,] <- c(0.35,0.29,0.25,0.23,0.21,0.20,0.18,0.18,0.15,0.14,0.15)
	mat[6,] <- c(0.39,0.33,0.27,0.24,0.24,0.25,0.21,0.17,0.19,0.19,0.19)
	mat[7,] <- c(0.44,0.37,0.33,0.32,0.28,0.27,0.25,0.25,0.21,0.21,0.23)
	mat[8,] <- c(0.50,0.42,0.41,0.41,0.42,0.37,0.37,0.39,  NA,  NA,  NA)
	mat[9,] <- c(0.53,0.50,0.47,0.49,0.51,0.47,0.53,0.56,  NA,  NA,  NA)
    }else if(freq == 2.8){
	mat[1,] <- c(0.16,0.12,0.11,0.10,0.10,0.10,0.08,0.08,0.08,0.08,0.08)
	mat[2,] <- c(0.22,0.18,0.17,0.15,0.15,0.14,0.12,0.11,0.11,0.11,0.10)
	mat[3,] <- c(0.28,0.22,0.20,0.19,0.19,0.18,0.17,0.16,0.14,0.13,0.14)
	mat[4,] <- c(0.30,0.26,0.23,0.23,0.21,0.21,0.17,0.16,0.16,0.14,0.13)
	mat[5,] <- c(0.35,0.29,0.26,0.23,0.22,0.21,0.19,0.18,0.16,0.14,0.15)
	mat[6,] <- c(0.39,0.33,0.27,0.25,0.24,0.26,0.23,0.20,0.21,0.22,0.20)
	mat[7,] <- c(0.43,0.38,0.35,0.33,0.31,0.29,0.27,0.28,0.25,0.21,0.22)
	mat[8,] <- c(0.49,0.43,0.41,0.40,0.40,0.38,0.38,0.38,0.36,  NA,  NA)
	mat[9,] <- c(0.50,0.47,0.42,0.43,0.43,0.37,0.26,0.29,0.17,0.25,  NA)
    }else{stop("{Fct 'uncertainty_Parsiv'}: Radar frequency not available.")}
  }else if(var_type == "ah"){
    #### ah ####
    if(freq == 9.4){
	class_inf <- c(0,0.05,0.1,0.2,0.3,0.4,0.5,0.6)		# X-band
	mat <- matrix(NA,nrow=length(class_inf)-1,ncol=n_tres)
	mat[1,] <- c(0.31,0.23,0.19,0.18,0.17,0.16,0.15,0.15,0.14,0.14,0.14)
	mat[2,] <- c(0.54,0.32,0.26,0.22,0.20,0.20,0.16,0.14,0.11,0.10,0.08)
	mat[3,] <- c(0.46,0.31,0.22,0.21,0.18,0.15,0.12,0.12,  NA,  NA,  NA)
	mat[4,] <- c(0.42,0.24,0.22,0.20,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[5,] <- c(0.33,0.24,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[6,] <- c(0.32,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[7,] <- c(0.31,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
    }else if(freq == 5.6){
 	class_inf <- c(0,0.01,seq(0.02,0.1,0.01))		# C-band
	mat <- matrix(NA,nrow=length(class_inf)-1,ncol=n_tres)
	mat[1,] <- c(0.26,0.18,0.15,0.14,0.13,0.13,0.12,0.11,0.11,0.12,0.11)
	mat[2,] <- c(0.36,0.29,0.24,0.30,0.28,0.29,0.26,0.25,0.17,0.20,  NA)
	mat[3,] <- c(0.71,0.39,0.66,0.52,0.47,0.44,  NA,  NA,  NA,  NA,  NA)
	mat[4,] <- c(0.64,0.88,0.48,0.39,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[5,] <- c(0.45,0.90,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[6,] <- c(0.52,0.76,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[7,] <- c(0.37,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[8,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[9,] <- c(1.06,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[10,] <- c(0.98,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
   }else if(freq == 2.8){
	class_inf <- c(0,0.01,seq(0.02,0.1,0.01))		# S-band
	mat <- matrix(NA,nrow=length(class_inf)-1,ncol=n_tres)
	mat[1,] <- c(0.22,0.16,0.13,0.12,0.11,0.11,0.10,0.10,0.09,0.09,0.09)
	mat[2,] <- c(0.25,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[3,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[4,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[5,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[6,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[7,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[8,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[9,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[10,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
   }else{stop("{Fct 'uncertainty_Parsiv'}: Radar frequency not available.")}
  }else if(var_type == "kdp"){
    #### kdp ####
    if(freq == 9.4){
	class_inf <- c(0,0.1,seq(0.2,1,0.1))			# X-band
	mat <- matrix(NA,nrow=length(class_inf)-1,ncol=n_tres)
	mat[1,] <- c(0.35,0.25,0.21,0.19,0.18,0.17,0.16,0.15,0.14,0.14,0.14)
	mat[2,] <- c(0.36,0.24,0.20,0.19,0.18,0.17,0.16,0.17,0.13,0.12,0.12)
	mat[3,] <- c(0.34,0.25,0.20,0.17,0.16,0.15,0.14,0.12,0.11,  NA,  NA)
	mat[4,] <- c(0.32,0.21,0.19,0.15,0.14,0.15,0.10,  NA,  NA,  NA,  NA)
	mat[5,] <- c(0.33,0.24,0.18,0.16,0.16,0.14,  NA,  NA,  NA,  NA,  NA)
	mat[6,] <- c(0.33,0.27,0.15,0.17,0.17,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[7,] <- c(0.33,0.23,0.18,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[8,] <- c(0.33,0.25,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[9,] <- c(0.50,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[10,] <- c(0.36,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
    }else if(freq == 5.6){
	class_inf <- c(0,0.1,seq(0.2,1,0.1))			# C-band
	mat <- matrix(NA,nrow=length(class_inf)-1,ncol=n_tres)
	mat[1,] <- c(0.35,0.25,0.21,0.19,0.18,0.17,0.16,0.15,0.15,0.15,0.14)
	mat[2,] <- c(0.40,0.29,0.22,0.19,0.19,0.18,0.15,0.14,0.10,0.11,0.09)
	mat[3,] <- c(0.45,0.26,0.20,0.17,0.16,0.16,0.12,  NA,  NA,  NA,  NA)
	mat[4,] <- c(0.45,0.30,0.20,0.20,0.18,0.13,  NA,  NA,  NA,  NA,  NA)
	mat[5,] <- c(0.40,0.21,0.21,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[6,] <- c(0.42,0.20,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[7,] <- c(0.38,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[8,] <- c(0.32,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[9,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[10,] <- c(0.31,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
    }else if(freq == 2.8){
	class_inf <- c(0,0.05,0.1,seq(0.15,0.5,0.05))		# S-band
	mat <- matrix(NA,nrow=length(class_inf)-1,ncol=n_tres)
	mat[1,] <- c(0.34,0.25,0.21,0.19,0.18,0.17,0.16,0.15,0.15,0.14,0.14)
	mat[2,] <- c(0.38,0.26,0.21,0.18,0.17,0.17,0.15,0.15,0.11,0.09,0.08)
	mat[3,] <- c(0.37,0.26,0.22,0.20,0.21,0.17,0.10,  NA,  NA,  NA,  NA)
	mat[4,] <- c(0.42,0.30,0.19,0.20,0.16,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[5,] <- c(0.39,0.30,0.22,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[6,] <- c(0.43,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[7,] <- c(0.73,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[8,] <- c(0.33,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[9,] <- c(0.37,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
	mat[10,] <- c(  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA,  NA)
    }else{stop("{Fct 'uncertainty_Parsiv'}: Radar frequency not available.")}
  }else{
    stop("No uncertainty values for the specified variable !")
  }
  
  if(fill.na == TRUE){		## fill NA values with the last available uncertainty value (only for Nt, R and kdp)
      if(sum(var_type == c("Nt","R","Z","kdp"))==1){
	  isna <- is.na(mat)
	  for(c in 1:n_tres){
	      pos <- max(which(isna[,c] == FALSE))
	      mat[isna[,c],c] <- mat[pos,c]
	  }
      }
  }

  if(sum(dt == list_tres) != 0){
    vect_mat <- mat[,pos_dt]
  }else{	## interpolate uncertainty values using least squared if specific time res.
    pos_1 <- max(which(dt > list_tres)) 
    pos_2 <- min(which(dt < list_tres))
    lambda1 <- 1/(dt - list_tres[pos_1])
    lambda2 <- 1/(list_tres[pos_2] - dt)
    int <- ((lambda1*mat[,pos_1]) + (lambda2*mat[,pos_2]))/(lambda1 + lambda2)
    vect_mat <- round(int,digits=2)
  }
  class_limits <- matrix(data=sort(c(class_inf,class_inf[2:(length(class_inf)-1)])),nrow=length(class_inf)-1,ncol=2,byrow=T)
  class_spread <- (class_limits[,2]-class_limits[,1])
  class_avg <- c(((class_spread/2) + class_limits[,1]),NA)
  number_class <- length(class_avg)
#   legend_classes <- c(apply(round(class_limits,digits=3),1,paste,collapse=","),paste(">",class_limits[number_class-1,2],sep=""))
  logical_class <- matrix(NA,nrow=n_meas,ncol=number_class)

  for (class_lim in 1:number_class){
    if (class_lim != number_class){
      cond_inf <- vect >= class_limits[class_lim,1]
      cond_sup <- vect <  class_limits[class_lim,2]
      logical_class[,class_lim] <- as.logical(cond_inf * cond_sup)
    }else{logical_class[,class_lim] <- as.logical(vect >= class_limits[(class_lim-1),2])}

    temp <- which(logical_class[,class_lim] == TRUE)
    uncer[temp,1] <- vect_mat[class_lim]
  }
  if(fill.na == TRUE){
      uncer[temp,1] <- vect_mat[class_lim-1]
      cond_inf <- vect < class_limits[1,1]
      temp <- which(cond_inf == TRUE)
      uncer[temp,1] <- vect_mat[1]
  }else{uncer[temp,1] <- NA}
  n_meas_classes <- colSums (logical_class, na.rm = TRUE)		# number of measurements in each class of rain rate.

  uncer[,2] <- uncer[,1] * vect

  return(uncer)
}