This repository has been archived by the owner on Jan 6, 2022. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 29
/
Copy pathprogramming.examples.Rmd
158 lines (130 loc) · 4.48 KB
/
programming.examples.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
---
title: "Programming-Notes"
author: "Naomi Tague"
date: "January 28, 2016"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Code from inclass examples
Working with data types and functions
```{r powergen}
# create a function that computes power generated from a reservoir
power_gen = function(height, flow, rho=1000, g=9.8, Keff=0.8) {
result = rho * height * flow * g * Keff
return(result)
}
# examples of function use
power_gen(20,1)
power_gen(height=20, flow=1)
power.estimate = power_gen(height=20,flow=1)
power.estimate
```
## Scoping
```{r scoping, echo=FALSE}
# this will fail
#power.estimate = power_gen(Keff=0.6)
power.estimate = power_gen(height=20, flow=0.6, Keff=0.6)
# this will fail
#Keff
```
## Nicely documented function
```{r speed }
#' Power Required by Speed
#'
#' This function determines the power required to keep a vehicle moving
#' a given speed
#' @param cdrag coefficient due to drag default=0.3
#' @param crolling coefficient due to rolling/friction default=0.015
#' @param v vehicle speed (m/2)
#' @param m vehicle mass (kg)
#' @param A area of front of vehicle (m2)
#' @param g acceleration due to gravity (m/s) default=9.8
#' @param pair (kg/m3) default =1.2
#' @return power (W)
power = function(cdrag=0.3, crolling=0.015,pair=1.2,g=9.8,V,m,A) {
P = crolling*m*g*V + 1/2*A*pair*cdrag*V**3
return(P)
}
v=seq(from=0, to=100, by=10)
plot(v, power(V=0.447*v, m=31752, A=25))
lines(v, power(V=0.447*v, m=61752, A=25))
```
## Data types in R (that you can exchange with functions)
```{r datafunctions }
mth.names=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct"
,"Nov","Dec")
reservoir.operation = data.frame(month=mth.names)
reservoir.operation$height = seq(from=32, to=10, by=-2)
reservoir.operation$flowrate = rnorm(n=12, mean=3, sd=0.25)
barplot( power_gen(height=reservoir.operation$height, flow=reservoir.operation$flow), names=reservoir.operation$month)
```
## Using other functions within a function
```{r morefunctions }
#'
#' This function computes total power generation from a reservoir given its height and flow rate into turbines and number of days (and secs) within those days that the turbines are in operation
#' @param rho Density of water (kg/m3) Default is 1000
#' @param g Acceleration due to gravity (m/sec2) Default is 9.8
#' @param Keff Turbine Efficiency (0-1) Default is 0.8
#' @param height height of water in reservoir (m)
#' @param flow flow rate (m3/sec)
#' @param number of days
#' @param secs in days Default is 86400
#' @author Naomi
#' @examples power_gen(20, 1, 10)
#' @return Power generation (MW)
power_gen_total = function(height, flow, days, secs=86400, rho=1000, g=9.8, Keff=0.8) {
result = rho * height * flow * g * Keff
result = result * days * secs
total = sum(result)/1e6
return(total)
}
power_gen_total(reservoir.operation$height, reservoir.operation$flowrate, days=30)
```
## Using lists to return multiple data types from a function
```{r listreturns }
#' computes profit from price for forest plot and Mg/C in that plot
#' @param price ($)
#' @param carbon (MgC)
#' @return list with mean, min, and max prices
compute_carbonvalue = function(price, carbon) {
cost.per.carbon = price/carbon
a = mean(cost.per.carbon)
b = max(cost.per.carbon)
c = min(cost.per.carbon)
result = list(avg=a, min=c, max=b)
return(result)
}
# example application, using some data that we 'generate'
obs=data.frame(prices=c(23,44,60,4,2,33,59), forestc=c(59,88,100,10,8,79,300) )
obs
forestvalue = compute_carbonvalue(obs$prices, obs$forestc)
forestvalue
```
## Function that uses factors (and has some wierd comments but ignore that for now
```{r diversity}
#' Describe diversity based on a list of species
#'
#' Compute a species diversity index
#' @param species list of species (names, or code)
#' @return list with the following items
#' \describe{
#' \item{num}{ Number of distinct species}
#' \item{simpson}{Value of simpson diversity index}
#' \item{dominant}{Name of the most frequently occuring species}
#' }
#' @examples
#' compute_diversity(c("butterfly","butterfly","mosquito","butterfly","ladybug","ladybug")))
#' @references
#' http://www.tiem.utk.edu/~gross/bioed/bealsmodules/simpsonDI.html
compute_diversity = function(species) {
species = as.factor(species)
tmp = (summary(species)/sum(summary(species))) ** 2
diversity = sum(tmp)
nspecies = length(summary(species))
tmp = which.max(summary(species))
dominant = names(summary(species)[tmp])
return(list(num=nspecies, simpson=diversity, dominant=dominant))
}
```