############################################################ #Created by: Stefan Metzger #Last change by: Stefan Metzger (2011-11-06) ############################################################ #WHAT FOR #calculates one footprint weight matrix at a time from analytical models #for the crosswind integrated footprint Kormann & Meixner, 2001 (KM01) or Kljun, 2004 (K04) model can be chosen #for the crosswind distribution the same Gaussian function is used in both cases (Pasquill, 1974; Kormann and Meixner, 2001) #the matrix extend is precalculated, and cut off where the footprint function reaches 1 % of its peak (Neftel, 2008) #matrix extends are uneven, with the aircraft / tower at its center #time consumption ~1/3 s per footprint for KM01 #K04 is 3-4 times faster #INPUT VARIABLES #angle #wind direction [°] to rotate the inertial footprint matrix; used for both models #Csize #cell size [m] of result grid; used for both models #hor #horizontal wind speed [m/s]; used in KM01, and to determine z0 for K04 #sigV #crosswind fluctuations [m/s]; used for both models to calculate the crosswind dispersion #sigW #vertical wind fluctuations [m/s]; used for K04 to characterize the vertical transport #ustar #friction velocity [m/s]; used a) for both models to characterize the shear stress, and b) to determine z0 for K04 #Lvirt #Obukhov length from local friction velocity and buoyancy flux [m]; used a) for KM01 to characterize the vertical transport, and b) for the universal functions to determine z0 for K04 #zmeas #height of measurement - displacement [m] ############################################################ #KM01 MODEL ############################################################ footKM01 <- function( angle, #wind direction [°] to rotate the inertial footprint matrix; used for both models Csize, #cell size [m] of result grid; used for both models hor, #horizontal wind speed [m/s]; used in KM01, and to determine z0 for K04 sigV, #crosswind fluctuations [m/s]; used for both models to calculate the crosswind dispersion ustar, #friction velocity [m/s]; used a) for both models to characterize the shear stress, and b) to determine z0 zmeas, #height of measurement - displacement [m] Lvirt, # cutoff=0.01 #cellwise contribution relative to peak at which footprint will be cut off ){ #----------------------------------------------------------- #PREPARATION OF INPUT PARAMETERS #calculation of local stability zeta [-] zeta <- zmeas / Lvirt #m exponents of the wind velocity power law if (zeta >= 0) { phi_m <- 1 + 5*zeta } else { phi_m <- (1-16*zeta)^(-0.25) } m <- ustar * phi_m / 0.4 / hor #n exponent of eddy diffusivity power law if (zeta >= 0) { n <- 1 /(1+5*zeta) } else { n <- (1-24*zeta)/(1-16*zeta) } #r shape factor r <- 2 + m - n #U constant in power law profile in the wind velocity U <- hor / (zmeas^m) #K vertical profile of the eddy diffusivity if (zeta >= 0) { phi_c <- 1 + 5*zeta } else { phi_c <- (1-16*zeta)^(-0.5) } K <- 0.4*ustar * zmeas /phi_c #kapa constant of the power law profile of the eddy diffusivity kapa <- K / zmeas^n #----------------------------------------------------------- #SIZE ESTIMATION OF FOOTPRINT MATRIX, SMALLEST CONTRIBUTION cutoff % OF PEAK CONTRIBUTION (NEFTEL, 2008) #alongwind (crosswind integrated) density distribution #position of maximum along x, 290.5334 m Xmax <- function(mue, xi) xi / (1+mue) xmax <- Xmax(mue=(1+m)/r, xi=(U*zmeas^r)/((r^2)*kapa)) #maximum value, 0.001369324 Fmax <- function(mue, xi) ( exp(-1-mue) * (1+mue)^(1+mue) ) / (gamma(mue) * xi) fmax <- Fmax(mue=(1+m)/r, xi=(U*zmeas^r)/((r^2)*kapa)) #crosswind integrated flux footprint FFPalong <- function(x, mue, xi) { f <- (1 / gamma(mue)) * ((xi^mue)/x^(1+mue)) * exp(-xi/x) return(f) } #length until contribution falls below 1% of fmax whri <- xmax; whro <- fmax #start from distribution peak while(whro > fmax * cutoff) { whri <- whri + Csize #use step width of landuse matrix whro <- FFPalong(whri, mue=(1+m)/r, xi=(U*zmeas^r)/((r^2)*kapa)) #calculate } whrx <- ceiling(whri / Csize) #cell length necessay in X direction #crosswind density distribution #crosswind distribution of footprint FFPcross <- function(x, y, sigmav, u) { sigma <- sigmav * x / u Dy <- (1 / (sqrt(2 * pi) * sigma)) * exp((-y^2) / (2 * (sigma^2))) return(Dy) } #maximum of crosswind density distribution ymax <- 0 fmax <- FFPcross(whrx * Csize, ymax, sigmav=sigV, u=hor) #length until contribution falls below 1 % fmax whri <- ymax; whro <- fmax #start from distribution peak while(whro > fmax * cutoff) { whri <- whri + Csize #use step width of landuse matrix whro <- FFPcross(whrx * Csize, whri, sigmav=sigV, u=hor) #calculate } whry <- ceiling(whri / Csize) #cell length necessay in Y direction #----------------------------------------------------------- #CELL ALLOCATION AND INTEGRATION #place aircraft in center cell of first row Xbou <- c(0, 1:whrx - 0.5)* Csize #alongwind integration boundaries with aicraft centered in 0 Ybou <- c(0, 1:whry - 0.5)* Csize #crosswind integration boundaries with aicraft centered in 0 Xcen <- sapply(1:(length(Xbou)-1), function(x) mean(Xbou[x:(x+1)])) #alongwind cell center coordinates # Ycen <- sapply(1:(length(Ybou)-1), function(y) mean(Xbou[y:(y+1)])) #crosswind cell center coordinates #integration of alongwind footprint PHIalong <- sapply(1:(length(Xbou)-1), function(xwhr) integrate(FFPalong, Xbou[xwhr], Xbou[xwhr+1], mue=(1+m)/r, xi=(U*zmeas^r)/((r^2)*kapa))$value) sumPHIalong <- sum(PHIalong) PHIalong <- PHIalong / sumPHIalong #normalisation to 1 #integration of crosswind footprint #function for crosswind dispersion FFPcrossY <- function(y, sigma) { Dy <- (1 / (sqrt(2 * pi) * sigma)) * exp((-y^2) / (2 * (sigma^2))) return(Dy) } #alongwind distance dependence of crosswind dispersion FFPcrossXY <- function(x, y, sigmav, u) { sigma <- sigmav * x / u PHIcross <- sapply(1:(length(y)-1), function(ywhr) integrate(FFPcrossY, y[ywhr], y[ywhr+1], sigma)$value) PHIcross <- PHIcross / (2 * sum(PHIcross)) #normalisation to 0.5 return(PHIcross) } #integration PHIcross <- t(sapply(1:length(Xcen), function(xwhr) FFPcrossXY(Xcen[xwhr], Ybou, sigmav=sigV, u=hor))) #----------------------------------------------------------- #COMBINE ALONG- AND CROSSWIND DENSITY DISTRIBUTIONS AND ROTATE INTO MEAN WIND #combine crosswind contributions on alongwind axis; will always yield uneven column number PHIcross <- cbind(fliplr(PHIcross[,2:ncol(PHIcross)]), 2*PHIcross[,1], PHIcross[,2:ncol(PHIcross)]) #YcenLR <- c(-rev(Ycen[2:length(Ycen)]), 0, Ycen[2:length(Ycen)]) #sapply(1:nrow(PHIcross), function(x) sum(PHIcross[x,])) #str(PHIcross) #combination of along- and cross component; along wind from up to down, crosswind from left to right PHI <- t(sapply(1:nrow(PHIcross), function(x) PHIalong[x] * PHIcross[x,])) #plot(PHIalong) #plot(PHIcross[5,]) #contour(rot90(PHI,3), levels=c(0.0001, 0.001, 0.01)) #center aicraft alongwind location in plot, pad with zeroes PHIc<-rbind(matrix(nrow=(nrow(PHI)-1), ncol=ncol(PHI), 0), PHI) #XcenUD <- c(-rev(Xcen[2:length(Xcen)]), 0, Xcen[2:length(Xcen)]) #contour(YcenLR, XcenUD, rot90(PHIc,1), levels=NIVo, labels=NIVi, col=colorRampPalette(c("black", "red"))(length(NIVo)), asp=1) #pad with zeroes if not rectangular matrix if(nrow(PHIc) == ncol(PHIc)) { PHIcp <- PHIc } else { if(nrow(PHIc) > ncol(PHIc)) { PHIcp <- cbind( matrix(nrow=nrow(PHIc), ncol=(nrow(PHIc)-ncol(PHIc)) / 2, 0), PHIc, matrix(nrow=nrow(PHIc), ncol=(nrow(PHIc)-ncol(PHIc)) / 2, 0) ) } if(nrow(PHIc) < ncol(PHIc)) { PHIcp <- rbind( matrix(ncol=ncol(PHIc), nrow=(ncol(PHIc)-nrow(PHIc)) / 2, 0), PHIc, matrix(ncol=ncol(PHIc), nrow=(ncol(PHIc)-nrow(PHIc)) / 2, 0) ) } } #YcenLRp <- XcenUD #contour(YcenLRp, XcenUD, rot90(PHIcp,1), levels=NIVo, labels=NIVi, col=colorRampPalette(c("black", "red"))(length(NIVo)), asp=1) #rotate image clockwise (align footprint in mean wind) PHIcpr <- rotate(PHIcp, 180-angle)@.Data #YcenLRpr <- ((1:ncol(PHIcpr)) - ceiling(ncol(PHIcpr) / 2)) * Csize #new X #XcenUDpr <- ((1:nrow(PHIcpr)) - ceiling(nrow(PHIcpr) / 2)) * Csize #new Y #contour(YcenLRpr, XcenUDpr, rot90(PHIcpr,3), levels=NIVo, labels=NIVi, col=colorRampPalette(c("black", "red"))(length(NIVo)), asp=1) #----------------------------------------------------------- #RETURN RESULTS #create result list export<-list(PHIcpr, sumPHIalong) names(export) <- c("PHI", "sum") #return list return(export) }