Rutils/maybe-not-useful/angle.rectangle.r

In manfredo89/ED2io: Manages Input/Output for Ecosystem Demography 2 Model

#==========================================================================================#
#==========================================================================================#
#     Function that finds the orientation angle.  If outmean is TRUE, then the mean        #
# orientation angle is provided instead.                                                   #
#------------------------------------------------------------------------------------------#
angle.rectangle <<- function(x,y,degree=TRUE,outmean=TRUE){
   #----- Check that x and y have the same length. ----------------------------------------#
   nx = length(x)
   ny = length(y)
   if (nx != ny){
      cat("--------------------------------------------------------------","\n",sep="")
      cat("  length(x) = ",nx,"\n",sep="")
      cat("  length(y) = ",ny,"\n",sep="")
      cat("--------------------------------------------------------------","\n",sep="")
      stop("x and y must have the same length")
   }#end if
   #---------------------------------------------------------------------------------------#


   #---------------------------------------------------------------------------------------#
   #      Find the difference in x and y.                                                  #
   #---------------------------------------------------------------------------------------#
   dx       = diff(x)
   dy       = diff(y)
   angle    = atan2(dy,dx) / pio180
   angle.12 = angle %% 180.
   angle.41 = ( ( angle + 90 ) %% 180 ) - 90
   dist     = sqrt(dx*dx+dy*dy)
   #---------------------------------------------------------------------------------------#




   #---------------------------------------------------------------------------------------#
   #      Check whether to return the mean angle or the angle for each segment.            #
   #---------------------------------------------------------------------------------------#
   if (outmean){
      #----- Find the mean,median, and standard deviation of the angles. ------------------#
      a12.median = weighted.quantile(x=angle.12,w=dist,na.rm=TRUE)$q
      a12.sdev   = weighted.sd      (x=angle.12,w=dist,na.rm=TRUE)
      a41.median = weighted.quantile(x=angle.41,w=dist,na.rm=TRUE)$q
      a41.sdev   = weighted.sd      (x=angle.41,w=dist,na.rm=TRUE)
      #------------------------------------------------------------------------------------#


      #------------------------------------------------------------------------------------#
      #     Decide which block to use based on the standard deviation.  The lowest is less #
      # likely to be near the endpoints.                                                   #
      #------------------------------------------------------------------------------------#
      if (a41.sdev < a12.sdev){
         ans = a41.median + 180
      }else{
         ans = a12.median
      }#end if
      #------------------------------------------------------------------------------------#
   }else{
      ans = angle.12
   }#end if
   #---------------------------------------------------------------------------------------#


   #----- Final answer. -------------------------------------------------------------------#
   if (! degree) ans = ans * pio180
   return(ans)
   #---------------------------------------------------------------------------------------#
}#end function angle.rectangle
#==========================================================================================#
#==========================================================================================#