#################### Estimating Gull and Tern Avoidance Rates
#######################
#######################
## load libraries

library(lubridate) ### needed to generate sequence of dates from start to end of each study period
library(truncnorm) ### truncated normal distributions for flight speeds etc in sCRM
library(foreach) ### to run code in parallel to speed up extended sCRM
library(doParallel) ### to run code in parallel to speed up extended sCRM
library(ggplot2) ### to plot outputs

#######################################################
############################################ Set working directory where data & outputs will be stored
#######################################################

setwd("C://Users/aonghais/Desktop/NE Avoidance Rates")
options("scipen"=100, "digits"=8)

#################################################
############################### read in data
#################################################

coll.dat = read.table("CollisionData.csv", header = T, sep = ",")  ### collated collision & survey data
flightheight = read.table("Flight Height Data.csv", header = T, sep = ",") ### modelled proportion of birds at collision risk height

BH.boot  = read.table("BH_Boot.csv", header = T, sep = ",") ### Black-headed Gull Flight height bootstraps
CG.boot  = read.table("CG_Boot.csv", header = T, sep = ",") ### Common Gull Flight height bootstraps
CT.boot  = read.table("CT_Boot.csv", header = T, sep = ",") ### Common Tern Flight height bootstraps
GB.boot  = read.table("GB_Boot.csv", header = T, sep = ",") ### Great Black-backed Gull Flight height bootstraps
HG.boot  = read.table("BH_Boot.csv", header = T, sep = ",") ### Herring Gull Flight height bootstraps
KI.boot  = read.table("KI_Boot.csv", header = T, sep = ",") ### Kittiwake Flight height bootstraps
LB.boot  = read.table("LB_Boot.csv", header = T, sep = ",") ### Lesser Black-backed Gull Flight height bootstraps
NE.boot  = read.table("NE_Boot.csv", header = T, sep = ",") ### Sandwich Tern Flight height bootstraps
LG.boot  = read.table("LG_Boot.csv", header = T, sep = ",") ### Little Gull Flight height bootstraps

########################################### FUNCTIONS

################ 
################ function for probability of collision from Basic Band model
################

### use k =1 & 3 blades as default

pcoll = function(k =1, no.blades = 3, max.chord, pitch, rotor.rad, rotor.period,wingspan, bird.length, aspect, flap.glide, bird.speed,prop.flights.upwind,radius,chord) {

	##### p(collision) data

	alpha = bird.speed * rotor.period / (radius * 2 * rotor.rad * pi)

	alphaLTaspect = which(alpha<aspect)
	
	### upwind p(collision)

	up_col_length = abs(max.chord * chord * sin(pitch*pi/180) + (alpha * max.chord * chord * cos(pitch*pi/180)))
	up_col_length2 = (wingspan * alpha) + up_col_length

	up_col_length2[alphaLTaspect] = bird.length+ up_col_length[alphaLTaspect]

	up_p_collision = (no.blades/rotor.period)*up_col_length2/bird.speed
	up_p_collision[up_p_collision > 1] = 1

	overall_up_p_collision = 2 * (sum(radius[1:19] * up_p_collision[1:19]) + up_p_collision[20]/2) * 0.05


	### downwind p(collision)

	down_col_length = abs(-max.chord * chord * sin(pitch*pi/180) + (alpha * max.chord * chord * cos(pitch*pi/180)))
	down_col_length2 = (wingspan * alpha) + down_col_length

	down_col_length2[alphaLTaspect] = bird.length + down_col_length[alphaLTaspect]

	down_p_collision = (no.blades/rotor.period)*down_col_length2/bird.speed
	down_p_collision[down_p_collision > 1] = 1

	overall_down_p_collision = 2 * (sum(radius[1:19] * down_p_collision[1:19]) + down_p_collision[20]/2) * 0.05

	return(((1-prop.flights.upwind) * overall_down_p_collision) + (prop.flights.upwind * overall_up_p_collision))

	}	

##########################################################
##########################################################
##########################################################
##### estimate the number of hours daylight & night during each study

nhoursday = function(start, end, latitude) {
			day.seq = seq(start, end, by = "day")
			day.seq = yday(day.seq)
			
			DayLength = vector()
			NightLength = vector()

			for(zz in 1:length(day.seq)){
		
			P = asin(0.39795 * cos (0.2163108 + 2 * atan (0.9671396 * tan(0.0086 * (day.seq[zz] - 186)))))
			DayLength[zz] = 24 - (24/pi) * acos((sin(0.8333*pi/180) + sin(latitude * pi / 180) * sin (P)) / (cos(latitude *pi /180 ) * cos(P)))
			NightLength[zz] = 24 - DayLength[zz]

				}

			DayLength = sum(DayLength)
			NightLength = sum(NightLength)

			return(list(DayLength, NightLength))

			}


##########################################################
######################################################
######################################################

######################################
###################################### Delta method for estimating SD around a ratio
######################################

std.dev = function(b,c){

var.b <- var(b)
var.c <- var(c)
cov.bc <- cov(b,c)

a <- sum(b)/sum(c)

pd.a.b <- 1/sum(c)
pd.a.c <- -sum(b)/sum(c)^2

pd.mat  <- matrix(c(pd.a.b, pd.a.c), nrow=1)
vcv.mat <- matrix(c(var.b, cov.bc, cov.bc, var.c), nrow = 2, byrow = TRUE)
mat1  <- pd.mat %*% vcv.mat
var.a <- mat1   %*% t(pd.mat)

return(sqrt(var.a))

}



############################## Functions for extended Band Collision Risk Model
############ Calculate collision probability at any given point along the rotor blade

pcollex = function (r,phi, updown, Rotor.Diameter, Blade.Width, Rotor.Speed, Pitch, flapglide, wingspan, birdlength, speed, Blades,coverC = coverC) {
 
if(r>1) cell_val =21
if(r<=1) which(coverC$rad == abs(round(ceiling(r*100)/5) * 5)/100) -> cell_val

 upper = coverC[cell_val,1]
 lower = coverC[cell_val - 1, 1]

 p = (r - lower) / (upper-lower)
 c = coverC[cell_val-1, 2] + p * (coverC[cell_val,2] - coverC[cell_val-1,2])


 radius = Rotor.Diameter/2 * r

if(radius ==0) return(1)

if(radius !=1) ({

 chord = Blade.Width * c
 omega = Rotor.Speed * 2 * pi /60
 pitch = Pitch * pi /180
 phi = phi * pi /180
 
if(updown == "up") direction = 1
if(updown != "up") direction = -1
 
 if(flapglide == "flap")  wingspan-> Wingspan2 
 if(flapglide != "flap") wingspan * abs(cos(phi))->  Wingspan2
 
  multiplier = Blades * omega /(2 * pi * speed)
 
 alpha = speed/(radius * omega)
 
 CollideLength_1 = abs(direction * chord * sin(pitch) + alpha * chord * cos(pitch))
 
 

if(birdlength > Wingspan2 * alpha) birdlength -> CollideLength2
 if(birdlength <= Wingspan2 * alpha) Wingspan2 * alpha -> CollideLength2
 
return(multiplier * (CollideLength_1 + CollideLength2))

})
 
 }

################# Finally calculate collision integral to account for risk and birds both varying with height

