library("Cairo") # http://www.acrim.com/RESULTS/Composite/composite_acrim_hdr.txt # http://www.acrim.com/RESULTS/Composite/composite_acrim.txt # http://www.pmodwrc.ch/pmod.php?topic=tsi/composite/SolarConstant # ftp://ftp.pmodwrc.ch/pub/data/irradiance/composite/DataPlots/composite_d41_62_1002.dat # http://adsabs.harvard.edu/full/1953AJ.....58..196G # http://www.docstoc.com/docs/9654876/Sun-Insolation # Colums 1: Year # Column 2: TSI (W/m2) # Column 3: TSI Uncertainty (W/m2) # Column 4: Epoch 1980 Day library(lattice) #acrim <- read.table("data/composite_acrim.txt", sep="", col.names=c("year", "tsi", "tsi_err", "epoch_day")) # http://en.wikipedia.org/wiki/File:InsolationTopOfAtmosphere.png # Q the theoretical daily-average insolation at the top of the atmosphere # Solar constant So = 1367 W m−2 # obliquity of ε = 23.4398° # longitude of perihelion of ϖ = 282.895° # eccentricity e = 0.016704. # Contour labels (green) are in units of W m−2. # Declination δ = ε * sin(θ) # (Ro/Re) = RR = 1 + e*cos(θ - ϖ) # θ = poisiton of orbit 0 at the vernal equinox. # φ = latitude of observation # Q = (So/pi)(Ro^2/Re^2)(h*sin(φ)*sin(δ) + cos(φ)*cos(δ)*sin(h)) # Insolation So <- 1367 # obliquity deg ε <- 23.4298*pi/180 # longitude of perihelion ϖ <- 282.895*pi/180 # eccentricity e <- 0.016704 # I like pie # so does R - this is predefined # pi <- 3.14159265 # position in orbit # spring equinox = 0 # θ <- 0 θ <- 270*pi/180 # Ro/Re distance of earth ratio RR <- 1 + e*cos(θ - ϖ) # declination # δ = ε * sin ((360/365)*(284 + day_of_year)) δ = ε * sin(θ) # latitude φ <- 90*pi/180 # hour of sunrise # convert to radians # m = -δ/(pi/2 - θ) if (tan(φ)*tan(δ) > 1) { h <- pi } if (tan(φ)*tan(δ) < -1 ) { h <- 0 } if (tan(φ)*tan(δ) <= 1 & tan(φ)*tan(δ) >= -1 ) { h <- acos(-tan(φ)*tan(δ)) } # calculated insolation daily insolation at latitude φ Q <- (So/pi)*RR*RR*(h*sin(φ)*sin(δ) + cos(φ)*cos(δ)*sin(h)) lat0 <- 121 # 30 N lat1 <- 161 # 70 N day_shift <- -90 # move summer soltice to center of graph north_insol <- array(0, dim=c(366,lat1-lat0)) south_insol <- array(0, dim=c(366,lat1-lat0)) #for (lat in 1:181) { #for (lat in lat0:lat1) { # north latitudes only for (lat in 1:(lat1-lat0)) { # north latitudes only latt <- lat+lat0-91 #for (deg in 1:361) { for (day in 1:366) { #θ <- (360/366)*(284 + day)) #θ <- (360/367)*(285-day)*pi/180 # move from spring equinox to jan 1 θ <- (360/367)*(day-84)*pi/180 # move from spring equinox to jan 1 #θ <- (deg-1)*pi/180 RR <- 1 + e*cos(θ - ϖ) δ = ε * sin(θ) φ <- (latt)*pi/180 if (tan(φ)*tan(δ) > 1) { h <- pi } if (tan(φ)*tan(δ) < -1 ) { h <- 0 } if (tan(φ)*tan(δ) <= 1 & tan(φ)*tan(δ) >= -1 ) { h <- acos(-tan(φ)*tan(δ)) } north_insol[day,lat] <- (So/pi)*RR*RR*(h*sin(φ)*sin(δ) + cos(φ)*cos(δ)*sin(h)) # lets do the S hemisphere too Sφ <- -(latt)*pi/180 if (tan(Sφ)*tan(δ) > 1) { h <- pi } if (tan(Sφ)*tan(δ) < -1 ) { h <- 0 } if (tan(Sφ)*tan(δ) <= 1 & tan(Sφ)*tan(δ) >= -1 ) { h <- acos(-tan(Sφ)*tan(δ)) } south_insol[day,lat] <- (So/pi)*RR*RR*(h*sin(Sφ)*sin(δ) + cos(Sφ)*cos(δ)*sin(h)) } } write.table(t(north_insol),"data/north_insol.txt",row.names=F) # shift insol by 90 days so summer soltice is centered #x <- north_insol[277:366,] #north_insol[91:366,] <- south_insol[1:276,] #north_insol[1:90,] <- x # south: shift insol by 90 days so summer soltice is centered x <- south_insol[1:183,] south_insol[1:183,] <- south_insol[184:366,] south_insol[184:366,] <- x # plot CairoPNG("img/huybers-2006b-fig1-b-rb.png",width=600,height=200) #ytics=c(0,0.167,0.33,0.5,0.667,0.833,1) #ylabels=c(-90,-60,-30,0,30,60,.90) ytics=c(0,1) ylabels=c(lat0,lat1) #xtics=c(0,0.25,0.5,0.75,1) #xlabels=c(0,90,180,270,360) xtics=seq(0,1,by=1/11) xlabels=c("J","F","M","A","M","J","J","A","S","O","N","D") lev=seq(0,600,by=50) colscale=c("darkblue","blue","lightblue","cyan","white","yellow","orange","red","firebrick") #filled.contour(north_insol, asp=1/4, levels=lev, color=colorRampPalette(colscale),ylab="φ", # key.axes=c(1), xlab="θ", plot.axes = { axis(1, xtics, xlabels); axis(2, ytics, ylabels)}) #levelplot(north_insol, cuts=15, ylab="Latitiude", #col.regions=c("black","darkblue","blue","lightblue","white","yellow","gold","orange","salmon","red","darkred")) contourplot(north_insol,cuts=22,region=T,pretty=T,scale=list(x=list(tick.number=0),y=list(labels=c("30","40","50","60","70"))), ylab="Latitude",xlab=xlabels,colorkey=list(space="bottom",tick.number=10),labels=F,main="Northern Hemisphere Daily Insolation", col.regions=colorRampPalette(colscale),aspect=1/5) # SOUTH CairoPNG("img/huybers-2006b-fig1-b-rb-south.png",width=600,height=200) ytics=c(0,1) ylabels=c(-lat0,-lat1) xtics=seq(0,1,by=1/11) xlabels=c("J","F","M","A","M","J","J","A","S","O","N","D") colscale=c("darkblue","blue","lightblue","cyan","white","yellow","orange","red","firebrick") contourplot(south_insol,cuts=22,region=T,pretty=T,scale=list(x=list(tick.number=0),y=list(labels=c("-30","-40","-50","-60","-70"))), ylab="Latitude",xlab=xlabels,colorkey=list(space="bottom",tick.number=10),labels=F,main="Southern Hemisphere Daily Insolation", col.regions=colorRampPalette(colscale),aspect=1/5) CairoPNG("img/huybers-2006b-fig1-b-rb-diff.png",width=600,height=200) contourplot(south_insol-north_insol,cuts=22,region=T,pretty=T,scale=list(x=list(tick.number=0),y=list(labels=c("30","40","50","60","70"))), ylab="Latitude",xlab=xlabels,colorkey=list(space="bottom",tick.number=10),labels=F,main="Difference of Southern and Northern Hemispheres", col.regions=colorRampPalette(colscale),aspect=1/5)