function(h_r_critical,SafetyMargin,fileinput,fileoutput){ # h_r_critical is a value you supply when you run the function
# safetyMargin <- 3
# reads from GopherusTb.csv into the data structure GopherusTb
GopherusTb <- read.csv(fileinput)
# MyExtinctionFn()
# These are local variable constants
# These are the slope and intercept relating h_r = intercept + slope*(Tmax-Tb)
# slope <- 0.535
# intercept <- 3.609
# Tb <- 35.025 # Gopherus Tb measured by Zimmermann et al. 1994
# file opens a file for writing "w+", and it asigns zz to point to it outside of [R] workshop
zz <- file(fileoutput,"w+") # kills the output file
# create header files for the file I will open in xls
cat(file=zz,"Lat",',',"Lon",',',"sumBefore",',',"sumAfter",',',"extb4",',',"extafter",'\n')
# the resolution is 10 arc minutes
resolution <- 10/60
# how many records are there in my GopherusTb data structure
print(length(GopherusTb$Lat)) # how many records were read in
# foreach record in GopherusTb execute the following instructions
for(i in 1: length(GopherusTb$Lat)){ # DO SOME LOOP, FOR
# only takes records above the equator! and west of Greenwich QUALITY CONTROL
#basic calculations for row lookups and column lookups
#coercing the floating point arithmetic into integer with as.integer()
r <- as.integer( ( 90 - GopherusTb$Lat[i])/resolution)
c <- as.integer( ( GopherusTb$Lon[i] + 180 ) / resolution )
# compute first tmax in May June and July which are long integers and divide by 30 to get degrees C
before <- (tmax_1975_11.10m[r,c]+tmax_1975_12.10m[r,c])/20
beforeMin <- (tmin_1975_11.10m[r,c]+tmin_1975_12.10m[r,c])/20
# do it again for 2080
after <- (tmax_2080_11.10m[r,c]+tmax_2080_12.10m[r,c])/20
afterMin <- (tmin_2080_11.10m[r,c]+tmin_2080_12.10m[r,c])/20
CTMAX <- GopherusTb$CTMAX[i]
sumBefore <- 0
sumAfter <- 0
for(h in 1:24){
for(i in 1:100){
if(CTMAX - SafetyMargin < (after-afterMin)/2*sin(pi/12*(h+i/100)-3*pi/4)+(after+afterMin)/2){
sumAfter <- sumAfter + 0.01
}
if(CTMAX - SafetyMargin < (before-beforeMin)/2*sin(pi/12*(h+i/100)-3*pi/4)+(before+beforeMin)/2){
sumBefore <- sumBefore + 0.01
}
}
}
# why divide by 30? because .bil files are long integers but they need to be converted to decimal
# results, read the readme in worldclim.org
if(sumBefore > h_r_critical){
extinct_b4 = 0
}
else{
extinct_b4 = 1
}
if(sumAfter > h_r_critical){
extinct_after = 0
}
else{
extinct_after = 1
}
cat(file=zz,GopherusTb$Lat[i],',',GopherusTb$Lon[i],',',sumBefore,',',sumAfter,',',extinct_b4,',',extinct_after,'\n')
# this cat correctly makes a comma delimited file with Lat, Lon, and h_restriction and extinction for before (1975) and after (2080) cases
}
close(zz)
GopherusTb
}