coll.int = function(HD.y,HD.d.y, FluxInt, Blades, HubHeight, Rotor.Diameter, Blade.Width, Rotor.Speed,Pitch,speed,flapglide,birdlength,wingspan, FH.dat, Prop_Upwind,coverC = coverC){

CollMinUP = HD.d.y[1] * xrisksum2(HD.y[1], 0.05, "up", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 ### risk at lowest point on rotor blade
CollIntUP = CollMinUP + HD.d.y[41] * xrisksum2(HD.y[41], 0.05, "up", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 ### risk at highest point on rotor blade

for (i in 2:40) {CollIntUP = CollIntUP + HD.d.y[i] * xrisksum2(HD.y[i], 0.05, "up", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) } #### Fill in intermediate heights
 
CollIntUP = CollIntUP * 0.05 * (2/pi)

## Down wind

CollMinDown = HD.d.y[1] * xrisksum2(HD.y[1], 0.05, "down", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 ### risk at lowest point on rotor blade
CollIntDown = CollMinDown + HD.d.y[41] * xrisksum2(HD.y[41], 0.05, "down", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 ### risk at highest point on rotor blade

for (i in 2:40) { CollIntDown = CollIntDown + HD.d.y[i] * xrisksum2(HD.y[i], 0.05, "down", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) } #### Fill in intermediate heights
 

CollIntDown = CollIntDown * 0.05 * (2/pi)

## Average Collision Integral

CollInt = (Prop_Upwind * CollIntUP) + ((1-Prop_Upwind) * CollIntDown)

return(CollInt)

}






###### Function 3


pcollxy = function(x,y,updown, Rotor.Diameter, Blade.Width, Rotor.Speed, Pitch, flapglide, wingspan, birdlength, speed, Blades,coverC = coverC) {

			r = (x * x + y * y) ^ 0.5

				if( y == 0 & x >= 0) pi /2 -> phi
				if( y == 0 & x < 0) -pi/2 -> phi
				if(y != 0)	atan (x/y) -> phi
				if( y <0) (phi + pi)*180/pi -> phi2
				if(y >= 0) phi*180/pi -> phi2
			
			pcollex(r,phi2,updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC)


}

#### Function 3


xrisksum2 = function (y, xinc,updown, Rotor.Diameter, Blade.Width, Rotor.Speed, Pitch, flapglide, wingspan, birdlength, speed, Blades,coverC = coverC) {

		xmax = (1- y * y) ^ 0.5
		imax = as.integer (xmax/xinc)

		risk = (pcollxy (imax * xinc, y, updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 + pcollxy(xmax, y, updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC)/2) * (xmax - imax * xinc)
		
		if(imax<=0) return(2*risk)		
		
		if(imax>0)({
		
		risk2 = risk + (pcollxy(0,y,updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 + pcollxy(imax * xinc, y, updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC)/2) * xinc
	
		for (i in 1: (imax - 1)) {
		
			risk2 = risk2 + pcollxy(i * xinc, y, updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) * xinc

					}
					
					 return(2*risk2)

					})

	


	}


###################################### Functions for Christie & Urquhart extension to Band model

#### Step 1

PrC = function(Radius, Blades, pitch, chord.width, angvel, length, wingspan, delta, BF10G11, WGamma, WBSpeed, RelRad, RelChord, theta, c, s, Tarc, RotorV) {

		## vector velocity of bird

		vbx = -WBSpeed * sin(delta)
		vby = 0
		vbz = WBSpeed * cos(delta)
		
		## vector of chord at theta = 0

		cx0 = RelChord[1] * chord.width * cos(pitch)
		cy0 = 0
		cz0 = RelChord[1] * chord.width * sin(pitch)

		## vector length of bird

		blx = -length * sin (WGamma)
		bly = 0 
		blz = length * cos(WGamma)

		## vector wingspan of gliding bird

		bwgx = wingspan * cos(WGamma)
		bwgy = 0
		bwgz = wingspan * sin (WGamma)

		## vector velocity of rotor		

		vrx = RotorV[1] * c
		vry = -RotorV[1] * s
		vrz = 0

		## vector velocity of bird relative to blade

		u = vbx-vrx
		v = vby-vry
		w= vbz - vrz
	
		if (w ==0 | is.na(w)) w = 0.000000001  ## so we don't divide by 0

		## vector of chord rotated through theta	

		cx = cx0 * c
		cy = -s * cx0
		cz = cz0


		l = cx+blx
		m = cy + bly
		n = cz + blz

		L1 = abs(((w*l-u*n)*c+(v*n-w*m)*s)/w)


		l2 = cx-blx
		m2= cy-bly
		n2 = cz-blz

		L2 = abs(((w*l2-u*n2)*c+(v*n2-w*m2)*s)/w)


		l3 = cx + bwgx
		m3 = cy+bwgy
		n3 = cz + bwgz

		L3 = abs(((w*l3-u*n3)*c+(v*n3-w*m3)*s)/w)

		l4 = cx-bwgx
		m4 = cy-bwgy
		n4 = cz-bwgz

		L4 = abs(((w*l4-u*n4)*c+(v*n4-w*m4)*s)/w)

		## wingspan of a flapping bird rotated

		p = c* cos(WGamma)
		q = -s * cos(WGamma)
		r = c*sin(WGamma)

		k = wingspan/(p^2+q^2+r^2)^0.5
		bwfx = k*p
		bwfy = k*q
		bwfz = k*r

		l5 = cx+bwfx
		m5 = cy+bwfy
		n5 = cz+bwfz

		L5 = abs(((w*l5-u*n5)*c+(v*n5-w*m5)*s)/w)

		l6 = cx-bwfx
		m6 = cy-bwfy
		n6 = cz-bwfz

		L6 = abs(((w*l6-u*n6)*c+(v*n6-w*m6)*s)/w)

## is bird flapping or gliding?
# 0 = flapping
# 1 = gliding

if(BF10G11 == 0) ColLength = pmax(L1, L2, L3, L4, L5, L6)
if(BF10G11 == 1) ColLength = pmax(L1, L2, L3, L4)


Ptot = min(ColLength[1]/Tarc[1], 1) + (2*min(ColLength[2]/Tarc[1], 1)) + (2*min(ColLength[3]/Tarc[1], 1)) + (2*min(ColLength[4]/Tarc[1], 1)) + (2*min(ColLength[5]/Tarc[1], 1)) + (2*min(ColLength[6]/Tarc[1], 1)) + (min(ColLength[7]/Tarc[1], 1))

PrC = Ptot/12
return(PrC)
} 

################ Step 2

mean.coll.risk = function(RelRad, RelChord, simpsons.weight, theta, c,s,Radius, Blades, pitch, chord.width, period, length, wingspan, speed, BF10G11, min.angle, max.angle, wind.speed){

angle.head = seq(min.angle,max.angle,5)
cor.pcoll = vector()

for(abc in 1:length(angle.head)) {

## transform parameters for model

delta = angle.head[abc]*pi/180
WGamma = delta - asin(wind.speed/speed*sin(delta))
WBSpeed = (wind.speed^2 + speed^2 - 2*wind.speed * speed *cos(WGamma))^0.5
angvel = 2*pi/period
Tarc = 2* pi * RelRad * Radius/Blades
RotorV = RelRad * Radius * angvel

col.vec = vector()

for(i in 1:20){ col.vec[i] = 2 * 0.05 *RelRad[i] *PrC(Radius, Blades, pitch, chord.width,angvel, length, wingspan,delta,BF10G11, WGamma, WBSpeed, RelRad[i],RelChord[i],theta,c, s, Tarc[i], RotorV[i])}

######## average collision risk along blade, weighted by Simpson's Rule to better approximate ara under the curve

col.entry = sum(simpsons.weight*col.vec)/3

######## Correct to reflect that area visible to bird as it approaches 
######## is smaller than area visible to bird approaching perpendicularly

## calculate frontal area

AFront = pi * Radius^2

## calculate area covered by side of the turbine

ASide = (mean(RelChord) * chord.width * sin(pitch) + wingspan) * 2 * Radius

## Calculate the area visible

AVis = AFront * abs(cos(delta)) + ASide * abs(sin(delta))

## now correct probability of collision

cor.pcoll[abc] = col.entry * AVis/AFront

}

return(cor.pcoll)

}

############################################
############################################






###############################################################
###############################################################
###############################################################
###### Start formatting, need to get the number of birds passing through a turbine rotor sweep over the survey period

## first  step is to get total number of hours daylight and night over each survey period
# format start date and end date into dates

coll.dat$start.date = as.Date(coll.dat$start.date, format = "%d/%m/%Y")
coll.dat$end.date = as.Date(coll.dat$end.date, format = "%d/%m/%Y")

##### create columns for number of ours daylight & night in each survey			

coll.dat$HoursDay = 0
coll.dat$HoursNight = 0

### now use this to estimate n hours daylight & night
for(i in 1:nrow(coll.dat)) { 

			if(!is.na(coll.dat$end.date[i]) )({ 
			hours = nhoursday(coll.dat$start.date[i], coll.dat$end.date[i],coll.dat$Latitude[i])
			coll.dat$HoursDay[i] = hours[1]
			coll.dat$HoursNight[i] = hours[2]
			 })
			}


#################################################################################
#################################################################################
#################################################################################
################# First work out avoidance rates for Basic Band (2012) model
#################################################################################
#################################################################################

#### need to get area of survey window at collision risk height


coll.dat$area.survey.window = coll.dat$width.of.survey.area * coll.dat$rotor.diameter
coll.dat$total.rotor.frontal.area = coll.dat$number.of.turbines * pi * (coll.dat$rotor.diameter/2)^2

########### For sites where we do not have a proportion of birds at risk height from site-specific survey data, we need to estimate this from the 
########### Johnston et al. modelled distributions

coll.dat$FH.modelled = 0

for(i in 1:nrow( coll.dat)) {


sp.fh.dat = flightheight[ , grepl( paste(coll.dat$Species.for.Flight.Height.Distribution[i], ".est", sep = ""), names( flightheight ) ) ]
min.rotor.height = coll.dat$hub.height[i] - (0.5 * coll.dat$rotor.diameter[i])
max.rotor.height = coll.dat$hub.height[i] + (0.5 * coll.dat$rotor.diameter[i])
coll.dat$FH.modelled[i] = sum(sp.fh.dat[min.rotor.height:max.rotor.height])

}

#### get proportion of birds at risk height for sites with no site-specific data

coll.dat$X..flights.at.rotor.height[is.na(coll.dat$X..flights.at.rotor.height)] = coll.dat$FH.modelled[is.na(coll.dat$X..flights.at.rotor.height)]

#### Now work out the probability of collision at each site
## vectors needed for the Band model

radius = seq(0.05, 1, 0.05)
chord = c(0.73, 0.79, 0.88, 0.96, 1, 0.98, 0.92, 0.85, 0.80, 0.75,0.7, 0.64,0.58,0.52,0.47,0.41,0.37,0.30,0.24, 0)

for(i in 1:nrow(coll.dat)){

coll.dat$Pcoll[i] = pcoll(k =1, no.blades = 3, max.chord = coll.dat$blade.width[i], pitch = coll.dat$rotor.pitch[i], rotor.rad = coll.dat$rotor.diameter[i]/2, rotor.period = 60/coll.dat$rotor.speed[i],wingspan = coll.dat$wingspan[i], bird.length = coll.dat$BirdLength[i], aspect = coll.dat$BirdLength[i]/coll.dat$wingspan[i], flap.glide = 0, bird.speed = coll.dat$flight.speed[i],prop.flights.upwind = 0.5,radius,chord)

}

########## Now need to add column to main dataset shown total number of birds passing at collision risk height
coll.dat$HoursDay = as.numeric(coll.dat$HoursDay)
coll.dat$HoursNight = as.numeric(coll.dat$HoursNight)

coll.dat$passage.rate[!is.na(coll.dat$Number.of.birds.recorded)] = coll.dat$Number.of.birds.recorded[!is.na(coll.dat$Number.of.birds.recorded)]/coll.dat$N.Hours.monitoring[!is.na(coll.dat$Number.of.birds.recorded)]

coll.dat$total.flux[!is.na(coll.dat$passage.rate)]  = (coll.dat$passage.rate[!is.na(coll.dat$passage.rate)]  * coll.dat$HoursNight[!is.na(coll.dat$passage.rate)]  * coll.dat$Nocturnal.Correction[!is.na(coll.dat$passage.rate)]) + 
					(coll.dat$passage.rate[!is.na(coll.dat$passage.rate)] * coll.dat$HoursDay[!is.na(coll.dat$passage.rate)] ) 

coll.dat$total.flux[!is.na(coll.dat$density)] = (3600/(coll.dat$width.of.survey.area[!is.na(coll.dat$density)]/coll.dat$flight.speed[!is.na(coll.dat$density)])) * coll.dat$HoursDay[!is.na(coll.dat$density)] * (coll.dat$density[!is.na(coll.dat$density)] * coll.dat$Survey.area[!is.na(coll.dat$density)])

coll.dat$flux.CRH = coll.dat$total.flux * coll.dat$X..flights.at.rotor.height * (coll.dat$total.rotor.frontal.area/coll.dat$area.survey.window)

coll.dat$PredColl.NoAvoid.Opt1 = coll.dat$Pcoll *  coll.dat$flux.CRH

#### plot site-specific avoidance rates

Band.Basic.plot = ggplot(coll.dat, aes(x = group, y = 1-(Estimated.Collisions/PredColl.NoAvoid.Opt1)))

### get avoidance rates summarised by species & group

spp = unique(coll.dat$Species)
Basic.Band.Spp.Avoid = vector()
Basic.Band.Spp.AvoidSD = vector()
Basic.Band.Spp.AvoidUpCI = vector()
Basic.Band.Spp.AvoidLowCI = vector()

for(i in 1:length(spp)) {

Basic.Band.Spp.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == spp[i]])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$Species == spp[i]]))
Basic.Band.Spp.AvoidSD[i] = std.dev(coll.dat$Estimated.Collisions[coll.dat$Species == spp[i]], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$Species == spp[i]])
Basic.Band.Spp.AvoidUpCI[i] = Basic.Band.Spp.Avoid[i] + (1.96 * Basic.Band.Spp.AvoidSD[i])
Basic.Band.Spp.AvoidLowCI[i] = Basic.Band.Spp.Avoid[i] - (1.96 * Basic.Band.Spp.AvoidSD[i])

}

Basic.Band.Avoid.Spp.Table = data.frame(Species = spp, Avoidance = Basic.Band.Spp.Avoid, AvoidanceSD = Basic.Band.Spp.AvoidSD, AvoidanceLCL = Basic.Band.Spp.AvoidLowCI, AvoidanceUCL = Basic.Band.Spp.AvoidUpCI)


group = unique(coll.dat$group)
Basic.Band.Group.Avoid = vector()
Basic.Band.Group.AvoidSD = vector()
Basic.Band.Group.AvoidUpCI = vector()
Basic.Band.Group.AvoidLowCI = vector()

for(i in 1:length(group)) {

Basic.Band.Group.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == group[i]])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group == group[i]]))
Basic.Band.Group.AvoidSD[i] = std.dev(coll.dat$Estimated.Collisions[coll.dat$group == group[i]], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group == group[i]])
Basic.Band.Group.AvoidUpCI[i] = Basic.Band.Group.Avoid[i] + (1.96 * Basic.Band.Group.AvoidSD[i])
Basic.Band.Group.AvoidLowCI[i] = Basic.Band.Group.Avoid[i] - (1.96 * Basic.Band.Group.AvoidSD[i])

}

