#to adapt change these eight first lines only
locationname="Burgauberg"
phi=47.162
lambda=16.130
out="" #""/png/jpg/svg/pdf
type="cart" #cart/polar
lang="de" #de/en/it/nl
locale="German_Austria" #month names German_Austria/English/Italian/Dutch
ieqot=F #F/T (do not take equation of time into account)
summert=0 #1/0 (summertime normal time)
if (lang=="de") {
ltext=c("N","O","S","W","Nord","Ost","Süd","West","Azimuth","Höhenwinkel","MEZ","MESZ","Uhr","Sonnenstandsdiagramm","Zenit","Horizont")
} else if (lang=="en") {
ltext=c("N","E","S","W","north","east","south","west","azimuth","elevation angle","CET","CEST","","sun chart","Zenith","Horizon")
} else if (lang=="it") {
ltext=c("N","E","S","O","Nord","Est","Sud","Ovest","azimut","altezza","CET","CEST","","Diagramma solare","Zenit","Orizzonte")
} else if (lang=="nl") {
ltext=c("N","O","Z","W","NB","OL","ZB","WL","Azimut","Hoogte","CET","CEST","","Zonnestanddiagram","Zenit","Horizon")
}
phi=phi*pi/180
lambda=lambda*pi/180
KAPPA=0.016700
EPSILON=23.4385*pi/180
ALPHA=(360-283.067)*pi/180 #77.6°
Sys.setlocale("LC_TIME", locale)
library("Cairo") #for nice plot into files
require("plotrix") #for polar plots (and other stuff?)
require("geonames") #to find the perfect timezone
TZ=as.numeric(GNtimezone(phi*180/pi,lambda*180/pi)[3]) #manual: round((lambda*180/pi)/15)
if (!exists(x="TZ",mode="numeric")) {
TZ=round((lambda*180/pi)/15)
}
h=function(t,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA) {
return( sin(t+delta(t,KAPPA)-ALPHA)*sin(EPSILON) )
}
r=function(t,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA) {
return( sqrt(1-h(t,KAPPA,EPSILON,ALPHA)^2) )
}
mu=function(t,ignore=F,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA) {
zaehl=-sin(delta(t,KAPPA))+tan(EPSILON/2)^2*sin(2*(t-ALPHA)+delta(t,KAPPA))
nenn = cos(delta(t,KAPPA))+tan(EPSILON/2)^2*cos(2*(t-ALPHA)+delta(t,KAPPA))
if (ignore) { res=0 } else { res=atan(zaehl/nenn) } #ignore the eq of time
return( res )
}
delta=function(t,KAPPA=KAPPA) {
return( 2*KAPPA*sin(t)+5/4*KAPPA^2*sin(2*t) )
}
coords=function(t,S,phi,lambda,ignore=ieqot,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA,tz=TZ) {
s=pi/180*15*(S+lambda*180/pi/15-tz)
x1= r(t,KAPPA,EPSILON,ALPHA)*sin(s+mu(t,ignore,KAPPA,EPSILON,ALPHA))
x2= r(t,KAPPA,EPSILON,ALPHA)*cos(s+mu(t,ignore,KAPPA,EPSILON,ALPHA))*sin(phi)+h(t,KAPPA,EPSILON,ALPHA)*cos(phi)
x3=-r(t,KAPPA,EPSILON,ALPHA)*cos(s+mu(t,ignore,KAPPA,EPSILON,ALPHA))*cos(phi)+h(t,KAPPA,EPSILON,ALPHA)*sin(phi)
return(cbind(x1,x2,x3))
}
#returns the angle for sperical coordinates
coord2angle=function(coord,filterout=TRUE) {
diag=sqrt(coord[,1]^2+coord[,2]^2)
psi=atan2(coord[,2],coord[,1])
phi=pi/2-atan(coord[,3]/diag)
azi=pi/2-psi
hoeh=pi/2-phi
azi[azi<0]=azi[azi<0]+2*pi
if (filterout) {
hoeh[hoeh<0]=NA
azi[hoeh<0]=NA
}
return( 180/pi*cbind(azi, hoeh) )
##hoeh=wasserglas(alpha=hoeh,D=wgD,d=wgd,n1=wgn1,n2=wgn2,n3=wgn3)
##return( cbind(azi*180/pi, hoeh) )
}
#transforms an angle to a formated date
angle2date=function(w){
day=(w/(2*pi))*365
dates=as.Date("2009-01-01")+day
if (lang=="de"){
fdate=format(dates, format="%d. %b")
} else {
fdate=format(dates, format="%d %b")
}
return( fdate )
}
#gets the coordinates and returns a nice desciptive string
location=function(long,lat,name) {
long=long*180/pi
if (long<0) { longStr=ltext[8] } else { longStr=ltext[6] }
lon1=trunc(long+0.00001) #to avoid output 14°60' for 15°0'
lon2=60*(long-lon1)
lat=lat*180/pi
if (lat <0) { latStr =ltext[7] } else { latStr=ltext[5] }
lat1=trunc(lat+0.00001) #to avoid output 14°60' for 15°0'
lat2=60*(lat-lat1)
latlon=paste(abs(lat1),"°",abs(round(lat2,digits=1)),as.character("' "),latStr,", ",abs(lon1),"°",abs(round(lon2,1)),as.character("' "),longStr,sep="")
string=paste(name,latlon,sep=", ")
return(string)
}
#Analemma for all hours
analemmabogen=function(suntype="cart",phi=phi,lambda=lambda,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA) {
jz=seq(0,2*pi,length=365*2)-9/365*2*pi
sz=seq(0,23,by=1)
for (si in 1:length(sz)) {
ko=coords(t=jz,S=sz[si],phi=phi,lambda=lambda,ignore=ieqot,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA,tz=TZ)
ang=coord2angle(ko,filterout=TRUE)
mitte=ceiling(length(jz)/2)
sunlines(suntype=type,x=ang[1:mitte,1],y=ang[1:mitte,2],type="l",lwd=3,col="darkgreen")
if (!ieqot) { sunlines(suntype=type,x=ang[mitte:length(jz),1],y=ang[mitte:length(jz),2],type="l",col="olivedrab",lwd=2,lty="dotdash") }
#labels of the full hour markers
stui=substitute(szsi^{h},list(szsi=sz[si]+summert))
yl=ang[mitte,2]+2.5
if(sz[si]==12) {
if(summert==1) { sozeText=ltext[12] } else { sozeText=ltext[11] }
stui=paste(12+summert,"h ",ltext[13]," ",sozeText,sep="")
if (suntype=="cart") {
yl=ang[mitte,2]+4
} else if (suntype=="polar") {
yl=ang[mitte,2]+7
}
}
if (suntype=="cart") {
text(x=1.025*ang[mitte,1]-4,y=yl,labels=stui,cex=1.7,font=2)
} else if (suntype=="polar") {
text(x=(90-yl)*sin((1.025*ang[mitte,1]-4)*pi/180),y=(90-yl)*cos((1.025*ang[mitte,1]-4)*pi/180),labels=stui,cex=1.3,font=2)
}
}
}
#arc of the sun at one date for different days within a year
tagesbogen=function(suntype=type,jahr,lty="solid",lwd=2,col="steelblue",phi=phi,lambda=lambda,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA) {
tag=seq(0,24,length=24*60+1)
for (jt in 1:length(jahr)) {
tagbogen=coords(t=jahr[jt],S=tag,phi=phi,lambda=lambda,ignore=ieqot,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA,tz=TZ)
tagb=coord2angle(tagbogen,filterout=TRUE)
sunlines(suntype=suntype,x=tagb[,1],y=tagb[,2],lty=lty,lwd=lwd,col=col)
}
}
#highlight the position of the sun on a specific date and time as a red point
highlightpoint=function(suntype="cart",timeofday,date="2009-05-01",phi=phi,lambda=lambda,shift=c(0,-3),col="red",KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA,summert=0) {
reddate=as.Date(date)
dayofyear=as.numeric(reddate-as.Date("2009-01-01"))
co=coords(t=dayofyear/365*2*pi,S=timeofday,phi=phi,lambda=lambda,ignore=ieqot,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA,tz=TZ)
winkel=coord2angle(co,filterout=TRUE)
sunlines(suntype=type,x=winkel[,1],y=winkel[,2],type="p",lwd=2,col=col)
if(summert==1) { timeofday=timeofday+1 }
minofday=formatC(round((timeofday-floor(timeofday))*60), digits=0, width=2, format="f", flag="0")
timestring=paste(floor(timeofday),":",minofday,sep="")
if (suntype=="cart") {
text(winkel+shift,paste(format(reddate,"%d. %b"),"\n",timestring),col=col,cex=1.4,adj=c(0.5,0.5))
} else if (suntype=="polar") {
pwinkel=winkel+shift
text(x=(90-pwinkel[2])*sin(pwinkel[1]*pi/180),y=(90-pwinkel[2])*cos(pwinkel[1]*pi/180),paste(format(reddate,"%d. %b"),"\n",timestring),col=col,cex=1,adj=c(0.5,0.5))
}
print(winkel)
return(winkel)
}
#function that handles both polar and cartesian coordinates
sunlines=function(suntype,x,y=NULL,...) {
if (suntype=="cart") {
lines(x,y, ...)
} else if (suntype=="polar") {
lines( x=(90-y)*sin(x*pi/180),y=(90-y)*cos(x*pi/180), ...)
}
}
filename=paste("Sunchart-",type,"-",lang,"-",locationname,".",out,sep="")
if (out=="png") {
CairoPNG(file=filename,width=1850,height=690)
#png(file=filename,width=1850,height=690)
#bitmap(file=filename,type="png16m",width=1850,height=690)
} else if (out=="pdf") {
pdf(file=filename,paper="a4r",width=0,height=0)
##pdf(file=filename,width=(wgD+2*wgd)*pi/2.54,height=5/2.54)
} else if (out=="jpg") {
#CairoJPEG(file=filename,width=1850,height=690)
jpeg(file=filename,width=1850,height=690,quality=98)
} else if (out=="svg") {
CairoSVG(file=filename,width=11,height=11)
}
xlab=c(paste(ltext[1],sep=""),paste(ltext[1],ltext[1],ltext[2],sep=""),paste(ltext[1],ltext[2],sep=""),paste(ltext[2],ltext[1],ltext[2],sep=""),paste(ltext[2],sep=""),paste(ltext[2],ltext[3],ltext[2],sep=""),paste(ltext[3],ltext[2],sep=""),paste(ltext[3],ltext[3],ltext[2],sep=""),paste(ltext[3],sep=""),paste(ltext[3],ltext[3],ltext[4],sep=""),paste(ltext[3],ltext[4],sep=""),paste(ltext[4],ltext[3],ltext[4],sep=""),paste(ltext[4],sep=""),paste(ltext[4],ltext[1],ltext[4],sep=""),paste(ltext[1],ltext[4],sep=""),paste(ltext[1],ltext[1],ltext[4],sep=""),paste(ltext[1],sep="") )
if (type=="cart") {
plot(c(0,0),c(0,0),type="n",xlim=c(45,315),ylim=c(0,min(90-phi*180/pi+31,90)),xaxs="i",yaxs="i", axes=FALSE, main=ltext[14], cex.main=2.3, font.main=2, xlab=ltext[9], ylab=ltext[10], cex.lab=1.33)
##par(mar=c(0,0,0,0),cex=0.33)
##plot(c(0,0),c(0,0),type="n",xlim=c(0,360),ylim=c(-5,0),xaxs="i",yaxs="i", axes=FALSE)
##box("figure", col="blue",lwd=2)
#par(mai=c(0.95625,0.76875,0.76875,0.39375))
#box("figure", col="blue")
legend(x=45,y=min(90-phi*180/pi+31,90),legend=location(lat=phi,long=lambda,name=locationname),box.col="white",bg="white",cex=1.7)
xticks=seq(par("usr")[1],par("usr")[2],by=22.5)
yticks=seq(par("usr")[3],par("usr")[4],by=10)
axis(1, at=xticks[1:length(xticks)], tick=F, cex.axis=1.33, labels=xlab[3:15],font=2)
##axis(1, at=xticks[1:length(xticks)], tck=0.01, mgp=c(0,-1.5,0), cex.axis=1.75, labels=xlab[1:length(xticks)],font=2)
axis(2, at=yticks, tick=F, cex.axis=1.5, labels=paste(yticks,"°",sep=""))
grid( nx=length(xticks)-1,ny=NA,col="gray75")
abline( h=yticks, col="gray75", lty="dotted")
abline( v=c(90,180,270), col="gray20",lty="dotted")
box(which="plot",col="grey50")
} else if (type=="polar") {
polar.plot(length=NA,polar.pos=NA,rp.type="p",clockwise=TRUE,label.pos=seq(45,360+22.5,by=22.5),start=135,labels=xlab,radlab=F,radial.labels="",radial.lim=c(0,30,60,90),show.centroid=F, main=ltext[14], cex.main=2.3, font.main=2 )
text(x=c(0,0,0,0),y=c(0,30,60,90),labels=c(ltext[15],"60°","30°",ltext[16]),col="darkgrey")
legend("topleft",legend=location(lat=phi,long=lambda,name=locationname),box.col="white",bg="white",cex=1)
}
#arc on a day: birthdays
jahr=c(113)/365*2*pi
tagesbogen(suntype=type,jahr,lty=4,lwd=2,col="orange",phi=phi,lambda=lambda,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA)
analemmabogen(suntype=type,phi=phi,lambda=lambda,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA)
#daily arc: first half of the year
jahr=cumsum(c(0,31,28,31,30,31))/365*2*pi #one line per month
tagesbogen(suntype=type,jahr,lty="solid",lwd=2,col="steelblue",phi=phi,lambda=lambda,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA)
#daily arc: second half of the year
jahr=cumsum(c(181,31,31,30,31,30))/365*2*pi #one line per month
tagesbogen(suntype=type,jahr,lty="dotted",lwd=2,col="skyblue3",phi=phi,lambda=lambda,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA)
#one point per month, each at 12 o'clock
jp=cumsum(c(0,31,28,31,30,31,30,31,31,30,31,30))/365*2*pi #one line per month
jahresbogen=coords(t=jp,S=12,phi=phi,lambda=lambda,ignore=ieqot,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA,tz=TZ)
angl=coord2angle(jahresbogen,filterout=TRUE)
sunlines(suntype=type,x=angl[,1],y=angl[,2],type="p",lwd=2)
adj=cbind(c(0.9,1,1.1,1.1,0,-0.2,1.2,1.1,-0.1,0,0,0.1),c(1.5,1.7,1.6,-0.2,1.3,1,-0.4,0.7,-0.3,-0.7,-0.7,-0.5))
for (i in 1:length(jp)) {
if (type=="cart") {
text(x=angl[i,1],y=angl[i,2],angle2date(jp[i]),adj=adj[i,],cex=1.33)
} else if (type=="polar") {
text(x=(90-angl[i,2])*sin(angl[i,1]*pi/180),y=(90-angl[i,2])*cos(angl[i,1]*pi/180),angle2date(jp[i]),adj=adj[i,],cex=0.7)
}
}
#sun location on a certain day and time
hipo=highlightpoint(suntype=type,timeofday=13+25/60,date="2009-04-23",phi=phi,lambda=lambda,shift=c(0,-3),col="brown4",KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA,summert=0)
if (out=="png" | out=="jpg" | out=="svg" | out=="pdf") {
dev.off()
print(paste("output written to ",getwd(),"/",filename,sep=""))
}
###
###
## ev paar Zeilen rauf schieben, dient für sonneneinstrahlung abzuschätzen
sunlines(suntype=type,x=c(113,113+180),y=c(0,0),type="l",lwd=3,col="red")
# utitlity function to integrate (is sum up) the sun intensity over one day (and then to sum it up over a year)
besonnung_t=function(jahr) {
tag=seq(0,24,length=60*24+1)
tagbogen=coords(t=jahr,S=tag,phi=phi,lambda=lambda,ignore=ieqot,KAPPA=KAPPA,EPSILON=EPSILON,ALPHA=ALPHA)
tagb=coord2angle(tagbogen,filterout=TRUE)
resu=tagb[,2]
resu[is.na(resu)]=0
resu[tagb[,1] < (113.5) ]=0 #SüdWand noch nicht besonnt
resu[tagb[,1] > (113.5+180) ]=0 #SüdWand nicht mehr besonnt
return( mean(resu) )
}
# utitlity function to integrate (is sum up) the sun intensity over the whole year)
besonnung_j=function() {
sonnensumme=0
jahr=seq(0,2*pi/365,length=365)
for (jt in 1:length(jahr)) {
sonnensumme=sonnensumme+besonnung_t(jahr[jt])
}
return(sonnensumme/length(jahr))
}
besonnung_t(jahr=c(94)/365*2*pi)
besonnung_j()
###
###