#NASA Diviner Data Source:
#http://pds-geosciences.wustl.edu/lro/lro-l-dlre-4-rdr-v1/lrodlr_1001/data/gcp/
setwd("~/sshen/climmath")
d19 <- read.table("data/tbol_snapshot.pbin4d-19.out-180-0.txt",header = FALSE)
dim(d19)
## [1] 259200 3
#[1] 259200 3 #259200 grid points at 0.5 lat-lon resolution
#259200=720*360, starting from (-179.75, -89.75) going north
#then back to south pole then going north until the end (179.75, 89.75)
m19 <- matrix(d19[,3],nrow = 360)
dim(m19)
## [1] 360 720
library(maps)
## Warning: package 'maps' was built under R version 4.1.3
Lat1 <- seq(-89.75,by = 0.5,len = 360)
Lon1 <- seq(-189.75,by = 0.5, len = 720)
mapmat <- t(m19)
# mapmat <- pmin(mapmat,10)
# mapmat <- mapmat[,seq(length(mapmat[1,]),1)]#, no flipping
#plot.new()
#png(filename = paste("Moon Surface Temperature, Snapshot=", 19,".png"),
# width = 800, height = 400)
int <- seq(0,400,length.out = 40)
rgb.palette <- colorRampPalette(c('skyblue', 'green', 'blue', 'yellow', 'orange',
'pink','red', 'maroon', 'purple', 'black'),interpolate = 'spline')
filled.contour(Lon1, Lat1, mapmat, color.palette = rgb.palette, levels = int,
plot.title = title("Moon Surface Temperature Observed by NASA Diviner,\n Snapshot 19",
xlab = "Longitude", ylab = "Latitude"
,cex.lab = 1.2,cex.axis = 1.2),
plot.axes = {axis(1); axis(2);grid()},
key.title = title(main = "K"))
#dev.off()
#plot.new()
#png(filename = paste("Moon's Equatorial Temperature at Snapshot", 19,".png"),
# width = 600, height = 400)
plot(Lon1,m19[180,],type = "l", col = "red",lwd = 2,
xlab = "Longitude", ylab = "Temperature [K]",
main = "Moon's Equatorial Temperature at Snapshot 19"
,cex.lab = 1.2,cex.axis = 1.2)
text(-100,250,"Nighttime",cex = 2)
text(80,250,"Daytime",cex = 2, col = "orange")
#dev.off()
#plot.new()
par(mar = c(3.5,4,2,0.5))
#png(filename = paste("Moon's Noon Time Meridional Temperature at Snapshot", 19,".png"),
# width = 600, height = 400)
plot(Lat1,m19[,540],type = "l", col = "red",lwd = 2,
xlab = "Latitude", ylab = "Temperature [K]",
main = "Moon's Noon Time Meridional Temperature at Snapshot 19",
cex.lab = 1.2,cex.axis = 1.2)
#dev.off()
bt <- d19[129601:259200,]
aw <- cos(bt[,2]*pi/180)
wbt <- bt[,3]*aw
bta <- sum(wbt)/sum(aw)
bta
## [1] 302.7653
dt <- d19[0:12960,]
aw <- cos(dt[,2]*pi/180)
wdt <- dt[,3]*aw
dta <- sum(wdt)/sum(aw)
dta
## [1] 124.7387
lat <- 0*pi/180
sigma <- 5.670367*10^(-8)
alpha <- 0.12
S <- 1368
ep <- 0.98
k <- 7.4*10^(-4)
h <- 0.4
T0 <- 260
fEBM <- function(T){(1-alpha)*S*cos(lat*pi/180) -(ep*sigma*T^4 + k*(T-T0)/h)}
#Numerically solve this EBM: fEBM = 0
uniroot(fEBM,c(100,420))
## $root
## [1] 383.6297
##
## $f.root
## [1] -0.0002452205
##
## $iter
## [1] 6
##
## $init.it
## [1] NA
##
## $estim.prec
## [1] 6.103516e-05
#Define a piecewise albedo function
a1 <- 0.7
a2 <- 0.3
T1 <- 250
T2 <- 280
ab <- function(T) {ifelse(T < T1, a1, ifelse(T < T2,((a1-a2)/(T1-T2))*(T-T2)+ a2, a2))}
#Define the range of temperature
T <- seq(200,350,len = 200)
#Plot the albedo as a nonlinear function of T
#setwd("/Users/sshen/climmath")
#png(file = "fig05-05.png", width = 400, height = 300)
par(mar = c(4,4,4,3))
plot(T, ab(T), type = "l", lwd = 2.0,
ylim = c(0,1), xlab = "Surface Temperature [K]",
ylab = "Albedo", main = "Nonlinear Albedo Feedback"
,cex.lab = 1.2,cex.axis = 1.2)