Basic.Band.Group.Avoid[4] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group %in% c("small gull", "large gull", "gull")]))
Basic.Band.Group.AvoidSD[4] = std.dev(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group  %in% c("small gull", "large gull", "gull")])
Basic.Band.Group.AvoidUpCI[4] = Basic.Band.Group.Avoid[4] + (1.96 * Basic.Band.Group.AvoidSD[4])
Basic.Band.Group.AvoidLowCI[4] = Basic.Band.Group.Avoid[4] - (1.96 * Basic.Band.Group.AvoidSD[4])



Basic.Band.Avoid.Group.Table = data.frame(Species = group, Avoidance = Basic.Band.Group.Avoid, AvoidanceSD = Basic.Band.Group.AvoidSD, AvoidanceLCL = Basic.Band.Group.AvoidLowCI, AvoidanceUCL = Basic.Band.Group.AvoidUpCI)

Basic.Band.Avoid.Group.Table = rbind(Basic.Band.Avoid.Spp.Table,Basic.Band.Avoid.Group.Table)[c(1,2,3,4,5,9,10,17,18,21,22,23,24),]

Basic.Band.Avoid.Group.Table$ARs = paste(round(Basic.Band.Avoid.Group.Table[,2],4), " (",round(Basic.Band.Avoid.Group.Table[,3],4), "; ", round(Basic.Band.Avoid.Group.Table[,4],4), " - ", round(Basic.Band.Avoid.Group.Table[,5],4),")", sep = "") 

##############################################################
##############################################################
##############################################################
############ Extended Band Model

##### try parallelizing to improve speed

comb <- function(...) {
  mapply('cbind', ..., SIMPLIFY=FALSE)
}


if(nrow(coll.dat) >= detectCores()) cl <- makeCluster(detectCores())
if(nrow(coll.dat) < detectCores()) cl = makeCluster(niter)

registerDoParallel(cl)

#################################################
#################################################

HD.y = round(seq(-1,1,0.05),2)
rad = round(seq(0,1,0.05),2) ###### very convoluted, but necessary due to the way R treats seq...
c = c(0.69,0.73,0.79,0.88,0.96,1,0.98,0.92,0.85,0.8,0.75,0.7,0.64,0.58,0.52,0.47,0.41,0.37,0.3,0.24,0)
coverC = data.frame(cbind(rad,c))

coll.int.out = foreach(i = 1:nrow(coll.dat), .combine = "comb") %dopar% {

FH.dat = flightheight[ , grepl( paste(coll.dat$Species.for.Flight.Height.Distribution[i], ".est", sep = ""), names( flightheight ) ) ]

height = coll.dat$hub.height[i] + (HD.y * (coll.dat$rotor.diameter[i]/2))
HD.d.y = (((height-floor(height)) * (FH.dat[ceiling(height)+1] - FH.dat[floor(height)+1])) + FH.dat[floor(height)+1]) * coll.dat$rotor.diameter[i]/2
FluxMin = HD.d.y[1] * (2 * ((1 - HD.y[1] * HD.y[1])^0.5)) / 2 ### risk at lowest point on rotor blade
FluxInt = FluxMin + HD.d.y[41] * (2 * ((1 - HD.y[41] * HD.y[41])^0.5)) / 2 ### risk at highest point on rotor blade
for (zz in 2:40) { FluxInt = FluxInt + HD.d.y[zz] * (2 * ((1 - HD.y[zz] * HD.y[zz])^0.5)) } #### Fill in intermediate heights
FluxInt = FluxInt * 0.05 * (2/pi)


out =coll.int(	HD.y = HD.y, HD.d.y = HD.d.y, FluxInt = FluxInt,
			Blades = 3, 
			HubHeight =  coll.dat$hub.height[i], 
			Rotor.Diameter = coll.dat$rotor.diameter[i], 
			Blade.Width = coll.dat$blade.width[i], 
			Rotor.Speed = coll.dat$rotor.speed[i],
			Pitch = coll.dat$rotor.pitch[i],
			speed = coll.dat$flight.speed[i],
			flapglide = "flap",
			birdlength = coll.dat$BirdLength[i],
			wingspan = coll.dat$wingspan[i], 
			FH.dat = FH.dat,
			Prop_Upwind = 0.5, coverC = coverC)

	list(out)

}


stopCluster(cl)

coll.dat$coll.int = coll.int.out[[1]][1,]

coll.dat$flux.All = coll.dat$total.flux * (coll.dat$total.rotor.frontal.area/coll.dat$area.survey.window)

coll.dat$PredColl.NoAvoid.Opt3 = coll.dat$flux.All * coll.dat$coll.int


Extended.Band.Spp.Avoid = vector()
Extended.Band.Spp.AvoidSD = vector()
Extended.Band.Spp.AvoidUpCI = vector()
Extended.Band.Spp.AvoidLowCI = vector()

for(i in 1:length(spp)) {

Extended.Band.Spp.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == spp[i]])/sum(coll.dat$PredColl.NoAvoid.Opt3[coll.dat$Species == spp[i]]))
Extended.Band.Spp.AvoidSD[i] = std.dev(coll.dat$Estimated.Collisions[coll.dat$Species == spp[i]], coll.dat$PredColl.NoAvoid.Opt3[coll.dat$Species == spp[i]])
Extended.Band.Spp.AvoidUpCI[i] = Extended.Band.Spp.Avoid[i] + (1.96 * Extended.Band.Spp.AvoidSD[i])
Extended.Band.Spp.AvoidLowCI[i] = Extended.Band.Spp.Avoid[i] - (1.96 * Extended.Band.Spp.AvoidSD[i])

}

Extended.Band.Avoid.Spp.Table = data.frame(Species = spp, Avoidance = Extended.Band.Spp.Avoid, AvoidanceSD = Extended.Band.Spp.AvoidSD, AvoidanceLCL = Extended.Band.Spp.AvoidLowCI, AvoidanceUCL = Extended.Band.Spp.AvoidUpCI)

Extended.Band.Group.Avoid = vector()
Extended.Band.Group.AvoidSD = vector()
Extended.Band.Group.AvoidUpCI = vector()
Extended.Band.Group.AvoidLowCI = vector()

for(i in 1:length(group)) {

Extended.Band.Group.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == group[i]])/sum(coll.dat$PredColl.NoAvoid.Opt3[coll.dat$group == group[i]]))
Extended.Band.Group.AvoidSD[i] = std.dev(coll.dat$Estimated.Collisions[coll.dat$group == group[i]], coll.dat$PredColl.NoAvoid.Opt3[coll.dat$group == group[i]])
Extended.Band.Group.AvoidUpCI[i] = Extended.Band.Group.Avoid[i] + (1.96 * Extended.Band.Group.AvoidSD[i])
Extended.Band.Group.AvoidLowCI[i] = Extended.Band.Group.Avoid[i] - (1.96 * Extended.Band.Group.AvoidSD[i])

}


Extended.Band.Group.Avoid[4] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt3[coll.dat$group %in% c("small gull", "large gull", "gull")]))
Extended.Band.Group.AvoidSD[4] = std.dev(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")], coll.dat$PredColl.NoAvoid.Opt3[coll.dat$group  %in% c("small gull", "large gull", "gull")])
Extended.Band.Group.AvoidUpCI[4] = Extended.Band.Group.Avoid[4] + (1.96 * Extended.Band.Group.AvoidSD[4])
Extended.Band.Group.AvoidLowCI[4] = Extended.Band.Group.Avoid[4] - (1.96 * Extended.Band.Group.AvoidSD[4])



Extended.Band.Avoid.Group.Table = data.frame(Species = group, Avoidance = Extended.Band.Group.Avoid, AvoidanceSD = Extended.Band.Group.AvoidSD, AvoidanceLCL = Extended.Band.Group.AvoidLowCI, AvoidanceUCL = Extended.Band.Group.AvoidUpCI)

Extended.Band.Avoid.Table = rbind(Extended.Band.Avoid.Spp.Table,Extended.Band.Avoid.Group.Table)[c(1,2,3,4,5,9,10,17,18,21,22,23,24),]
Extended.Band.Avoid.Table$ARs = paste(round(Extended.Band.Avoid.Table[,2],4), " (",round(Extended.Band.Avoid.Table[,3],4), "; ", round(Extended.Band.Avoid.Table[,4],4), " - ", round(Extended.Band.Avoid.Table[,5],4),")", sep = "") 

###########################
###########################


#########################################################################
#########################################################################
#########################################################################

######################### Now for the stochastic avoidance rates
######################### First step is to simultate the number of birds passing through the rotor swept areas

######  Now need to work out the number of birds passing through each wind farm during the study period
nsim = 1000

flux.dat = matrix(nrow = nrow(coll.dat), ncol = nsim)

##### for sites where we already have total flux (e.g. Oosterbierum)

flux.dat[!is.na(coll.dat$total.flux),] = coll.dat$total.flux[!is.na(coll.dat$total.flux)]

### For sites where we have  a number of birds recorded and a number of hours monitoring

coll.dat$passage.rate[!is.na(coll.dat$Number.of.birds.recorded)] = coll.dat$Number.of.birds.recorded[!is.na(coll.dat$Number.of.birds.recorded)]/coll.dat$N.Hours.monitoring[!is.na(coll.dat$Number.of.birds.recorded)]

for(i in 1:nrow(coll.dat[!is.na(coll.dat$passage.rate),])) {

nbirds.pass = rpois(nsim,coll.dat$passage.rate[!is.na(coll.dat$passage.rate)][i]) 

flux.dat[!is.na(coll.dat$passage.rate),] [i,]=  (nbirds.pass * coll.dat$HoursNight[!is.na(coll.dat$passage.rate)][[i]]  * coll.dat$Nocturnal.Correction[!is.na(coll.dat$passage.rate)] [i]) + (nbirds.pass * coll.dat$HoursDay[!is.na(coll.dat$passage.rate)] [[i]]) 

}

#### get total number of flights through rotor sweep NOT adjusted for proportion at risk height

for(i in 1:ncol(flux.dat)) { flux.dat[,i] = flux.dat[,i] / coll.dat$area.survey.window * coll.dat$total.rotor.frontal.area }

##### Do the same for sites where bird data are availabile as densities rather than counts

for(i in 1:nrow(flux.dat[!is.na(coll.dat$density),])) {

flux.dat[!is.na(coll.dat$density),][i,] = (3600/(coll.dat$width.of.survey.area[!is.na(coll.dat$density)][i]/rtruncnorm(n = nsim, a = 4, coll.dat$flight.speed[!is.na(coll.dat$density)][i]))) * coll.dat$HoursDay[!is.na(coll.dat$density)][i] * (rpois(nsim, coll.dat$density[!is.na(coll.dat$density)][i]) * coll.dat$Survey.area[!is.na(coll.dat$density)][i])

				
}

#### We now have 1000 estimates of the number of birds passing through turbine rotors, need to correct for PCH and probability of collision
#### Now work out the probability of collision at each site
## vectors needed for the Band model

##### now simulate Pcoll

Pcoll.dat = matrix(nrow = nrow(coll.dat), ncol = nsim)

for(i in 1:nrow(Pcoll.dat)) {
for(y in 1:ncol(Pcoll.dat)){

simRotorSpeed = rnorm(1,coll.dat$rotor.speed[i],0.5)
simRotorPitch = rnorm(1,coll.dat$rotor.pitch[i],0.1)
simBodyLength = rnorm(1, coll.dat$BirdLength[i], coll.dat$BirdLengthSD[i])
simWingspan = rnorm(1, coll.dat$wingspan[i], coll.dat$wingspanSD[i])
simFlightSpeed = rtruncnorm(n = 1,a = 4, mean = coll.dat$flight.speed[i], sd = coll.dat$flightspeedSD[i])

Pcoll.dat[i,y] = pcoll(k =1, no.blades = 3, max.chord = coll.dat$blade.width[i], pitch = simRotorPitch, rotor.rad = coll.dat$rotor.diameter[i]/2, rotor.period = 60/simRotorSpeed,wingspan = simWingspan, bird.length = simBodyLength, aspect = simBodyLength/simWingspan, flap.glide = 0, bird.speed = simFlightSpeed,prop.flights.upwind = 0.5,radius,chord)

}}

### now estimate collisions in the absence of avoidance

coll.no.avoid = matrix(nrow = nrow(coll.dat), ncol = nsim)

for( i in 1:ncol(coll.no.avoid)){  coll.no.avoid[,i] = Pcoll.dat[,i] * flux.dat[,i] * coll.dat$X..flights.at.rotor.height  }

## now bootstrapped avoidance rates

HG.Basic.sCRM.Avoid = vector()
LB.Basic.sCRM.Avoid = vector()
GB.Basic.sCRM.Avoid = vector()
KI.Basic.sCRM.Avoid = vector()
BH.Basic.sCRM.Avoid = vector()
CG.Basic.sCRM.Avoid = vector()
LT.Basic.sCRM.Avoid = vector()
ST.Basic.sCRM.Avoid = vector()
CT.Basic.sCRM.Avoid = vector()
SG.Basic.sCRM.Avoid = vector()
LG.Basic.sCRM.Avoid = vector()
AG.Basic.sCRM.Avoid = vector()
AT.Basic.sCRM.Avoid = vector()

for (i in 1:ncol(coll.no.avoid)){

HG.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "herring gull"])/sum(coll.no.avoid[coll.dat$Species == "herring gull", i]))
LB.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "lesser black-backed gull"])/sum(coll.no.avoid[coll.dat$Species == "lesser black-backed gull", i]))
GB.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "great black-backed gull"])/sum(coll.no.avoid[coll.dat$Species == "great black-backed gull", i]))
KI.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "kittiwake"])/sum(coll.no.avoid[coll.dat$Species == "kittiwake", i]))
BH.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "black-headed gull"])/sum(coll.no.avoid[coll.dat$Species == "black-headed gull", i]))
CG.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common gull"])/sum(coll.no.avoid[coll.dat$Species == "common gull", i]))
LT.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern"])/sum(coll.no.avoid[coll.dat$Species == "little tern", i]))
ST.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern"])/sum(coll.no.avoid[coll.dat$Species == "sandwich tern", i]))
CT.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern"])/sum(coll.no.avoid[coll.dat$Species == "common tern", i]))
SG.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "small gull"])/sum(coll.no.avoid[coll.dat$group == "small gull", i]))
LG.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "large gull"])/sum(coll.no.avoid[coll.dat$group == "large gull", i]))
AG.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("large gull", "small gull", "gull")])/sum(coll.no.avoid[coll.dat$group %in% c("gull","large gull", "small gull"), i]))
AT.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern"])/sum(coll.no.avoid[coll.dat$group == "tern", i]))

}




