## Libraries ##
require(geosphere)  # Generates regular coordinates

plotResults <- function(){
    ylim  = c(0, 1.05*max(Temps[1:iterB, ]))
    lclHr = (lon + t)*12/pi

    par(mfrow = c(2, 2))
    plot(nSave*(1:iterB)/length(tod),
         rowMeans(Temps[1:iterB,, drop = F]),
         panel.first = grid(),
         xlab        = "Days Elapsed",
         ylab        = "Mean Temperature (K)",
         type        = "l",
         main        = paste("Day: ", d, "  ",
                            "Rot: ", 180 + round(t*180/pi), "  ",
                            "Tav: ", round(mean(T), 2)))
    plot(lat,   T, ylim = ylim, cex = .95, panel.first = grid())
    plot(lclHr, T, ylim = ylim, cex = .95, xaxt  =  "n", panel.first = grid())
    adjTicks = axTicks(1)
    adjTicks = ifelse(adjTicks > 0, adjTicks, 24 + adjTicks)
    axis(side = 1, at = axTicks(1), labels = round(adjTicks))
    hist(T, col = "Blue", breaks = 20, main = round(range(T), 2))
}



## Model for a single point ##
# C*dT/dt = I*(1-a) - s*e*T^4

## Simulation Parameters ##
T0    = 0    # Initial Temperature (K)
dt    = 1    # Timestep (min)
nDays = 2e2  # Num days to run simulation
nPrts = 5e3  # Num subdivisions of surface
nSave = 5e3  # Save every nSave iterations


## Constant Surface Parameters ##
S0 = 1370           # Solar Constant (W/m^2)
s  = 5.670373e-8    # SB Constant (W*K^-4*m^-2)
C  = 1.3e8          # Heat Capacity (W*s*K^-1*m^-2)
e  = 1              # Emissivity


## Orbital Parameters ##
#LonLat = as.data.frame(regularCoordinates(nPrts))  # Degrees lat/lon
LonLat = as.data.frame(randomCoordinates(nPrts))    # Degrees lat/lon
We     = 7.2921150e-5                               # Radians/s (Earth)
Wm     = We*.0339                                   # Radians/s (Moon)


## Choose rotation rate ##
W = c(We, Wm)[2]


## Derived Values ##
dt    = dt*60                                                   # Timestep (sec)
tod   = seq(-pi, pi, by = dt*W)                                 # Radians/timestep
nIter = nDays*length(tod)                                       # Num of total steps
lon   = LonLat$lon*pi/180                                       # Longitude (radians)
lat   = LonLat$lat*pi/180                                       # Latitude (radians)
a     = 0.045 + 0.045*(abs(lat)/45)^3 + 0.14*(abs(lat)/90)^8    # Albedo



## Simulation ##
T      = rep(0, nrow(LonLat))
Temps  = matrix(T0, nrow = ceiling(nIter/nSave), ncol = nrow(LonLat))
iterA  = iterB  = 0
for(d in 1:nDays){
    for(t in tod){
        iterA = iterA + 1
        I     = pmax(S0*cos(lon + t)*cos(lat), 0)
        dT    = (dt/C)*(I*(1-a) - s*e*T^4)
        T     = T + dT

        if(iterA %% nSave == 0 | iterA == 1){
            iterB = iterB + 1
            Temps[iterB, ] = T

            # png(paste0("moon", iterB,".png"), width = 1920,  height = 1080)
            plotResults()
            # dev.off()
        }
    }
}
# head(Temps)