Nothing
# solar power functions
I0.orbit <- function(leap=FALSE,plot=TRUE){
if(leap==TRUE){ny=366; npx=4} else {ny=365; npx=3}
np=182+npx # July 4 perihelion
na=1+2 # January 3 aphelion
# Solar constant
SC = 1.377 # kW/m2
# ET solar radiation I0
n <- seq(1,ny)
I0 <- SC*(1+0.033*cos((2*pi/ny)*(n-np))) # kW/m2
ylabel <- "I0 ET radiation (kW/m2)"
if(plot==TRUE){
y <- I0;x <- n
ymin <- min(y); ymax <- max(y)
yr <- (ymax - ymin)
shift <- 0.02*yr
yseq <- round(seq(ymin,ymax,0.1*yr),2)
ypad <- c(ymin-shift, ymax+shift)
plot(x,y,type="l",xlab="Number of Day",ylab=ylabel,lwd=2,ylim=ypad)
abline(h=SC,lty=2);text(20,SC+shift,"Average",cex=0.7)
abline(v=np,lty=2);text(np+2,ymin-shift,"Perihelion",cex=0.7)
abline(v=na,lty=2);text(na+2,ymin-shift,"Aphelion",cex=0.7)
}
return(list(I0=I0))
}
days.mo <- function(day,leap=FALSE){
# day for the month
if (leap==FALSE)daymo <- c(31,28,31,30,31,30,31,31,30,31,30,31)
else daymo <- c(31,29,31,30,31,30,31,31,30,31,30,31)
day.mo <-array()
day.mo[1] <- day
for(i in 2:12) day.mo[i] <- day.mo[i-1]+ daymo[i-1]
mid.mo = c(15,46,74,105,135,166,196,227,258,288,319,349)
return(list(day.mo=day.mo))
}
declination <- function(leap=FALSE,plot=TRUE){
# declination
if(leap==TRUE){ny=366; ne=81} else {ny=365; ne=80}
n=seq(1,ny)
dec <- 23.45*sin((2*pi/ny)*(n-ne))
if(plot==TRUE){
plot(n,dec,type="l",xlab="Day number",ylab=expression("Declination [ "*degree*"]"),lwd=2)
abline(h=0,lty=2);text(20,0,"Average",cex=0.7)
abline(v=ne,lty=2);text(90,-20,"March Equinox",cex=0.7)
abline(v=ne+183,lty=2);text(260,-20,"September Equinox",cex=0.7)
}
return(list(ne=ne,dec=dec))
}
sun.elev <- function(xdec,lat,plot=TRUE){
elev <- xdec$dec + 90 - lat
me <- round(mean(elev),2)
mine <- min(elev); maxe <- max(elev)
shift <- 0.02*(maxe - mine)
ypad <- c(mine, 1.05*maxe)
n <- seq(1,length(elev),1)
if(plot==TRUE){
plot(n,elev,type="l",xlab="Day number",ylab=expression("Sun Angle Elev [ "*degree*"]"),lwd=2,ylim=ypad)
abline(h=me,lty=2);text(15,me+shift,as.expression(bquote("Average= "*.(me)*degree)),cex=0.7)
abline(v=xdec$ne,lty=2);text(xdec$ne,mine,"Equinox",cex=0.7)
abline(v=xdec$ne+183,lty=2);text(xdec$ne+183,mine,"Equinox",cex=0.7)
mtext(side=3,line=-1,as.expression(bquote("Lat= "*.(lat)*degree)),cex=0.7)
}
return(list(ne=xdec$ne,elev=elev,lat=lat))
}
read.tau <- function(file){
loc<- read.csv(file,header=FALSE,sep=",",nrows=1)
lat.long.elev <- read.table(file,sep=",",skip=1,header=TRUE,nrows=1)
tau <- read.table(file,skip=3,header=TRUE,sep=",")
return(list(loc=loc, lat.long.elev=lat.long.elev, tau=tau))
}
beam.diffuse <- function(dat,plot=TRUE){
kab <- c(1.454,-0.406,-0.268,0.021)
kad <- c(0.507,0.205,-0.08,-0.19)
lat <- as.numeric(dat$lat.long.elev[1])
tau <- dat$tau
y <- sun.elev(declination(plot=FALSE),lat=lat)
z <- I0.orbit(plot=FALSE)
ab <- array(); ad <- array()
Ib <- array(); Id <- array()
day21.mo <- days.mo(21)$day.mo
elev <- y$elev[day21.mo]
I0 <- z$I0[day21.mo]
air.mass0 <- 1/sin(elev*2*pi/360)
air.mass1 <- sqrt((708*sin(elev*2*pi/360))^2+1417)-708*sin(elev*2*pi/360)
#air.mass2 <- 1/(sin(elev*2*pi/360)+0.50572*(6.07995+elev)^-1.6364)
#matplot(day21.mo,cbind(air.mass0,air.mass1,air.mass2),type='b')
air.mass <- air.mass1
for(i in 1:12){
ab[i] <- kab[1]+kab[2]*tau[i,2]+kab[3]*tau[i,3]+kab[4]*tau[i,2]*tau[i,3]
ad[i] <- kad[1]+kad[2]*tau[i,3]+kad[3]*tau[i,2]+kad[4]*tau[i,2]*tau[i,3]
Ib[i] <- I0[i]*exp(-tau[i,2]*air.mass[i]^ab[i])
Id[i] <- I0[i]*exp(-tau[i,3]*air.mass[i]^ad[i])
}
Id.Ib <- round(100*(Id/Ib),2)
if(plot==TRUE){
panels(7,7,3,2,pty="m")
plot(day21.mo,tau[,2],xlab="Day Number",ylab="Beam Pseudo Optical Depth",type='b')
plot(day21.mo,I0,xlab="Day Number",ylab="I0 (kW/m2)",type='b')
plot(day21.mo,tau[,3],xlab="Day Number",ylab="Diffuse Pseudo Optical Depth",type='b')
plot(day21.mo,Ib,xlab="Day Number",ylab="Direct Beam Ibn (kW/m2)",type='b')
plot(day21.mo,air.mass,xlab="Day Number",ylab="Air Mass",type='b')
plot(day21.mo,Id,xlab="Day Number",ylab="Diffuse Id (kW/m2)",type='b')
}
return(list(day21.mo=day21.mo,tau=tau,air.mass=air.mass,lat=lat,Ib=Ib,Id=Id,Id.Ib=Id.Ib))
}
I0.blackbody <- function(T.sun,wl.nm,plot=TRUE){
# solar spectrum from Plancks law
lambda <- wl.nm
lambda.m <- lambda*10^-9 # m
c.light <- 2.99*10^8 # m/s
nu <- c.light/lambda.m
h.planck <- 6.62606957*10^-34 # Js
k.boltz <- 1.3806488*10^-23 # J/K
radius.earth <- 1.496*10^8; radius.sun <- 6.95*10^5
f.rad <- (radius.sun/radius.earth)^2
hnu.kT <- (h.planck*nu)/(k.boltz*T.sun)
mult <- 2*pi*h.planck*c.light^2/lambda.m^5
I.sun.m <- f.rad* mult /(exp(hnu.kT)-1)
I.sun.nm <- I.sun.m *10^-9
I0 <- sum(I.sun.nm)
ymax <- max(I.sun.nm)
if(plot==TRUE){
plot(lambda,I.sun.nm,type="l",
xlab="Wavelength (nm)",ylab="Spectral irradiance (W/m2/nm)",
xaxs="r",yaxs="r",ylim=c(0,ymax),lwd=2,cex.axis=0.8)
abline(v=c(400,700,1400),lty=2)
text(c(300,550,900,1500),rep(0.05,4),c("UV","Visible","NIR","IR"),cex=0.8)
legend('topright', legend=c(paste("Blackbody ",T.sun,"K")),cex=0.8)
}
return(list(lambda=lambda,I.sun.nm=I.sun.nm,I0=I0))
}
spectral <- function(X,label,wl.lim.nm,T.sun,plot.surf=FALSE){
lambda <- X[,1]
I0.bb.list <- I0.blackbody(T.sun=T.sun,wl.nm=lambda,plot=FALSE)
I0.bb <- I0.bb.list$I.sun.nm
I0 <- X[,2]
I <- X[,3]
ymax <- max(X[,2:4])
labelI0 <- paste(label, "I0"); labelI <- paste(label, "I");
label.bb <- paste("Blackbody ",T.sun,"K")
matplot(lambda,cbind(I0.bb,I0),type="l", col=1, lty=c(1,1),
xlab="Wavelength (nm)",ylab="Spectral irradiance (W/m2/nm)",
xaxs="r",yaxs="r",ylim=c(0,ymax),xlim=wl.lim.nm,lwd=c(2,0.5),cex.axis=0.8)
legend('topright', legend=c(label.bb,labelI0),
col=1, lty=c(1,1),lwd=c(2,0.5),cex=0.8)
abline(v=c(400,700,1400),lty=2)
text(c(300,550,900,1500),rep(0.05,4),c("UV","Visible","NIR","IR"),cex=0.8)
if(plot.surf==TRUE){
matplot(lambda,cbind(I0.bb,I0,I),type="l", col=c('gray','gray','black'), lty=1,
xlab="Wavelength (nm)",ylab="Spectral irradiance (W/m2/nm)",
xaxs="r",yaxs="r",ylim=c(0,ymax),xlim=wl.lim.nm,lwd=c(2,0.5,0.5),cex.axis=0.8)
legend('topright', legend=c(label.bb,labelI0,labelI),col=c('gray','gray','black'),
lty=1,lwd=c(2,0.5,0.5),cex=0.8)
abline(v=c(400,700,1400),lty=2)
text(c(300,550,900,1500),rep(0.05,4),c("UV","Visible","NIR","IR"),cex=0.8)
}
}
useful.waste <- function(I0.bb){
I0 <- I0.bb
I.sun.nm <- I0$I.sun.nm; lambda <- I0$lambda
lambda.m <- lambda*10^-9
c.light <- 2.99*10^8 # m/s
h.planck <- 6.62606957*10^-34 # Js
k.boltz <- 1.3806488*10^-23 # J/K
radius.earth <- 1.496*10^8; radius.sun <- 6.95*10^5
f.rad <- (radius.sun/radius.earth)^2
y.1100 <- I.sun.nm[which(lambda==1100)]
x.y.1100 <- lambda[min(which(round(I.sun.nm,2)== round(y.1100,2)))]
lines(c(1100,1100),c(0,y.1100),lty=2)
#lines(c(x.y.1100,1100),c(y.1100,y.1100),lty=2)
# eV/J
eV.J <- 1.6* 10^(-19)
# Energy of photon
Ep <- (h.planck*c.light/lambda.m)/eV.J
plot(lambda,Ep, type="l",xlab="Wavelength(nm)",ylab="Photon energy (eV)",lwd=2)
abline(h=1.12,lty=2); text(400,1.3,"Eg=1.12 eV",cex=0.7)
abline(v=1110,lty=2); text(1300,2,"1110 nm",cex=0.7)
arrows(1110,3,600,3,length=0.1); text(800,3.2,"Enough energy",cex=0.7)
arrows(1110,0.4,2000,0.4,length=0.1); text(1500,0.6,"Not Enough energy",cex=0.7)
y <- Ep - 1.12
for(i in 1: length(y)) if(y[i] <0) y[i] <- 0
waste <- sum(y)/sum(Ep)
}
sun.path <- function(lat,nday,plot=TRUE){
# declination
dec <- 23.45*sin((2*pi/365)*(nday-81))
# hr angle and trig
hr.noon <- seq(-6,+6,1)
hr.angle <- 15*hr.noon
sin.H <- sin((pi/180)*hr.angle)
cos.H <- cos((pi/180)*hr.angle)
# elevation and angle
sin.elev <- cos((pi/180)*lat)*cos((pi/180)*dec)*cos.H +
sin((pi/180)*lat)*sin((pi/180)*dec)
elev <- asin(sin.elev)*180/pi
# azimuth
sin.azi <- cos((pi/180)*dec) * sin.H/cos((pi/180)*elev)
# value for testing azimuth
test.tan <- tan(dec*pi/180)/tan(lat*pi/180)
kn <- which(hr.noon==0)
nn <- length(cos.H)
azi <- array()
for(k in 1:nn){
if(cos.H[k] >= test.tan){
azi[k] <- asin(sin.azi[k])*180/pi
}
else{
if(k <= kn) {corr = -180} else {corr = 180}
azi[k] <- corr - asin(sin.azi[k])*180/pi
}
}
if(plot==TRUE){
wd=7; ht=7
panels(wd,ht,2,1,pty="m")
plot(hr.noon,elev,type="b",xlab="Hour before noon",ylab=expression("Sun angle ["*degree*"]"),lwd=2)
text(0,10,as.expression(bquote("Lat= "*.(lat)*degree*" Day= "*.(nday))),cex.main=0.9)
plot(hr.noon,azi,type="b",xlab="Hour before noon",ylab=expression("Azimuth ["*degree*"]"),lwd=2)
text(-4,10,as.expression(bquote("Lat= "*.(lat)*degree*" Day= "*.(nday))),cex.main=0.9)
plot(azi,elev,type="b",xlab=expression("Azimuth ["*degree*"]"),ylab=expression("Sun Elevation ["*degree*"]"),lwd=2)
text(0,35,"Noon");text(-40,12,"-4h < noon"); text(40,12,"+4h > noon");
title(as.expression(bquote("Lat="*.(lat)*degree*" Day= "*.(nday))),cex.main=0.9)
}
return(list(nday=nday, hr.noon=hr.noon, azi=azi, elev=elev))
}
sun.diagram <- function(lat){
# allocate memory for lists
# later we will put results in lists
azi.path <- list()
elev.path <- list()
# number of days in each month
daymo = c(31,28,31,30,31,30,31,31,30,31,30,31)
# select days to display in diagram
# jan 21, mar 21, june 21
display = c(0,sum(daymo[1:2]),sum(daymo[1:5])) + 21
dayname <- c("Jan 21", "Mar 21", "Jun 21")
# loop for selected days
for(id in 1:length(display)){
# day in the set
nday <- display[id]
# declination
dec <- 23.45*sin((2*pi/365)*(nday-81))
# hr angle and trig
hr.noon <- seq(-12,+12,1)
hr.angle <- 15*hr.noon
sin.H <- sin((pi/180)*hr.angle)
cos.H <- cos((pi/180)*hr.angle)
# elevation and angle
sin.elev <- cos((pi/180)*lat)*cos((pi/180)*dec)*cos.H +
sin((pi/180)*lat)*sin((pi/180)*dec)
elev <- asin(sin.elev)*180/pi
# azimuth
sin.azi <- cos((pi/180)*dec) * sin.H/cos((pi/180)*elev)
# value for testing azimuth
test.tan <- tan(dec*pi/180)/tan(lat*pi/180)
kn <- which(hr.noon==0)
nn <- length(cos.H)
azi <- array()
#loop for hours
for(k in 1:nn){
if(cos.H[k] >= test.tan){
azi[k] <- asin(sin.azi[k])*180/pi
}
else{
if(k <= kn) {corr = -180} else {corr = 180}
azi[k] <- corr - asin(sin.azi[k])*180/pi
}
} # end of hour loop
# store results in list
azi.path[[id]] <- azi
elev.path[[id]] <- elev
} # end of day loop
wd=7; ht=7
panels(wd,ht,1,1,pty="m")
plot(azi.path[[1]],elev.path[[1]],type="b",
xlab=expression("Azimuth ["*degree*"]"),ylab=expression("Elevation ["*degree*"]"),lwd=2,
ylim=c(0,90),xlim=c(-120,120))
for(id in 2:length(display)){
lines(azi.path[[id]],elev.path[[id]],type="b",lwd=2)
}
text(0,30,dayname[1])
text(0,50,dayname[2])
text(0,75,dayname[3])
}
collector <- function(Ibd,sunpath,tilt,azi.c,fr,label=""){
elev <- sunpath$elev
azi <- sunpath$azi
hr.noon <- sunpath$hr.noon
nday <- sunpath$nday
cos.inci <- cos(elev*pi/180)*cos((azi-azi.c)*pi/180)*sin(tilt*pi/180)+
sin(elev*pi/180)*cos(tilt*pi/180)
daymo = c(31,28,31,30,31,30,31,31,30,31,30,31)
sum.mo <-array()
sum.mo[1] <- daymo[1]; if (nday <= sum.mo[1]) n.month <- 1
for(i in 2:12) {
sum.mo[i] <- sum.mo[i-1]+ daymo[i-1]
if (nday <= sum.mo[i] && nday > sum.mo[i-1]) n.month <- i
}
Ib <- Ibd$Ib[n.month]
Id <- Ibd$Id[n.month]
Ibc <- round(Ib*cos.inci,3)
Idc <- round(Id*(1+cos(tilt*pi/180))/2,3)
# reflected
IBH <- Ib*sin(elev*pi/180)
IDH <- Id
Irc <- round(fr*(IBH+IDH)*(1-cos(tilt*pi/180))/2,3)
# total
for(i in 1:length(Ibc)) if(Ibc[i]<=0) Ibc[i] <-0
for(i in 1:length(Irc)) if(Irc[i]<=0) Irc[i] <-0
for(i in 1:length(Idc)) if(Idc[i]<=0) Idc[i] <-0
Ic <- Ibc+Idc+Irc
Ibc.h <- sum(Ibc)
Idc.h <- sum(Idc)
Irc.h <- sum(Irc)
Ic.h <- sum(Ic)
I.h <- c(Ibc.h,Idc.h,Irc.h,Ic.h)
ymax <- 1.2*max(Ic)
par(mar = c(4,4,1,4))
matplot(hr.noon,cbind(Ibc,Ic),type="l",xlab="Hour noon",
ylab="Direct and Total (kW/m2)",ylim=c(0,ymax),lwd=1,col=1,lty=1:2,cex.axis=0.8)
mtext(side=3,line=-1,label)
# plot y2 right hand side axes
par(new=T)
matplot(hr.noon,cbind(Idc,Irc), axes=F, type="l", ylim=c(0,max(Idc+Irc)), xlab=NA,
ylab=NA,lty=3:4,col=1,cex.lab=0.8,cex.axis=0.7,lwd=1.8)
axis(side=4,cex.axis=0.8); mtext(side=4, line=2,"Diffuse and Reflected (kW/m2)")
legend('topright',legend=c("Direct","Total","Diffuse","Reflected"),
lty=1:4,col=1,lwd=1,cex=0.7)
return(list(Ib=Ib,Id=Id,Ibc=Ibc,Idc=Idc,Irc=Irc,Ic=Ic,I.h=I.h))
}
tilt.adj <- function(lat,days,labels){
tilt <- lat - round(declination(plot=F)$dec[days],2)
barplot(tilt,names.arg=labels,ylab=expression("Tilt ["*degree*"]"),ylim=c(0,1.2*max(tilt)),
legend=c(as.expression(bquote("Lat= "*.(lat)*degree))),args.legend=list(x='top'))
return(list(days=days,tilt=tilt))
}
month.prod <- function(dat){
months <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
Ibd <- beam.diffuse(dat,plot=FALSE)
days <- days.mo(21)$day.mo
tilt <- array();Ep <- array(); Ea <- array()
for(i in 1:12){
wd=7; ht=7; panels(wd,ht,2,1,pty="m")
sunpath <- sun.path(lat=Ibd$lat,nday=days[i],plot=FALSE)
label.day <- paste("Day",days[i], ", 21 of",months[i])
Ic <- collector(Ibd, sunpath, tilt=Ibd$lat,
azi.c=0,fr=0.2,label=paste("Polar Mount",label.day))
Ep[i] <- Ic$I.h[4]
tilt[i] <- 90 - max(sunpath$elev)
Ic <- collector(Ibd, sunpath, tilt=tilt[i],
azi.c=0,fr=0.2,label=paste("Tilt Adjusted",label.day))
Ea[i] <- Ic$I.h[4]
}
E <- rbind(Ep,Ea); Emax=1.5*max(E)
wd=7; ht=7; panels(wd,ht,2,1,pty="m")
barplot(tilt,names.arg=months,ylab=expression("Tilt ["*degree*"]"),ylim=c(0,1.2*max(tilt)))
barplot(E,ylim=c(0,Emax),names.arg=months,beside=TRUE,ylab="E (kWh/m2/d)",
legend=c("Polar Mount","Tilt_Adjusted"))
Eadd <- round(100*(Ea-Ep)/Ep,2)
return(list(tilt=tilt, Ea=Ea,Ep=Ep, Eadd=Eadd))
}
two.axis.tracking <- function(dat){
Ibd <- beam.diffuse(dat,plot=FALSE)
days <- days.mo(21)$day.mo
months <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
tilt <- array();Ep <- array(); Ea <- array()
for(i in 1:12){
wd=7; ht=7; panels(wd,ht,2,1,pty="m")
sunpath <- sun.path(lat=Ibd$lat,nday=days[i],plot=FALSE)
label.day <- paste("Day",days[i], ", 21 of",months[i])
Ic <- collector(Ibd, sunpath, tilt=Ibd$lat,
azi.c=0,fr=0.2,label=paste("Polar Mount",label.day))
Ep[i] <- Ic$I.h[4]
Ic <- collector(Ibd, sunpath, tilt=90-sunpath$elev,
azi.c=sunpath$azi,fr=0.2,label=paste("Two-axis Tracking",label.day))
Ea[i] <- Ic$I.h[4]
}
E <- rbind(Ep,Ea); Emax=1.5*max(E)
months <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
wd=7; ht=3.5; panels(wd,ht,1,1,pty="m")
barplot(E,ylim=c(0,Emax),names.arg=months,beside=TRUE,ylab="E (kWh/m2/d)",
legend=c("Polar Mount","Two-axis Tracking"))
Eadd <- round(100*(Ea-Ep)/Ep,2)
mean(Eadd)
return(list(tilt=tilt, Ea=Ea,Ep=Ep, Eadd=Eadd))
}
one.axis.tracking <- function(dat,mode='PNS'){
Ibd <- beam.diffuse(dat,plot=FALSE)
days <- days.mo(21)$day.mo
months <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
tilt <- array();Ep <- array(); Ea <- array()
for(i in 1:12){
wd=7; ht=7; panels(wd,ht,2,1,pty="m")
sunpath <- sun.path(lat=Ibd$lat,nday=days[i],plot=FALSE)
label.day <- paste("Day",days[i], ", 21 of",months[i])
Ic <- collector(Ibd, sunpath, tilt=Ibd$lat,
azi.c=0,fr=0.2,label=paste("Polar Mount",label.day))
Ep[i] <- Ic$I.h[4]
Ic <- collector(Ibd, sunpath, tilt=90-max(sunpath$elev),
azi.c=sunpath$azi,fr=0.2,label=paste("One-axis PNS Tracking",label.day))
Ea[i] <- Ic$I.h[4]
}
E <- rbind(Ep,Ea); Emax=1.5*max(E)
months <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
wd=7; ht=3.5; panels(wd,ht,1,1,pty="m")
barplot(E,ylim=c(0,Emax),names.arg=months,beside=TRUE,ylab="E (kWh/m2/d)",
legend=c("Polar Mount","One-axis PNS Tracking"))
Eadd <- round(100*(Ea-Ep)/Ep,2)
mean(Eadd)
return(list(tilt=tilt, Ea=Ea,Ep=Ep, Eadd=Eadd))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.