Basic.sCRM.Avoid.Table = data.frame(Species = c("Herring Gull", "Lesser Black-backed Gull", "Great Black-Backed Gull", "Kittiwake", "Black-headed Gull", "Common Gull", "Little Tern", "Sandwich Tern", "Common Tern", "Small Gull", "Large Gull", "All Gull", "All Tern"), 
							Avoidance = c(median(HG.Basic.sCRM.Avoid),median(LB.Basic.sCRM.Avoid),median(GB.Basic.sCRM.Avoid),median(KI.Basic.sCRM.Avoid[KI.Basic.sCRM.Avoid >0]),median(BH.Basic.sCRM.Avoid), median(CG.Basic.sCRM.Avoid),median(LT.Basic.sCRM.Avoid), median(ST.Basic.sCRM.Avoid),median(CT.Basic.sCRM.Avoid),median(SG.Basic.sCRM.Avoid),median(LG.Basic.sCRM.Avoid),median(AG.Basic.sCRM.Avoid),median(AT.Basic.sCRM.Avoid)),
							AvoidanceSD = c(sd(HG.Basic.sCRM.Avoid),sd(LB.Basic.sCRM.Avoid),sd(GB.Basic.sCRM.Avoid),sd(KI.Basic.sCRM.Avoid[KI.Basic.sCRM.Avoid >0]),sd(BH.Basic.sCRM.Avoid), sd(CG.Basic.sCRM.Avoid),sd(LT.Basic.sCRM.Avoid), sd(ST.Basic.sCRM.Avoid),sd(CT.Basic.sCRM.Avoid),sd(SG.Basic.sCRM.Avoid),sd(LG.Basic.sCRM.Avoid),sd(AG.Basic.sCRM.Avoid),sd(AT.Basic.sCRM.Avoid)),
							AvoidanceLCL = c(quantile(HG.Basic.sCRM.Avoid, 0.025),quantile(LB.Basic.sCRM.Avoid, 0.025),quantile(GB.Basic.sCRM.Avoid, 0.025),quantile(KI.Basic.sCRM.Avoid[KI.Basic.sCRM.Avoid >0], 0.025),quantile(BH.Basic.sCRM.Avoid, 0.025),quantile(CG.Basic.sCRM.Avoid, 0.025),quantile(LT.Basic.sCRM.Avoid, 0.025),quantile(ST.Basic.sCRM.Avoid, 0.025),quantile(CT.Basic.sCRM.Avoid, 0.025),quantile(SG.Basic.sCRM.Avoid, 0.025),quantile(LG.Basic.sCRM.Avoid, 0.025),quantile(AG.Basic.sCRM.Avoid, 0.025),quantile(AT.Basic.sCRM.Avoid, 0.025)), 
							AvoidanceUCL = c(quantile(HG.Basic.sCRM.Avoid, 0.975),quantile(LB.Basic.sCRM.Avoid, 0.975),quantile(GB.Basic.sCRM.Avoid, 0.975),quantile(KI.Basic.sCRM.Avoid[KI.Basic.sCRM.Avoid >0], 0.975),quantile(BH.Basic.sCRM.Avoid, 0.975),quantile(CG.Basic.sCRM.Avoid, 0.975),quantile(LT.Basic.sCRM.Avoid, 0.975),quantile(ST.Basic.sCRM.Avoid, 0.975),quantile(CT.Basic.sCRM.Avoid, 0.975),quantile(SG.Basic.sCRM.Avoid, 0.975),quantile(LG.Basic.sCRM.Avoid, 0.975),quantile(AG.Basic.sCRM.Avoid, 0.975),quantile(AT.Basic.sCRM.Avoid, 0.975)))


Basic.sCRM.Avoid.Table$ARs = paste(round(Basic.sCRM.Avoid.Table[,2],4), " (",round(Basic.sCRM.Avoid.Table[,3],4), "; ", round(Basic.sCRM.Avoid.Table[,4],4), " - ", round(Basic.sCRM.Avoid.Table[,5],4),")", sep = "") 

###### Now  Collision Rates for FCM

HG.FCM.CollRate = vector()
LB.FCM.CollRate = vector()
GB.FCM.CollRate = vector()
KI.FCM.CollRate = vector()
BH.FCM.CollRate = vector()
CG.FCM.CollRate = vector()
LT.FCM.CollRate = vector()
ST.FCM.CollRate = vector()
CT.FCM.CollRate = vector()
SG.FCM.CollRate = vector()
LG.FCM.CollRate = vector()
AG.FCM.CollRate = vector()
AT.FCM.CollRate = vector()

for (i in 1:ncol(coll.no.avoid)){

HG.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "herring gull"])/sum(flux.dat[coll.dat$Species == "herring gull", i]))
LB.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "lesser black-backed gull"])/sum(flux.dat[coll.dat$Species == "lesser black-backed gull", i]))
GB.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "great black-backed gull"])/sum(flux.dat[coll.dat$Species == "great black-backed gull", i]))
KI.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "kittiwake"])/sum(flux.dat[coll.dat$Species == "kittiwake", i]))
BH.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "black-headed gull"])/sum(flux.dat[coll.dat$Species == "black-headed gull", i]))
CG.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common gull"])/sum(flux.dat[coll.dat$Species == "common gull", i]))
LT.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern"])/sum(flux.dat[coll.dat$Species == "little tern", i]))
ST.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern"])/sum(flux.dat[coll.dat$Species == "sandwich tern", i]))
CT.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern"])/sum(flux.dat[coll.dat$Species == "common tern", i]))
SG.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$group == "small gull"])/sum(flux.dat[coll.dat$group == "small gull", i]))
LG.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$group == "large gull"])/sum(flux.dat[coll.dat$group == "large gull", i]))
AG.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("large gull", "small gull", "gull")])/sum(flux.dat[coll.dat$group %in% c("large gull", "small gull", "gull"), i]))
AT.FCM.CollRate[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern"])/sum(flux.dat[coll.dat$group == "tern", i]))

}


FCM.CollRate.Table = data.frame(Species = c("Herring Gull", "Lesser Black-backed Gull", "Great Black-Backed Gull", "Kittiwake", "Black-headed Gull", "Common Gull", "Little Tern", "Sandwich Tern", "Common Tern", "Small Gull", "Large Gull", "All Gull", "All Tern"), 
							Avoidance = c(median(HG.FCM.CollRate),median(LB.FCM.CollRate),median(GB.FCM.CollRate),median(KI.FCM.CollRate[KI.FCM.CollRate <1]),median(BH.FCM.CollRate), median(CG.FCM.CollRate),median(LT.FCM.CollRate), median(ST.FCM.CollRate),median(CT.FCM.CollRate),median(SG.FCM.CollRate),median(LG.FCM.CollRate),median(AG.FCM.CollRate),median(AT.FCM.CollRate)),
							AvoidanceSD = c(sd(HG.FCM.CollRate),sd(LB.FCM.CollRate),sd(GB.FCM.CollRate),sd(KI.FCM.CollRate[KI.FCM.CollRate <1]),sd(BH.FCM.CollRate), sd(CG.FCM.CollRate),sd(LT.FCM.CollRate), sd(ST.FCM.CollRate),sd(CT.FCM.CollRate),sd(SG.FCM.CollRate),sd(LG.FCM.CollRate),sd(AG.FCM.CollRate),sd(AT.FCM.CollRate)),
							AvoidanceLCL = c(quantile(HG.FCM.CollRate, 0.025),quantile(LB.FCM.CollRate, 0.025),quantile(GB.FCM.CollRate, 0.025),quantile(KI.FCM.CollRate[KI.FCM.CollRate <1], 0.025),quantile(BH.FCM.CollRate, 0.025),quantile(CG.FCM.CollRate, 0.025),quantile(LT.FCM.CollRate, 0.025),quantile(ST.FCM.CollRate, 0.025),quantile(CT.FCM.CollRate, 0.025),quantile(SG.FCM.CollRate, 0.025),quantile(LG.FCM.CollRate, 0.025),quantile(AG.FCM.CollRate, 0.025),quantile(AT.FCM.CollRate, 0.025)), 
							AvoidanceUCL = c(quantile(HG.FCM.CollRate, 0.975),quantile(LB.FCM.CollRate, 0.975),quantile(GB.FCM.CollRate, 0.975),quantile(KI.FCM.CollRate[KI.FCM.CollRate <1], 0.975),quantile(BH.FCM.CollRate, 0.975),quantile(CG.FCM.CollRate, 0.975),quantile(LT.FCM.CollRate, 0.975),quantile(ST.FCM.CollRate, 0.975),quantile(CT.FCM.CollRate, 0.975),quantile(SG.FCM.CollRate, 0.975),quantile(LG.FCM.CollRate, 0.975),quantile(AG.FCM.CollRate, 0.975),quantile(AT.FCM.CollRate, 0.975)))


