# Packages to uniformly distribute lat/lon on a sphere and plot map
require(geosphere)
require(rworldmap)


# Steffan-Boltzmann Law
SBlaw <- function(I, alpha = 0.3, epsilon = 1, sigma = 5.670373e-8){
	(I*(1-alpha)/(epsilon*sigma))^.25
}

# Calculate intensity of sunlight at each lat/lon ##
#  -The light is brightest at lat = 0, lon =0 (max 1367 W/m^2)
#   We need to convert lat/lon to radians for R's cos function
#   Irrad cannot be negative, so a lower bound is set at zero
set.seed(1234)
Imax   = 1367
Npts   = 100
LonLat = randomCoordinates(Npts)*pi/180
Irrad  = pmax(0, Imax*cos(LonLat[, "lon"])*cos(LonLat[, "lat"]))
res    = as.data.frame(cbind(LonLat, Irrad, Temp = SBlaw(Irrad)))

# Mean Temperature
mean(SBlaw(Irrad))

# Equilibrium Temperature
SBlaw(mean(Irrad))