• API
• FAQ
• Tools
• Archive
SHARE
TWEET # Untitled a guest Jun 10th, 2017 69 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. ## Libraries ##
2. require(geosphere)  # Generates regular coordinates
3.
4. plotResults <- function(){
5.     ylim  = c(0, 1.05*max(Temps[1:iterB, ]))
6.     lclHr = (lon + t)*12/pi
7.
8.     par(mfrow = c(2, 2))
9.     plot(nSave*(1:iterB)/length(tod),
10.          rowMeans(Temps[1:iterB,, drop = F]),
11.          panel.first = grid(),
12.          xlab        = "Days Elapsed",
13.          ylab        = "Mean Temperature (K)",
14.          type        = "l",
15.          main        = paste("Day: ", d, "  ",
16.                             "Rot: ", 180 + round(t*180/pi), "  ",
17.                             "Tav: ", round(mean(T), 2)))
18.     plot(lat,   T, ylim = ylim, cex = .95, panel.first = grid())
19.     plot(lclHr, T, ylim = ylim, cex = .95, xaxt  =  "n", panel.first = grid())
22.     axis(side = 1, at = axTicks(1), labels = round(adjTicks))
23.     hist(T, col = "Blue", breaks = 20, main = round(range(T), 2))
24. }
25.
26.
27.
28. ## Model for a single point ##
29. # C*dT/dt = I*(1-a) - s*e*T^4
30.
31. ## Simulation Parameters ##
32. T0    = 0    # Initial Temperature (K)
33. dt    = 1    # Timestep (min)
34. nDays = 2e1  # Num days to run simulation
35. nPrts = 5e3  # Num subdivisions of surface
36. nSave = 5e1  # Save every nSave iterations
37.
38.
39. ## Constant Surface Parameters ##
40. S0 = 1370           # Solar Constant (W/m^2)
41. s  = 5.670373e-8    # SB Constant (W*K^-4*m^-2)
42. C  = 1.3e5          # Heat Capacity (W*s*K^-1*m^-2)
43. e  = 1              # Emissivity
44.
45.
46. ## Orbital Parameters ##
47. #LonLat = as.data.frame(regularCoordinates(nPrts))  # Degrees lat/lon
48. LonLat = as.data.frame(randomCoordinates(nPrts))    # Degrees lat/lon
49. We     = 7.2921150e-5                               # Radians/s (Earth)
50. Wm     = We*.0339                                   # Radians/s (Moon)
51.
52.
53. ## Choose rotation rate ##
54. W = c(We, Wm)
55.
56.
57. ## Derived Values ##
58. dt    = dt*60                                                   # Timestep (sec)
59. tod   = seq(-pi, pi, by = dt*W)                                 # Radians/timestep
60. nIter = nDays*length(tod)                                       # Num of total steps
61. lon   = LonLat\$lon*pi/180                                       # Longitude (radians)
62. lat   = LonLat\$lat*pi/180                                       # Latitude (radians)
63. a     = 0.045 + 0.045*(abs(lat)/45)^3 + 0.14*(abs(lat)/90)^8    # Albedo
64.
65.
66.
67. ## Simulation ##
68. T      = rep(0, nrow(LonLat))
69. Temps  = matrix(T0, nrow = ceiling(nIter/nSave), ncol = nrow(LonLat))
70. iterA  = iterB  = 0
71. for(d in 1:nDays){
72.     for(t in tod){
73.         iterA = iterA + 1
74.         I     = pmax(S0*cos(lon + t)*cos(lat), 0)
75.         dT    = (dt/C)*(I*(1-a) - s*e*T^4)
76.         T     = T + dT
77.
78.         if(iterA %% nSave == 0 | iterA == 1){
79.             iterB = iterB + 1
80.             Temps[iterB, ] = T
81.
82.             # png(paste0("moon", iterB,".png"), width = 1920,  height = 1080)
83.             plotResults()
84.             # dev.off()
85.         }
86.     }
87. }