FCM.CollRate.Table$ARs = paste(round(FCM.CollRate.Table[,2],8), " (",round(FCM.CollRate.Table[,3],8), "; ", round(FCM.CollRate.Table[,4],8), " - ", round(FCM.CollRate.Table[,5],8),")", sep = "") 


###########################
###########################
### Bootstrapped extended

if(nrow(coll.dat) >= detectCores()) cl <- makeCluster(detectCores())
if(nrow(coll.dat) < detectCores()) cl = makeCluster(niter)

registerDoParallel(cl)

coll.int.out2 = foreach(z = 1:nsim, .combine = "comb") %dopar% {

out = vector()

for(i in 1:nrow(coll.dat)){

fh.dist.sample = sample(seq(2,201,1),1)

if(coll.dat$Species.for.Flight.Height.Distribution[i] == "BlackHeadedGull") FH.dat = BH.boot[,fh.dist.sample]
if(coll.dat$Species.for.Flight.Height.Distribution[i] == "CommonGull") FH.dat = CG.boot[,fh.dist.sample]
if(coll.dat$Species.for.Flight.Height.Distribution[i] == "LesserBlackBackedGull") FH.dat = LB.boot[,fh.dist.sample]
if(coll.dat$Species.for.Flight.Height.Distribution[i] == "HerringGull") FH.dat = HG.boot[,fh.dist.sample]
if(coll.dat$Species.for.Flight.Height.Distribution[i] == "GreatBlackBackedGull") FH.dat = GB.boot[,fh.dist.sample]
if(coll.dat$Species.for.Flight.Height.Distribution[i] == "Kittiwake") FH.dat = KI.boot[,fh.dist.sample]
if(coll.dat$Species.for.Flight.Height.Distribution[i] == "SandwichTern") FH.dat = NE.boot[,fh.dist.sample]
if(coll.dat$Species.for.Flight.Height.Distribution[i] == "LittleGull") FH.dat = LG.boot[,fh.dist.sample]
if(coll.dat$Species.for.Flight.Height.Distribution[i] == "CommonTern") FH.dat = CT.boot[,fh.dist.sample]

height = coll.dat$hub.height[i] + (HD.y * (coll.dat$rotor.diameter[i]/2))
HD.d.y = (((height-floor(height)) * (FH.dat[ceiling(height)+1] - FH.dat[floor(height)+1])) + FH.dat[floor(height)+1]) * coll.dat$rotor.diameter[i]/2
FluxMin = HD.d.y[1] * (2 * ((1 - HD.y[1] * HD.y[1])^0.5)) / 2 ### risk at lowest point on rotor blade
FluxInt = FluxMin + HD.d.y[41] * (2 * ((1 - HD.y[41] * HD.y[41])^0.5)) / 2 ### risk at highest point on rotor blade
for (zz in 2:40) { FluxInt = FluxInt + HD.d.y[zz] * (2 * ((1 - HD.y[zz] * HD.y[zz])^0.5)) } #### Fill in intermediate heights
FluxInt = FluxInt * 0.05 * (2/pi)


out[i] =coll.int(	HD.y = HD.y, HD.d.y = HD.d.y, FluxInt = FluxInt,
			Blades = 3, 
			HubHeight =  coll.dat$hub.height[i], 
			Rotor.Diameter = coll.dat$rotor.diameter[i], 
			Blade.Width = coll.dat$blade.width[i], 
			Rotor.Speed = rnorm(1,coll.dat$rotor.speed[i],0.5),
			Pitch = rnorm(1, coll.dat$rotor.pitch[i],0.1),
			speed = rnorm(1, coll.dat$flight.speed[i],coll.dat$flightspeedSD[i]),
			flapglide = "flap",
			birdlength = rnorm(1, coll.dat$BirdLength[i],coll.dat$BirdLengthSD[i]),
			wingspan = rnorm(1, coll.dat$wingspan[i],coll.dat$wingspanSD[i]), 
			FH.dat = FH.dat,
			Prop_Upwind = 0.5, coverC = coverC)

			}
	list(out)

}

write.csv(coll.int.out2[[1]], "CollIntSim.csv")


### now estimate collisions in the absence of avoidance

coll.no.avoid.Opt3 = matrix(nrow = nrow(coll.dat), ncol = nsim)

for( i in 1:ncol(coll.no.avoid.Opt3)){  coll.no.avoid.Opt3[,i] = coll.int.out2[[1]][,i] * flux.dat[,i]  }

HG.Extended.sCRM.Avoid = vector()
LB.Extended.sCRM.Avoid = vector()
GB.Extended.sCRM.Avoid = vector()
KI.Extended.sCRM.Avoid = vector()
BH.Extended.sCRM.Avoid = vector()
CG.Extended.sCRM.Avoid = vector()
LT.Extended.sCRM.Avoid = vector()
ST.Extended.sCRM.Avoid = vector()
CT.Extended.sCRM.Avoid = vector()
SG.Extended.sCRM.Avoid = vector()
LG.Extended.sCRM.Avoid = vector()
AG.Extended.sCRM.Avoid = vector()
AT.Extended.sCRM.Avoid = vector()

for (i in 1:ncol(coll.no.avoid.Opt3)){

HG.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "herring gull"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "herring gull", i]))
LB.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "lesser black-backed gull"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "lesser black-backed gull", i]))
GB.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "great black-backed gull"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "great black-backed gull", i]))
KI.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "kittiwake"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "kittiwake", i]))
BH.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "black-headed gull"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "black-headed gull", i]))
CG.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common gull"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "common gull", i]))
LT.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "little tern", i]))
ST.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "sandwich tern", i]))
CT.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "common tern", i]))
SG.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "small gull"])/sum(coll.no.avoid.Opt3[coll.dat$group == "small gull", i]))
LG.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "large gull"])/sum(coll.no.avoid.Opt3[coll.dat$group == "large gull", i]))
AG.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("large gull", "small gull", "gull")])/sum(coll.no.avoid.Opt3[coll.dat$group %in% c("gull","large gull", "small gull"), i]))
AT.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern"])/sum(coll.no.avoid.Opt3[coll.dat$group == "tern", i]))

}




Extended.sCRM.Avoid.Table = data.frame(Species = c("Herring Gull", "Lesser Black-backed Gull", "Great Black-Backed Gull", "Kittiwake", "Black-headed Gull", "Common Gull", "Little Tern", "Sandwich Tern", "Common Tern", "Small Gull", "Large Gull", "All Gull", "All Tern"), 
							Avoidance = c(median(HG.Extended.sCRM.Avoid),median(LB.Extended.sCRM.Avoid),median(GB.Extended.sCRM.Avoid),median(KI.Extended.sCRM.Avoid[KI.Extended.sCRM.Avoid >0]),median(BH.Extended.sCRM.Avoid), median(CG.Extended.sCRM.Avoid),median(LT.Extended.sCRM.Avoid), median(ST.Extended.sCRM.Avoid),median(CT.Extended.sCRM.Avoid),median(SG.Extended.sCRM.Avoid),median(LG.Extended.sCRM.Avoid),median(AG.Extended.sCRM.Avoid),median(AT.Extended.sCRM.Avoid)),
							AvoidanceSD = c(sd(HG.Extended.sCRM.Avoid),sd(LB.Extended.sCRM.Avoid),sd(GB.Extended.sCRM.Avoid),sd(KI.Extended.sCRM.Avoid[KI.Extended.sCRM.Avoid >0]),sd(BH.Extended.sCRM.Avoid), sd(CG.Extended.sCRM.Avoid),sd(LT.Extended.sCRM.Avoid), sd(ST.Extended.sCRM.Avoid),sd(CT.Extended.sCRM.Avoid),sd(SG.Extended.sCRM.Avoid),sd(LG.Extended.sCRM.Avoid),sd(AG.Extended.sCRM.Avoid),sd(AT.Extended.sCRM.Avoid)),
							AvoidanceLCL = c(quantile(HG.Extended.sCRM.Avoid, 0.025),quantile(LB.Extended.sCRM.Avoid, 0.025),quantile(GB.Extended.sCRM.Avoid, 0.025),quantile(KI.Extended.sCRM.Avoid[KI.Extended.sCRM.Avoid >0], 0.025),quantile(BH.Extended.sCRM.Avoid, 0.025),quantile(CG.Extended.sCRM.Avoid, 0.025),quantile(LT.Extended.sCRM.Avoid, 0.025),quantile(ST.Extended.sCRM.Avoid, 0.025),quantile(CT.Extended.sCRM.Avoid, 0.025),quantile(SG.Extended.sCRM.Avoid, 0.025),quantile(LG.Extended.sCRM.Avoid, 0.025),quantile(AG.Extended.sCRM.Avoid, 0.025),quantile(AT.Extended.sCRM.Avoid, 0.025)), 
							AvoidanceUCL = c(quantile(HG.Extended.sCRM.Avoid, 0.975),quantile(LB.Extended.sCRM.Avoid, 0.975),quantile(GB.Extended.sCRM.Avoid, 0.975),quantile(KI.Extended.sCRM.Avoid[KI.Extended.sCRM.Avoid >0], 0.975),quantile(BH.Extended.sCRM.Avoid, 0.975),quantile(CG.Extended.sCRM.Avoid, 0.975),quantile(LT.Extended.sCRM.Avoid, 0.975),quantile(ST.Extended.sCRM.Avoid, 0.975),quantile(CT.Extended.sCRM.Avoid, 0.975),quantile(SG.Extended.sCRM.Avoid, 0.975),quantile(LG.Extended.sCRM.Avoid, 0.975),quantile(AG.Extended.sCRM.Avoid, 0.975),quantile(AT.Extended.sCRM.Avoid, 0.975)))

Extended.sCRM.Avoid.Table$ARs = paste(round(Extended.sCRM.Avoid.Table[,2],4), " (",round(Extended.sCRM.Avoid.Table[,3],4), "; ", round(Extended.sCRM.Avoid.Table[,4],4), " - ", round(Extended.sCRM.Avoid.Table[,5],4),")", sep = "") 


############################# Finally, the Christie & Urquhart model

##### Input parameters
## generic vectors
RelRad = seq(0.05, 1, 0.05)
RelChord = c(0.73, 0.79, 0.88, 0.96, 1, 0.98, 0.92, 0.85, 0.80, 0.75,0.7, 0.64,0.58,0.52,0.47,0.41,0.37,0.30,0.24, 0)
simpsons.weight = c(4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,1)
theta = seq(0,180,30)
c = cos(theta * pi /180)
s = sin(theta * pi/180)		
angle.head = seq(0,180,5) ### vector of bird headings in relation to headwind

