Guest User

Untitled

a guest
Jun 16th, 2017
105
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())
  20. adjTicks = axTicks(1)
  21. adjTicks = ifelse(adjTicks > 0, adjTicks, 24 + adjTicks)
  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 = 2e2 # Num days to run simulation
  35. nPrts = 5e3 # Num subdivisions of surface
  36. nSave = 5e3 # 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.3e8 # 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)[2]
  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. }
  88. # head(Temps)
RAW Paste Data