xyz.plot <- function (data, z, x, y, 
                        residual.plot = FALSE,
                        lower.limit,
                        upper.limit,
                        origin,
                        maxsize = 0.075, 
                        seg.col = "red", 
                        seg.lty = par("lty"),
                        seg.lwd = par("lwd"), 
                        seg.azm = 0,
                        pch.pnt = 16,
                        plot.points = FALSE,
                        equal.scale = TRUE, 
                        messages.screen = FALSE, ... ) {

# purpose:          function plots the x,y coordinates of a point and
#                   draws line segments through the points whereby the
#                   length of the segments is linearly related to z
#                   
#                   
# arguments:
#
# data              dataframe with x, y, and z-components or 
#                   list of class "geodata" from which missing 
#                   x, y, z are taken                   
# z                 value of target variable
# x                 x-coordinate of sampling location
# y                 y-coordinate of sampling location
# residual.plot     TRUE: residual type plot, 
#                       sign(z-lower.limit) 45 degrees slanted segments, 
#                   FALSE: normal type plot, vertical or horizontal segments
# lower.limit       value for lower.limit of target variable, 
#                   default value: lower.limit=mean(z),
#                   there are three possibilities for drawing segments
#                   1)  lower.limit < min(z): normal type plot, 
#                       segment length proportional to (z - lower.limit)
#                   2)  lower.limit in range(z) and residual.plot = TRUE: residual type plot;
#                       segment length proportional to abs(z - lower.limit)
#                   3)  else: normal type plot, 
#                       segment length proportional to (z - min(z))
# upper.limit       value for upper.limit of target varible
#                   default value: upper.limit=max(z)
#                   upper.limit > max(z) then the segment length is given by
#                   (z - lower.limit) / (upper.limit - lower.limit)
# origin            the same as lower.limit (compatibility)
# maxsize           maximum length of segment in y-units 
# equal.scale       TRUE: calls eqscplot() of library MASS
#                   FALSE: calls basic plot()
# seg.col           color attribute for segments
# seg.lty           line type attribute for segments
# seg.lwd           line width attribute for segments
# seg.azm           azimuth angle of segments in degrees (angle clockwise from north)
# pch.pnt           plotting character for plot()
# plot.points       TRUE: sampling locations are shown
#                   FALSE: sampling locations are not shown 
# messages.screen   write diagnostics to stdout
# ...               further arguments passed to plot(), eqscplot(), points()
#               
# author:           A. Papritz
# date:             2003-04-01


# setup x, y, z vectors

    if ( !missing(data) ) {
        if( missing(x) ) {
            if( is.data.frame(data) ) {
                if( any( names(data) == 'x' ) ) x <- data$x
                else stop('dataframe without "x"-component ')
            }
            else {
                if( is.list(data) && class(data) == "geodata")
                    x <- data$coords[,1]
                else stop('object given for argument "data" not supported')
            }
        }
        if( missing(y) ) {
            if( is.data.frame(data) ){
                if( any( names(data) == 'y' ) ) y <- data$y
                else stop('dataframe without "y"-component ')
            }
            else {
                if( is.list(data) && class(data) == "geodata")
                    y <- data$coords[,2]
                else stop('object given for argument "data" not supported')
            }
        }
        if( missing(z) ) {
            if( is.data.frame(data) ){
                if( any( names(data) == 'z' ) ) z <- data$z
                else stop('dataframe without "z"-component ')
            }
            else {
                if( is.list(data) && class(data) == "geodata")
                    z <- data$data
                else stop('object given for argument "data" not supported')
            }
        }
    }
    else {
        if ( missing(x) ) stop('"x" is a required argument')
        if ( missing(y) ) stop('"y" is a required argument')
        if ( missing(z) ) stop('"z" is a required argument')
    }
    
    
# check mode of important arguments
    
    if( !is.numeric( x ) ) stop( "x-coordinate must be of mode numeric" )
    if( !is.numeric( y ) ) stop( "y-coordinate must be of mode numeric" )
    if( !is.numeric( z ) ) stop( "z-value must be of mode numeric" )
    if( !( is.character( seg.col ) | is.numeric( seg.col ) ) 
        ) stop("seq.col must be either of mode character or numeric")
    if( !is.numeric( seg.azm )) stop( "seg.azm must be of mode numeric" )
        

# omit NAs in x, y or z

    n <- min( length(x), length(y), length(z) )
    
    if( length(seg.col) == 1 ) (seg.col <- rep(seg.col, n) )
    if( length(seg.azm) == 1 ) (seg.azm <- rep(seg.azm, n) )
    
    df <- data.frame(x = x[1:n],
        y = y[1:n], 
        z = z[1:n], 
        seg.col = I(seg.col[1:n]),
        seg.azm = seg.azm[1:n]
    )
    df <- na.omit(df)
    x <- df$x; y <- df$y; z <- df$z; 
    seg.col <- df$seg.col; seg.azm <- df$seg.azm
    
    if( messages.screen )
        cat( '  number of omitted cases with NAs: ', n-length(x), '\n' )
    rm(df)
    
 
# set default value for lower.limit

    if ( missing(lower.limit) ) {
        if (missing(origin))
            lower.limit <- mean(z)
        else lower.limit <- origin
    }


# plot sample locations

    if( plot.points ) 
        pt <- "p"
    else 
        pt <- "n"
    
    if( equal.scale )
        t.asp <- 1
    else
        t.asp <- NA
        
    plot(x, y, pch=pch.pnt, type = pt, asp = t.asp, ... )
    

# compute length of line segments

    wdim <- par('usr'); wx <- diff(wdim[1:2]); wy <- diff(wdim[3:4])
    if( messages.screen )
        cat('  extremes of user coordinates in plotting window ', 
            format(wdim), '\n')
        
    maxsize <- wy * maxsize * 0.5
    
    rz <- range(z)
    
    if( messages.screen )
        cat( '  range of z: ', rz, '\n' )
        
    if ( missing(upper.limit) )
        upper.limit <- rz[2]
    else {
        if (upper.limit <= rz[2])
            upper.limit <- rz[2]
    }

    if ( lower.limit < rz[1] ) {
        dy <- (z - lower.limit) / (upper.limit - lower.limit)
        if( messages.screen ) {
            cat( '  value used for lower.limit: ', lower.limit, '\n' )
            cat( '  length of segments proportional to ( z - lower.limit )\n' )
        }
    }
    else {
        if ( lower.limit < rz[2] && residual.plot ) {
            dy <- (z - lower.limit) / max( abs( c( rz[1]-lower.limit, rz[2]-lower.limit ) ) )
            if( messages.screen ) {
                cat( '  value used for lower.limit: ', lower.limit, '\n' )
                cat( '  length of segments proportional to ( z - lower.limit )\n' )
            }
        }
        else {
            dy <- (z-rz[1]) / (upper.limit - rz[1])
            if( messages.screen )
                cat( '  length of segments proportional to ( z - min(z) )\n' )
        }
    }
    dy <- dy * maxsize

# draw segments

    oldpar <- par(xpd=TRUE)
    pwin <- par('pin')
    
#     print(dy)
#     print(seg.col)
#     print(seg.azm)

    if ( all( dy >= 0 ) ) {
        dx <- sin(seg.azm/180*pi) * dy
        dy <- cos(seg.azm/180*pi) * dy
        
        segments( x-dx, y-dy, x+dx, y+dy, col=seg.col, lty=seg.lty, lwd=seg.lwd )
    }
    else 
        segments( x-abs(dy/wy*pwin[2]*wx/pwin[1]), y-dy, 
                       x+abs(dy/wy*pwin[2]*wx/pwin[1]), y+dy, 
                       col=seg.col, lty=seg.lty, lwd=seg.lwd )

    par(oldpar)

}



azimuth.angle <- function(x, y) {
    
# purpose:          function computes azimuth angle after centering
#                   on (center.x, center.y)
#                         
# arguments:
#
# x                 x-coordinates
# y                 y-coordinates
#               
# author:           A. Papritz
# date:             2001-06-05

    x <- x[!is.na(x)]
    y <- y[!is.na(y)]
    
    if(length(x) != length(y)) stop('lengths of "x" and "y" differ')
    
    angle <- -Arg( complex( real=x, imag=y ))
    angle[angle<0] <- 2*pi + angle[angle<0]
    angle <- angle + pi/2
    angle[angle>2*pi] <- angle[angle>2*pi] - 2*pi
    angle
}