if(nrow(coll.dat) >= detectCores()) cl <- makeCluster(detectCores())
if(nrow(coll.dat) < detectCores()) cl = makeCluster(niter)

registerDoParallel(cl)

clusterEvalQ(cl, library(truncnorm))

CuPColl = foreach(z = 1:nsim, .combine = "comb") %dopar% {

out = vector()

for(i in 1:nrow(coll.dat)){

simRotorSpeed = rnorm(1,coll.dat$rotor.speed[i],0.5)
simRotorPitch = rnorm(1,coll.dat$rotor.pitch[i],0.1)
simBodyLength = rnorm(1, coll.dat$BirdLength[i], coll.dat$BirdLengthSD[i])
simWingspan = rnorm(1, coll.dat$wingspan[i], coll.dat$wingspanSD[i])
simFlightSpeed = rtruncnorm(n = 1,a = 4, mean = coll.dat$flight.speed[i], sd = coll.dat$flightspeedSD[i])
sim.wind.speed = rtruncnorm(n = 1, a = 0, mean = 5, sd = 2)
period = 60/simRotorSpeed
pitch = simRotorPitch * pi/180 ## convert pitch to radians

if(sim.wind.speed > simFlightSpeed) sim.wind.speed = simFlightSpeed ### flight speed must be greater than wind speed 

out[i]= mean(mean.coll.risk (RelRad = RelRad, 
			RelChord = RelChord, 
			simpsons.weight = simpsons.weight, 
			theta = theta, 
			c = c,
			s = s,
			Radius = coll.dat$rotor.diameter[i]/2, 
			Blades = 3, 
			pitch = pitch, 
			chord.width = coll.dat$blade.width[i], 
			period = period, 
			length = simBodyLength, 
			wingspan = simWingspan, 
			speed = simFlightSpeed, 
			BF10G11 = 0, 
			min.angle = 0, 
			max.angle = 180, 
			wind.speed = sim.wind.speed))

			}
	list(out)

}


coll.no.avoid.CU = matrix(nrow = nrow(coll.dat), ncol = nsim)

for( i in 1:ncol(coll.no.avoid.CU)){  coll.no.avoid.CU[,i] = CuPColl[[1]][,i] * flux.dat[,i]  * coll.dat$X..flights.at.rotor.height }

HG.CU.Avoid = vector()
LB.CU.Avoid = vector()
GB.CU.Avoid = vector()
KI.CU.Avoid = vector()
BH.CU.Avoid = vector()
CG.CU.Avoid = vector()
LT.CU.Avoid = vector()
ST.CU.Avoid = vector()
CT.CU.Avoid = vector()
SG.CU.Avoid = vector()
LG.CU.Avoid = vector()
AG.CU.Avoid = vector()
AT.CU.Avoid = vector()

for (i in 1:ncol(coll.no.avoid.CU)){

HG.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "herring gull"])/sum(coll.no.avoid.CU[coll.dat$Species == "herring gull", i]))
LB.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "lesser black-backed gull"])/sum(coll.no.avoid.CU[coll.dat$Species == "lesser black-backed gull", i]))
GB.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "great black-backed gull"])/sum(coll.no.avoid.CU[coll.dat$Species == "great black-backed gull", i]))
KI.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "kittiwake"])/sum(coll.no.avoid.CU[coll.dat$Species == "kittiwake", i]))
BH.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "black-headed gull"])/sum(coll.no.avoid.CU[coll.dat$Species == "black-headed gull", i]))
CG.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common gull"])/sum(coll.no.avoid.CU[coll.dat$Species == "common gull", i]))
LT.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern"])/sum(coll.no.avoid.CU[coll.dat$Species == "little tern", i]))
ST.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern"])/sum(coll.no.avoid.CU[coll.dat$Species == "sandwich tern", i]))
CT.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern"])/sum(coll.no.avoid.CU[coll.dat$Species == "common tern", i]))
SG.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "small gull"])/sum(coll.no.avoid.CU[coll.dat$group == "small gull", i]))
LG.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "large gull"])/sum(coll.no.avoid.CU[coll.dat$group == "large gull", i]))
AG.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("large gull", "small gull", "gull")])/sum(coll.no.avoid.CU[coll.dat$group %in% c("gull","large gull", "small gull"), i]))
AT.CU.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern"])/sum(coll.no.avoid.CU[coll.dat$group == "tern", i]))

}




CU.Avoid.Table = data.frame(Species = c("Herring Gull", "Lesser Black-backed Gull", "Great Black-Backed Gull", "Kittiwake", "Black-headed Gull", "Common Gull", "Little Tern", "Sandwich Tern", "Common Tern", "Small Gull", "Large Gull", "All Gull", "All Tern"), 
							Avoidance = c(median(HG.CU.Avoid),median(LB.CU.Avoid),median(GB.CU.Avoid),median(KI.CU.Avoid[KI.CU.Avoid >0]),median(BH.CU.Avoid), median(CG.CU.Avoid),median(LT.CU.Avoid), median(ST.CU.Avoid),median(CT.CU.Avoid),median(SG.CU.Avoid),median(LG.CU.Avoid),median(AG.CU.Avoid),median(AT.CU.Avoid)),
							AvoidanceSD = c(sd(HG.CU.Avoid),sd(LB.CU.Avoid),sd(GB.CU.Avoid),sd(KI.CU.Avoid[KI.CU.Avoid >0]),sd(BH.CU.Avoid), sd(CG.CU.Avoid),sd(LT.CU.Avoid), sd(ST.CU.Avoid),sd(CT.CU.Avoid),sd(SG.CU.Avoid),sd(LG.CU.Avoid),sd(AG.CU.Avoid),sd(AT.CU.Avoid)),
							AvoidanceLCL = c(quantile(HG.CU.Avoid, 0.025),quantile(LB.CU.Avoid, 0.025),quantile(GB.CU.Avoid, 0.025),quantile(KI.CU.Avoid[KI.CU.Avoid >0], 0.025),quantile(BH.CU.Avoid, 0.025),quantile(CG.CU.Avoid, 0.025),quantile(LT.CU.Avoid, 0.025),quantile(ST.CU.Avoid, 0.025),quantile(CT.CU.Avoid, 0.025),quantile(SG.CU.Avoid, 0.025),quantile(LG.CU.Avoid, 0.025),quantile(AG.CU.Avoid, 0.025),quantile(AT.CU.Avoid, 0.025)), 
							AvoidanceUCL = c(quantile(HG.CU.Avoid, 0.975),quantile(LB.CU.Avoid, 0.975),quantile(GB.CU.Avoid, 0.975),quantile(KI.CU.Avoid[KI.CU.Avoid >0], 0.975),quantile(BH.CU.Avoid, 0.975),quantile(CG.CU.Avoid, 0.975),quantile(LT.CU.Avoid, 0.975),quantile(ST.CU.Avoid, 0.975),quantile(CT.CU.Avoid, 0.975),quantile(SG.CU.Avoid, 0.975),quantile(LG.CU.Avoid, 0.975),quantile(AG.CU.Avoid, 0.975),quantile(AT.CU.Avoid, 0.975)))


CU.Avoid.Table$ARs = paste(round(CU.Avoid.Table[,2],4), " (",round(CU.Avoid.Table[,3],4), "; ", round(CU.Avoid.Table[,4],4), " - ", round(CU.Avoid.Table[,5],4),")", sep = "") 



out.table = merge(Basic.Band.Avoid.Group.Table[,c(1,6)],Extended.Band.Avoid.Table[,c(1,6)], by = "Species")
names(out.table)[c(2,3)] = c("Basic Band Model", "Extended Band Model")
out.table = out.table[c(11,3,10,13, 7,1,2,9,6,4,12,8,5),]
out.table = cbind(out.table, Basic.sCRM.Avoid.Table[c(8,9,7,13,4,5,6,2,1,3,10,11,12),6], Extended.sCRM.Avoid.Table[c(8,9,7,13,4,5,6,2,1,3,10,11,12),6], CU.Avoid.Table[c(8,9,7,13,4,5,6,2,1,3,10,11,12),6], FCM.CollRate.Table [c(8,9,7,13,4,5,6,2,1,3,10,11,12),6])
names(out.table)[c(4,5,6,7)] = c("Basic sCRM", "Extended sCRM", "Christie & Urquhart Model", "Flux Collision Model")




###############################
############################### Look at implications of some of the assumptions (Basic Band model only)

### firstly recorded rather than estimated collisions - what happens if we don't correct for scavengers etc


Basic.Band.Unc.Spp.Avoid = vector()
Basic.Band.Unc.Spp.AvoidSD = vector()
Basic.Band.Unc.Spp.AvoidUpCI = vector()
Basic.Band.Unc.Spp.AvoidLowCI = vector()

for(i in 1:length(spp)) {

Basic.Band.Unc.Spp.Avoid[i] = 1-(sum(coll.dat$Recorded.Collisions[coll.dat$Species == spp[i]])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$Species == spp[i]]))
Basic.Band.Unc.Spp.AvoidSD[i] = std.dev(coll.dat$Recorded.Collisions[coll.dat$Species == spp[i]], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$Species == spp[i]])
Basic.Band.Unc.Spp.AvoidUpCI[i] = Basic.Band.Unc.Spp.Avoid[i] + (1.96 * Basic.Band.Unc.Spp.AvoidSD[i])
Basic.Band.Unc.Spp.AvoidLowCI[i] = Basic.Band.Unc.Spp.Avoid[i] - (1.96 * Basic.Band.Unc.Spp.AvoidSD[i])

}

Basic.Band.Unc.Avoid.Spp.Table = data.frame(Species = spp, Avoidance = Basic.Band.Unc.Spp.Avoid, AvoidanceSD = Basic.Band.Unc.Spp.AvoidSD, AvoidanceLCL = Basic.Band.Unc.Spp.AvoidLowCI, AvoidanceUCL = Basic.Band.Unc.Spp.AvoidUpCI)


group = unique(coll.dat$group)
Basic.Band.Unc.Group.Avoid = vector()
Basic.Band.Unc.Group.AvoidSD = vector()
Basic.Band.Unc.Group.AvoidUpCI = vector()
Basic.Band.Unc.Group.AvoidLowCI = vector()

for(i in 1:length(group)) {

Basic.Band.Unc.Group.Avoid[i] = 1-(sum(coll.dat$Recorded.Collisions[coll.dat$group == group[i]])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group == group[i]]))
Basic.Band.Unc.Group.AvoidSD[i] = std.dev(coll.dat$Recorded.Collisions[coll.dat$group == group[i]], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group == group[i]])
Basic.Band.Unc.Group.AvoidUpCI[i] = Basic.Band.Unc.Group.Avoid[i] + (1.96 * Basic.Band.Unc.Group.AvoidSD[i])
Basic.Band.Unc.Group.AvoidLowCI[i] = Basic.Band.Unc.Group.Avoid[i] - (1.96 * Basic.Band.Unc.Group.AvoidSD[i])

}


Basic.Band.Unc.Avoid.Group.Table = data.frame(Species = group, Avoidance = Basic.Band.Unc.Group.Avoid, AvoidanceSD = Basic.Band.Unc.Group.AvoidSD, AvoidanceLCL = Basic.Band.Unc.Group.AvoidLowCI, AvoidanceUCL = Basic.Band.Unc.Group.AvoidUpCI)


