So now I want to the turn the Daisyworld model into a Monte Carlo Metropolis style calculation. This means a model that you use random events to model reality instead of solving differential equations. This is some exciting because we can get a nice visual graphic for how it works.
Time to have some fun.
The first thing I want to do is to look how the probability of daisies growing depends on temperature. As we discussed in the original Daisyworld model. There is this probability function beta that depends on the local temperature.
In a new python window, I created a function that determines the probability beta as a function of temperature.
import numpy as np
import matplotlib.pyplot as plt
p_k = 1/(17.5)**2 # K - parameter
p_Topt = 295.5 # optimal daisy growth
p_q = 2.1e9 # q factor in model
p_aw = 0.75 #albeto for white daisy
p_ab = 0.25 #albeto for black daisy
p_ag = 0.50 #albeto for ground
# Get the probability of growing a daisy of a certain type.
def Get_Daisy_Prob (aA, aTemperature):
#White Daisy Factor
tt = np.arange(220, 370, 0.1)
TLocW = (p_q * (aA-p_aw) + tt**4)**0.25
WhiteProb = 1 - p_k * (TLocW-p_Topt)**2
mm = WhiteProb < 0
WhiteProb[mm] = 0
NoWhitePRob = 1 - WhiteProb
#Black Daisy Factor
TLocB = (p_q * (aA-p_ab) + tt**4)**0.25
BlackProb = 1 - p_k * (TLocB-p_Topt)**2
mm = BlackProb < 0
BlackProb[mm] = 0
NoBlackPRob = 1 - BlackProb
GroundProb = NoWhitePRob * NoBlackPRob
#Normalize Daisy Growth prob
mag = np.sqrt(WhiteProb**2 + BlackProb**2) + .0001
WP = WhiteProb/mag / np.sqrt(2)
BP = BlackProb/mag / np.sqrt(2)
plt.plot (tt,(1-GroundProb)*WP)
plt.plot (tt,(1-GroundProb)*BP)
plt.xlabel ('Temperature (K)')
plt.ylabel ('Probability of growth')
Prob = np.array([(1-GroundProb)*WP, (1-GroundProb)*BP])
return Prob
A = 0.5
print (Get_Daisy_Prob (A, 300))
Here the orange line is the probability for the black daisies and blue is the white daisies.
[discussion question 7] Compare and contrast this calculation to what we did in the original Daisyworld model. Is there agreement, are there questions? In addition to your homework, please add the questions to the comments of this page.
There are a few things to unpack here. The line of code that says WhiteProb = 1 – p_k * (TLocW-p_Topt)**2 determines the probability that just white daisies would grow. But what about the black daisies? To solve this I did some probability math. First I calculated what the probability that there would be no white daisies (NoWhitePRob) and no black daisies (BlackProb). So to determine the probability that there were any daisies at all, I multiplied them together to get GroundProb.
So the probability of getting a daisy is just (1-GroundProb).
Next I had to figure out if it was a white daisy or a black daisy. I toyed with a few different ideas but in the end I stuck with trig. Determine the magnitude and the factor of each part that for each type of daisy.
mag = np.sqrt(WhiteProb**2 + BlackProb**2) + .0001
WP = WhiteProb/mag / np.sqrt(2)
BP = BlackProb/mag / np.sqrt(2)Next I had to convert this function so that it only took one temperature. That was fairly straightforward. See the work below:
# Get the probability of growing a daisy of a certain type.
def Get_Daisy_Prob (aA, aTemperature):
#White Daisy Factor
TLocW = (p_q * (aA-p_aw) + aTemperature**4)**0.25
WhiteProb = 1 - p_k * (TLocW-p_Topt)**2
if WhiteProb < 0:
WhiteProb = 0
NoWhitePRob = 1 - WhiteProb
#Black Daisy Factor
TLocB = (p_q * (aA-p_ab) + aTemperature**4)**0.25
BlackProb = 1 - p_k * (TLocB-p_Topt)**2
if BlackProb < 0:
BlackProb = 0
NoBlackPRob = 1 - BlackProb
GroundProb = NoWhitePRob * NoBlackPRob
#Normalize Daisy Growth prob
mag = np.sqrt(WhiteProb**2 + BlackProb**2) + .0001
WP = WhiteProb/mag / np.sqrt(2)
BP = BlackProb/mag / np.sqrt(2)
ratio = WP / (WP + BP)
Prob = np.array([ratio, GroundProb])
return Prob[discussion question 8] Compare and contrast this function with the one above. What are the major differences?
So now we had to choose whether a daisy lives or dies. To do this we wrote the following block of code:
#flag for seeing if daisy should be recalculated
GrowNewDaisy = 0
#Is there a daisy alive?
if daisy == 0:
GrowNewDaisy = 1
#Now there is a daisy, but does it die?
else:
#Does Daisy die?
pickNum = np.random.rand()
if pickNum < p_DeathRate:
GrowNewDaisy = 1 There are three states, daisy = -1 (white), daisy = +1 (black), and daisy = 0 (no daisy). If there is no daisy on the ground then the ground is available for daisy growth. Else if there is a daisy there we have to determine if it dies so a new daisy can grow. To do this, we pick a random number and compare it to the death rate. If the daisy dies, this it is available for growth.
[discussion question 9] Walk through the logic for a couple of scenerios and show that it gives the correct answer.
So now if there is available ground for the daisy to grow, then it should grow a new daisy. Right?
#Grow New Daisy
if GrowNewDaisy:
daisy = 0
BWGrowRates = Get_Daisy_Prob (aA, Temp)
pickNum = np.random.rand()
if pickNum > BWGrowRates[1]:
pickNum = np.random.rand()
if pickNum < BWGrowRates[0]:
daisy = -1
else:
daisy = +1BWGrowRates come from the Get_Daisy_Prob function discussed above. The first term BWGrowRates[1] is the probability anything will grow. Whereas the BWGrowRates[0] is the probability that white daisies will grow vs. black daisies.
[discussion question 10]Translate the code block above to your own words. What is going on there?
Putting it all together!
In the first attempt at the Metropolis Daisyworld Model, I created a grid of points that were different kinds of daisies at different places.
[discussion question 11] Compare this final code (below) to the code we did last time in the first attempt at the Metropolis Daisyworld Model. What is happening in line 160?
#implement daisy model for one point on earth
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import gaussian_filter1d
p_k = 1/(17.5)**2 # K - parameter
p_Topt = 295.5 # optimal daisy growth
p_q = 2.1e9 # q factor in model
p_SL0 = 1100 #solar output factor
p_sig = 5.67e-8 #Stefan-Boltzmann constant
p_q = 2.1e9 # q factor in model
p_aw = 0.75 #albeto for white daisy
p_ab = 0.25 #albeto for black daisy
p_ag = 0.50 #albeto for ground
p_GrowRate = 1
p_DeathRate = 0.1
xLength = 100
yLength = 100
numTimes = 100
time = np.arange(0,numTimes)
SL = p_SL0 * (1 + 0.2 * np.cos(time/20))
Daisy = np.zeros((yLength,xLength))
Temp = np.zeros((yLength,xLength))
EarthTemp = np.zeros(numTimes)
# Get the probability of growing a daisy of a certain type.
def Get_Daisy_Prob (aA, aTemperature):
#White Daisy Factor
TLocW = (p_q * (aA-p_aw) + aTemperature**4)**0.25
WhiteProb = 1 - p_k * (TLocW-p_Topt)**2
if WhiteProb < 0:
WhiteProb = 0
NoWhitePRob = 1 - WhiteProb
#Black Daisy Factor
TLocB = (p_q * (aA-p_ab) + aTemperature**4)**0.25
BlackProb = 1 - p_k * (TLocB-p_Topt)**2
if BlackProb < 0:
BlackProb = 0
NoBlackPRob = 1 - BlackProb
GroundProb = NoWhitePRob * NoBlackPRob
#Normalize Daisy Growth prob
mag = np.sqrt(WhiteProb**2 + BlackProb**2) + .0001
WP = WhiteProb/mag / np.sqrt(2)
BP = BlackProb/mag / np.sqrt(2)
ratio = WP / (WP + BP)
Prob = np.array([ratio, GroundProb])
return Prob
def getDaisySitutation (aA, Temp, daisy):
#flag for seeing if daisy should be recalculated
GrowNewDaisy = 0
#Is there a daisy alive?
if daisy == 0:
GrowNewDaisy = 1
#Now there is a daisy, but does it die?
else:
#Does Daisy die?
pickNum = np.random.rand()
if pickNum < p_DeathRate:
GrowNewDaisy = 1
#Grow New Daisy
if GrowNewDaisy:
daisy = 0
BWGrowRates = Get_Daisy_Prob (aA, Temp)
pickNum = np.random.rand()
if pickNum > BWGrowRates[1]:
pickNum = np.random.rand()
if pickNum < BWGrowRates[0]:
daisy = -1
else:
daisy = +1
return daisy
def getSmoothedArray (Y):
sY = gaussian_filter1d(Y,3)
return sY
for tt in range(numTimes):
TLat = np.zeros(yLength)
Albeto = np.zeros(yLength)
for yy in range(yLength):
#Dtermine the Albeto for each lat.
white = 0
black = 0
ground = 0
SolarRad = (0.5*(1 - (2*np.abs(yy - yLength/2))**2 / yLength**2) + 0.5) * SL[tt]
for xx in range(xLength):
if Daisy[yy,xx] < 0:
white = white + 1
elif Daisy[yy,xx] > 0:
black = black + 1
else:
ground = ground + 1
Albeto[yy] = white/xLength * p_aw + black/xLength * p_ab + ground/xLength * p_ag
#EarthTemperature note power to the fourth
TLat[yy] = (SolarRad * (1-Albeto[yy]) / p_sig)**0.25
#Albeto = getSmoothedArray(Albeto)
TLat = getSmoothedArray(TLat)
for yy in range(yLength):
for xx in range(xLength):
Temp[yy,xx] = TLat[yy]
Daisy[yy,xx] = getDaisySitutation (Albeto[yy], TLat[yy], Daisy[yy,xx])
plt.imshow(Daisy)
plt.show()
EarthTemp[tt] = np.mean(Temp)
#plt.plot(TLat)
plt.plot(EarthTemp)
Check out, some of the pictures at different times? The sun’s energy is changing by 20% over this cycle.
You can defintely see changes. Note that the green is ground, the yellow spots are black daisies and blue spots are white daisies. And yes, I should have spent more time figuring out the other scheme.
So even though the sun’s energy density was changing by 20%, the temperature only changes by a little.
[discussion question 12] Did I include the latitude effects? If so how so? Do you agree?
[discussion question 13] Which model is better original Daisyworld model or this one? Discuss. There is no right or wrong answer.