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())
  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 = 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)[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
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top