unc.ARs = rbind(Basic.Band.Unc.Avoid.Spp.Table[c(1:5, 9, 10,17,18),], Basic.Band.Unc.Avoid.Group.Table)

corr.ARs = rbind(Basic.Band.Avoid.Spp.Table[c(1:5, 9, 10,17,18),], Basic.Band.Avoid.Group.Table)
corr.ARs = unique(corr.ARs)
par(mar = c(4,0,0,0))

plot(unc.ARs[,2], seq(1,13,1), pch = 19, col = "red", bty = "n", xaxt = "n", xlab = "Avoidance Rate", yaxt = "n", ylab = "", xlim = c(0.85, 1))
points(corr.ARs[,2], seq(1,13,1), pch = 19, col = "blue")
text (0.85, seq(1,13,1), corr.ARs[,1], pos = 4)
axis(1, at = c(0.9, 0.95, 1))

for(i in 1:nrow(unc.ARs)) { arrows (unc.ARs[i,2], i, corr.ARs[i,2], i,  length = 0.05, lwd = 2)} 

###### Sensitivity to other parameters

coll.dat$PcollFlightSpeed  = 0
coll.dat$PcollWingSpan = 0
coll.dat$PcollBodyLength = 0
coll.dat$PcollRotorSpeed = 0
coll.dat$PcollRotorPitch = 0

for(i in 1:nrow(coll.dat)){

coll.dat$PcollFlightSpeed[i] = pcoll(k =1, no.blades = 3, max.chord = coll.dat$blade.width[i], pitch = coll.dat$rotor.pitch[i], rotor.rad = coll.dat$rotor.diameter[i]/2, rotor.period = 60/coll.dat$rotor.speed[i],wingspan = coll.dat$wingspan[i], bird.length = coll.dat$BirdLength[i], aspect = coll.dat$BirdLength[i]/coll.dat$wingspan[i], flap.glide = 0, bird.speed = coll.dat$flight.speed[i]*1.1,prop.flights.upwind = 0.5,radius,chord)
coll.dat$PcollWingSpan[i] = pcoll(k =1, no.blades = 3, max.chord = coll.dat$blade.width[i], pitch = coll.dat$rotor.pitch[i], rotor.rad = coll.dat$rotor.diameter[i]/2, rotor.period = 60/coll.dat$rotor.speed[i],wingspan = coll.dat$wingspan[i]*1.1, bird.length = coll.dat$BirdLength[i], aspect = coll.dat$BirdLength[i]/(coll.dat$wingspan[i]*1.1), flap.glide = 0, bird.speed = coll.dat$flight.speed[i],prop.flights.upwind = 0.5,radius,chord)
coll.dat$PcollBodyLength[i] = pcoll(k =1, no.blades = 3, max.chord = coll.dat$blade.width[i], pitch = coll.dat$rotor.pitch[i], rotor.rad = coll.dat$rotor.diameter[i]/2, rotor.period = 60/coll.dat$rotor.speed[i],wingspan = coll.dat$wingspan[i], bird.length = coll.dat$BirdLength[i]*1.1, aspect = (coll.dat$BirdLength[i]*1.1)/coll.dat$wingspan[i], flap.glide = 0, bird.speed = coll.dat$flight.speed[i],prop.flights.upwind = 0.5,radius,chord)
coll.dat$PcollRotorSpeed[i] = pcoll(k =1, no.blades = 3, max.chord = coll.dat$blade.width[i], pitch = coll.dat$rotor.pitch[i], rotor.rad = coll.dat$rotor.diameter[i]/2, rotor.period = 60/coll.dat$rotor.speed[i]*1.1,wingspan = coll.dat$wingspan[i], bird.length = coll.dat$BirdLength[i], aspect = coll.dat$BirdLength[i]/coll.dat$wingspan[i], flap.glide = 0, bird.speed = coll.dat$flight.speed[i],prop.flights.upwind = 0.5,radius,chord)
coll.dat$PcollRotorPitch[i] = pcoll(k =1, no.blades = 3, max.chord = coll.dat$blade.width[i], pitch = coll.dat$rotor.pitch[i]*1.1, rotor.rad = coll.dat$rotor.diameter[i]/2, rotor.period = 60/coll.dat$rotor.speed[i],wingspan = coll.dat$wingspan[i], bird.length = coll.dat$BirdLength[i], aspect = coll.dat$BirdLength[i]/coll.dat$wingspan[i], flap.glide = 0, bird.speed = coll.dat$flight.speed[i],prop.flights.upwind = 0.5,radius,chord)

}

coll.dat$flux.CRHPCH = coll.dat$total.flux * (coll.dat$X..flights.at.rotor.height*1.1) * (coll.dat$total.rotor.frontal.area/coll.dat$area.survey.window)
coll.dat$flux.CRHFLUX = (coll.dat$total.flux*1.1) * coll.dat$X..flights.at.rotor.height*1.1 * (coll.dat$total.rotor.frontal.area/coll.dat$area.survey.window)

coll.dat$PredColl.NoAvoid.Opt1FlightSpeed = coll.dat$PcollFlightSpeed *  coll.dat$flux.CRH
coll.dat$PredColl.NoAvoid.Opt1WingSpan = coll.dat$PcollWingSpan *  coll.dat$flux.CRH
coll.dat$PredColl.NoAvoid.Opt1BodyLength = coll.dat$PcollBodyLength *  coll.dat$flux.CRH
coll.dat$PredColl.NoAvoid.Opt1RotorSpeed = coll.dat$PcollRotorSpeed *  coll.dat$flux.CRH
coll.dat$PredColl.NoAvoid.Opt1RotorPitch = coll.dat$PcollRotorPitch *  coll.dat$flux.CRH
coll.dat$PredColl.NoAvoid.Opt1UpPCH = coll.dat$Pcoll *  coll.dat$flux.CRHPCH
coll.dat$PredColl.NoAvoid.Opt1FLUX = coll.dat$PcollUp *  coll.dat$flux.CRHFLUX


Basic.Band.Group.Avoid = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group %in% c("small gull", "large gull", "gull")]))
Basic.Band.Group.AvoidFlightSpeed = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt1FlightSpeed[coll.dat$group %in% c("small gull", "large gull", "gull")]))
Basic.Band.Group.AvoidWingSpan = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt1WingSpan[coll.dat$group %in% c("small gull", "large gull", "gull")]))
Basic.Band.Group.AvoidBodyLength = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt1BodyLength[coll.dat$group %in% c("small gull", "large gull", "gull")]))
Basic.Band.Group.AvoidRotorSpeed = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt1RotorSpeed[coll.dat$group %in% c("small gull", "large gull", "gull")]))
Basic.Band.Group.AvoidRotorPitch = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt1RotorPitch[coll.dat$group %in% c("small gull", "large gull", "gull")]))
Basic.Band.Group.AvoidPCH = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt1UpPCH[coll.dat$group %in% c("small gull", "large gull", "gull")]))
Basic.Band.Group.AvoidFlux = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt1FLUX[coll.dat$group %in% c("small gull", "large gull", "gull")]))

xmin = min(floor(c(Basic.Band.Group.Avoid,Basic.Band.Group.AvoidPCH,Basic.Band.Group.AvoidFlux,Basic.Band.Group.AvoidFlightSpeed,Basic.Band.Group.AvoidWingSpan,Basic.Band.Group.AvoidBodyLength,Basic.Band.Group.AvoidRotorSpeed,Basic.Band.Group.AvoidRotorPitch)*100))/100
xmax = min(ceiling(c(Basic.Band.Group.Avoid,Basic.Band.Group.AvoidPCH,Basic.Band.Group.AvoidFlux,Basic.Band.Group.AvoidFlightSpeed,Basic.Band.Group.AvoidWingSpan,Basic.Band.Group.AvoidBodyLength,Basic.Band.Group.AvoidRotorSpeed,Basic.Band.Group.AvoidRotorPitch)*100))/100

plot(c(c(Basic.Band.Group.AvoidPCH,Basic.Band.Group.AvoidFlux,Basic.Band.Group.AvoidFlightSpeed,Basic.Band.Group.AvoidWingSpan,Basic.Band.Group.AvoidBodyLength,Basic.Band.Group.AvoidRotorSpeed,Basic.Band.Group.AvoidRotorPitch)), seq(1,7,1), pch = 19, col = "blue", yaxt = "n", xlim = c(0.984,0.990), bty = "n", xlab = "Avoidance Rate", ylab = "")
points(rep(Basic.Band.Group.Avoid,7), seq(1,7,1), pch = 19, col = "red")

text(0.984, seq(1,8,1), c("Height", "Flux", "Flight Speed", "Wingspan", "Bird Length", "Rotor Speed", "Rotor Pitch"), pos = 4)


arrows (Basic.Band.Group.Avoid, 1, Basic.Band.Group.AvoidPCH, 1, length = 0.05, lwd = 2)
arrows (Basic.Band.Group.Avoid, 2, Basic.Band.Group.AvoidFlux, 2, length = 0.05, lwd = 2)
arrows (Basic.Band.Group.Avoid, 3, Basic.Band.Group.AvoidFlightSpeed, 3, length = 0.05, lwd = 2)
arrows (Basic.Band.Group.Avoid, 4, Basic.Band.Group.AvoidWingSpan, 4, length = 0.05, lwd = 2)
arrows (Basic.Band.Group.Avoid, 5, Basic.Band.Group.AvoidBodyLength, 5, length = 0.05, lwd = 2)
arrows (Basic.Band.Group.Avoid, 6, Basic.Band.Group.AvoidRotorSpeed, 6, length = 0.05, lwd = 2)
arrows (Basic.Band.Group.Avoid, 7, Basic.Band.Group.AvoidRotorPitch, 7, length = 0.05, lwd = 2)

################# Sandwich Terns/Terns without Zeebrugge


Basic.Band.Spp.AvoidST = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"]))
Basic.Band.Spp.AvoidSDST = std.dev(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])
Basic.Band.Spp.AvoidUpCIST = Basic.Band.Spp.AvoidST + (1.96 * Basic.Band.Spp.AvoidSDST)
Basic.Band.Spp.AvoidLowCIST = Basic.Band.Spp.AvoidST - (1.96 * Basic.Band.Spp.AvoidSDST)

Extended.Band.Spp.AvoidST = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.dat$PredColl.NoAvoid.Opt3[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"]))
Extended.Band.Spp.AvoidSDST = std.dev(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"], coll.dat$PredColl.NoAvoid.Opt3[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])
Extended.Band.Spp.AvoidUpCIST = Extended.Band.Spp.AvoidST + (1.96 * Extended.Band.Spp.AvoidSDST)
Extended.Band.Spp.AvoidLowCIST = Extended.Band.Spp.AvoidST - (1.96 * Extended.Band.Spp.AvoidSDST)

Basic.Band.Spp.AvoidCT = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"]))
Basic.Band.Spp.AvoidSDCT = std.dev(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])
Basic.Band.Spp.AvoidUpCICT = Basic.Band.Spp.AvoidCT + (1.96 * Basic.Band.Spp.AvoidSDCT)
Basic.Band.Spp.AvoidLowCICT = Basic.Band.Spp.AvoidCT - (1.96 * Basic.Band.Spp.AvoidSDCT)

Extended.Band.Spp.AvoidCT = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.dat$PredColl.NoAvoid.Opt3[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"]))
Extended.Band.Spp.AvoidSDCT = std.dev(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"], coll.dat$PredColl.NoAvoid.Opt3[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])
Extended.Band.Spp.AvoidUpCICT = Extended.Band.Spp.AvoidCT + (1.96 * Extended.Band.Spp.AvoidSDCT)
Extended.Band.Spp.AvoidLowCICT = Extended.Band.Spp.AvoidCT - (1.96 * Extended.Band.Spp.AvoidSDCT)

Basic.Band.Spp.AvoidAT = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"]))
Basic.Band.Spp.AvoidSDAT = std.dev(coll.dat$Estimated.Collisions[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])
Basic.Band.Spp.AvoidUpCIAT = Basic.Band.Spp.AvoidAT + (1.96 * Basic.Band.Spp.AvoidSDAT)
Basic.Band.Spp.AvoidLowCIAT = Basic.Band.Spp.AvoidAT - (1.96 * Basic.Band.Spp.AvoidSDAT)

Extended.Band.Spp.AvoidAT = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.dat$PredColl.NoAvoid.Opt3[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"]))
Extended.Band.Spp.AvoidSDAT = std.dev(coll.dat$Estimated.Collisions[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"], coll.dat$PredColl.NoAvoid.Opt3[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])
Extended.Band.Spp.AvoidUpCIAT = Extended.Band.Spp.AvoidAT + (1.96 * Extended.Band.Spp.AvoidSDAT)
Extended.Band.Spp.AvoidLowCIAT = Extended.Band.Spp.AvoidAT - (1.96 * Extended.Band.Spp.AvoidSDAT)


LT.Basic.sCRM.AvoidNoZee = vector()
ST.Basic.sCRM.AvoidNoZee = vector()
CT.Basic.sCRM.AvoidNoZee = vector()
AT.Basic.sCRM.AvoidNoZee = vector()

LT.Extended.sCRM.AvoidNoZee = vector()
ST.Extended.sCRM.AvoidNoZee = vector()
CT.Extended.sCRM.AvoidNoZee = vector()
AT.Extended.sCRM.AvoidNoZee = vector()

LT.CU.AvoidNoZee = vector()
ST.CU.AvoidNoZee = vector()
CT.CU.AvoidNoZee = vector()
AT.CU.AvoidNoZee = vector()

LT.FCM.CollRateNoZee = vector()
ST.FCM.CollRateNoZee = vector()
CT.FCM.CollRateNoZee = vector()
AT.FCM.CollRateNoZee = vector()

for (i in 1:ncol(coll.no.avoid)){
#basic
LT.Basic.sCRM.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid[coll.dat$Species == "little tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
ST.Basic.sCRM.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
CT.Basic.sCRM.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
AT.Basic.sCRM.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))

#extended
LT.Extended.sCRM.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "little tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
ST.Extended.sCRM.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
CT.Extended.sCRM.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
AT.Extended.sCRM.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid.Opt3[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))

#Christie & Urquhart
LT.CU.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid.CU[coll.dat$Species == "little tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
ST.CU.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid.CU[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
CT.CU.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid.CU[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
AT.CU.AvoidNoZee[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(coll.no.avoid.CU[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))

#FCM
LT.FCM.CollRateNoZee[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(flux.dat[coll.dat$Species == "little tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
ST.FCM.CollRateNoZee[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(flux.dat[coll.dat$Species == "sandwich tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
CT.FCM.CollRateNoZee[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(flux.dat[coll.dat$Species == "common tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))
AT.FCM.CollRateNoZee[i] = (sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge"])/sum(flux.dat[coll.dat$group == "tern" & coll.dat$Wind.Farm.2 != "zeebrugge", i]))

}

median(ST.Basic.sCRM.AvoidNoZee)
sd(ST.Basic.sCRM.AvoidNoZee)
quantile(ST.Basic.sCRM.AvoidNoZee, 0.025)
quantile(ST.Basic.sCRM.AvoidNoZee, 0.975)

median(ST.Extended.sCRM.AvoidNoZee)
sd(ST.Extended.sCRM.AvoidNoZee)
quantile(ST.Extended.sCRM.AvoidNoZee, 0.025)
quantile(ST.Extended.sCRM.AvoidNoZee, 0.975)

median(ST.CU.AvoidNoZee)
sd(ST.CU.AvoidNoZee)
quantile(ST.CU.AvoidNoZee, 0.025)
quantile(ST.CU.AvoidNoZee, 0.975)

median(ST.FCM.CollRateNoZee)
sd(ST.FCM.CollRateNoZee)
quantile(ST.FCM.CollRateNoZee, 0.025)
quantile(ST.FCM.CollRateNoZee, 0.975)


median(CT.Basic.sCRM.AvoidNoZee)
sd(CT.Basic.sCRM.AvoidNoZee)
quantile(CT.Basic.sCRM.AvoidNoZee, 0.025)
quantile(CT.Basic.sCRM.AvoidNoZee, 0.975)

median(CT.Extended.sCRM.AvoidNoZee)
sd(CT.Extended.sCRM.AvoidNoZee)
quantile(CT.Extended.sCRM.AvoidNoZee, 0.025)
quantile(CT.Extended.sCRM.AvoidNoZee, 0.975)

median(CT.CU.AvoidNoZee)
sd(CT.CU.AvoidNoZee)
quantile(CT.CU.AvoidNoZee, 0.025)
quantile(CT.CU.AvoidNoZee, 0.975)

median(CT.FCM.CollRateNoZee)
sd(CT.FCM.CollRateNoZee)
quantile(CT.FCM.CollRateNoZee, 0.025)
quantile(CT.FCM.CollRateNoZee, 0.975)

median(LT.Basic.sCRM.AvoidNoZee)
sd(LT.Basic.sCRM.AvoidNoZee)
quantile(LT.Basic.sCRM.AvoidNoZee, 0.025)
quantile(LT.Basic.sCRM.AvoidNoZee, 0.975)

median(LT.Extended.sCRM.AvoidNoZee)
sd(LT.Extended.sCRM.AvoidNoZee)
quantile(LT.Extended.sCRM.AvoidNoZee, 0.025)
quantile(LT.Extended.sCRM.AvoidNoZee, 0.975)

median(LT.CU.AvoidNoZee)
sd(LT.CU.AvoidNoZee)
quantile(LT.CU.AvoidNoZee, 0.025)
quantile(LT.CU.AvoidNoZee, 0.975)

median(LT.FCM.CollRateNoZee)
sd(LT.FCM.CollRateNoZee)
quantile(LT.FCM.CollRateNoZee, 0.025)
quantile(LT.FCM.CollRateNoZee, 0.975)

median(AT.Basic.sCRM.AvoidNoZee)
sd(AT.Basic.sCRM.AvoidNoZee)
quantile(AT.Basic.sCRM.AvoidNoZee, 0.025)
quantile(AT.Basic.sCRM.AvoidNoZee, 0.975)

median(AT.Extended.sCRM.AvoidNoZee)
sd(AT.Extended.sCRM.AvoidNoZee)
quantile(AT.Extended.sCRM.AvoidNoZee, 0.025)
quantile(AT.Extended.sCRM.AvoidNoZee, 0.975)

median(AT.CU.AvoidNoZee)
sd(AT.CU.AvoidNoZee)
quantile(AT.CU.AvoidNoZee, 0.025)
quantile(AT.CU.AvoidNoZee, 0.975)

median(AT.FCM.CollRateNoZee)
sd(AT.FCM.CollRateNoZee)
quantile(AT.FCM.CollRateNoZee, 0.025)
quantile(AT.FCM.CollRateNoZee, 0.975)

#################### data for species accounts

unique(coll.dat$Wind.Farm.2[coll.dat$Species == "sandwich tern"])
 min(coll.dat$passage.rate[coll.dat$Species == "sandwich tern"])
 max(coll.dat$passage.rate[coll.dat$Species == "sandwich tern"])
sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern"])

unique(coll.dat$Wind.Farm.2[coll.dat$Species == "common tern"])
 min(coll.dat$passage.rate[coll.dat$Species == "common tern"])
 max(coll.dat$passage.rate[coll.dat$Species == "common tern"])
sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern"])

unique(coll.dat$Wind.Farm.2[coll.dat$Species == "little tern"])
 min(coll.dat$passage.rate[coll.dat$Species == "little tern"])
 max(coll.dat$passage.rate[coll.dat$Species == "little tern"])
sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern"])

unique(coll.dat$Wind.Farm.2[coll.dat$group == "tern"])
sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern"])

aggregate(passage.rate~Wind.Farm.2 + Year, data = coll.dat[coll.dat$group == "tern",], FUN = sum)

unique(coll.dat$Wind.Farm.2[coll.dat$Species == "kittiwake"])
 min(coll.dat$passage.rate[coll.dat$Species == "kittiwake"])
 max(coll.dat$passage.rate[coll.dat$Species == "kittiwake"])
sum(coll.dat$Estimated.Collisions[coll.dat$Species == "kittiwake"])

unique(coll.dat$Wind.Farm.2[coll.dat$Species == "black-headed gull"])
 min(coll.dat$passage.rate[coll.dat$Species == "black-headed gull"])
 max(coll.dat$passage.rate[coll.dat$Species == "black-headed gull"])
sum(coll.dat$Estimated.Collisions[coll.dat$Species == "black-headed gull"])

unique(coll.dat$Wind.Farm.2[coll.dat$Species == "common gull"])
 min(coll.dat$passage.rate[coll.dat$Species == "common gull"])
 max(coll.dat$passage.rate[coll.dat$Species == "common gull"])
sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common gull"])


unique(coll.dat$Wind.Farm.2[coll.dat$Species == "lesser black-backed gull"])
 min(coll.dat$passage.rate[coll.dat$Species == "lesser black-backed gull"])
 max(coll.dat$passage.rate[coll.dat$Species == "lesser black-backed gull"])
sum(coll.dat$Estimated.Collisions[coll.dat$Species == "lesser black-backed gull"])

unique(coll.dat$Wind.Farm.2[coll.dat$Species == "herring gull"])
 min(coll.dat$passage.rate[coll.dat$Species == "herring gull"])
 max(coll.dat$passage.rate[coll.dat$Species == "herring gull"])
sum(coll.dat$Estimated.Collisions[coll.dat$Species == "herring gull"])


unique(coll.dat$Wind.Farm.2[coll.dat$Species == "great black-backed gull"])
 min(coll.dat$passage.rate[coll.dat$Species == "great black-backed gull"])
 max(coll.dat$passage.rate[coll.dat$Species == "great black-backed gull"])
sum(coll.dat$Estimated.Collisions[coll.dat$Species == "great black-backed gull"])

unique(coll.dat$Wind.Farm.2[coll.dat$group == "small gull"])
sum(coll.dat$Estimated.Collisions[coll.dat$group == "small gull"])

aggregate(passage.rate~Wind.Farm.2 + Year, data = coll.dat[coll.dat$group == "small gull",], FUN = sum)

unique(coll.dat$Wind.Farm.2[coll.dat$group %in% c("gull", "small gull", "large gull")])
sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("gull", "small gull", "large gull")])

aggregate(passage.rate~Wind.Farm.2 + Year, data = coll.dat[coll.dat$group %in% c("gull", "small gull", "large gull"),], FUN = sum)
