On Fri, 2008-05-02 at 08:33 -0400, Michael Droettboom wrote:
> I converted your example to a standalone one (attached), and it does
> work for me. It's hard to say what may be going wrong for you without
> knowing what a_steady and L_range are and what the limits of the axes
> are getting set to etc. Can you provide a fully self-contained script
> that exhibits this behavior? Also, what version of matplotlib and
> platform are you using?
>
> Cheers,
> Mike
Hello:
I am running Debian sid on Powerpc and I have the Debian sid
python-matplotlib version 0.91.2-2. Your stand alone script also works
for me. Strange.
I have included the script I am using set to plot 5 lines. The script is
designed to solve the Daisyworld model proposed by James Lovelock in the
early 80's. Thank you for your help!
Zoho
#!/usr/bin/env python
from __future__ import division
import pdb
import numpy as np
from numpy import sum, zeros, less, where
import scipy as sp
#import matplotlib
#matplotlib.use('PS')
import matplotlib.pyplot as plt
from matplotlib import collections
from matplotlib.colors import colorConverter
def main():
"""plot the results of daisyworld equations"""
S = 952.56 # solar constant
q = 2.06425e9 # relates daisy temp to its albedo
gamma = 0.3 # Death rate for daisy
sigma = 5.67e-8 # Stefan-Boltzman constant
rho = 3.265e-3 # Modifier for the growth curve
kappa = 295.5 # Temperature for optimal growth in Kelvin
Ax = 0.50 # Albedo for bare ground
types = 5 # number of daisy types
A = np.linspace(0.25, 0.75, types) # equally spaced albedo's for daisies
N = 300 # Limit for the L range
M = 10 # Limit for powerseries loop
t = 1 # step size
L_initial = 0.5 # starting level for L
L_final = 1.7 # finish level for L
L_range = np.linspace(L_initial, L_final, N+1).tolist() # the range for L as a list
# cyclic L to illustrate hysterisis
# L0 = linspace(L_initial, L_final, N+1)
# L1 = arange(L_final, L_initial, N+1)
# L = np.hstack((L0,L1)).tolist()
# create the arrays for the powerseries
a = zeros((4,types))
x = zeros((3))
beta = zeros((3,types))
T = zeros((3,types))
Ap = zeros((3))
Tp = zeros((3))
# lists to hold steady state
a_steady = []
Tp_steady = []
# The main loop that iterates and stores pertinent values
for L in L_range:
# reset the intial conditions for phenotypes which have died.
indices = less(a[0], 0.001)
a[0] = where(indices, 0.001, a[0])
for iteration in range(1, M):
x[0] = 1 - sum(a[0,:])
Ap[0] = x[0]*Ax + sum(a[0,:]*A)
Tp[0] = (S*L/sigma*(1 - Ap[0]))**(1./4)
T[0,0:types] = (q*(Ap[0] - A[0:types]) + Tp[0]**4)**(1./4)
beta[0,0:types] = np.e**(-rho*(kappa - T[0,0:types])**2)
a[1,0:types] = a[0,0:types]*(x[0]*beta[0,0:types] - gamma)
x[1] = -sum(a[1,:])
Ap[1] = x[1]*Ax + sum(a[1,:]*A)
Tp[1] = -1./(4*Tp[0]**3)*(S*L/sigma)*Ap[1]
T[1,0:types] = 1./(4*T[0,0:types]**3)*(q*Ap[1] + 4*Tp[0]**3*Tp[1])
beta[1,0:types] = beta[0,0:types]*2*rho*(kappa - T[0,0:types])*T[1,0:types]
a[2,0:types] = 1./2*(a[0,0:types]*(x[0]*beta[1,0:types] + x[1]*beta[0,0:types]) + a[1,0:types]*(x[0]*beta[0,0:types] - gamma))
x[2] = -sum(a[2,:])
Ap[2] = x[2]*Ax + sum(a[2,:]*A)
Tp[2] = -1./(4*Tp[0]**3)*((S*L/sigma)*Ap[2] + 6*(Tp[0]*Tp[1])**2)
T[2,0:types] = 1./(4*T[0,0:types]**3)*(q*Ap[2] + 4*Tp[0]**3*Tp[2] + 6*(Tp[0]*Tp[1])**2 - 6*(T[0,0:types]*T[1,0:types])**2)
beta[2,0:types] = rho*beta[0,0:types]*(2*(kappa - T[0,0:types])*T[2,0:types] + 1./(2*rho)*(beta[1,0:types]/beta[0,0:types])**2 - T[1,0:types]**2)
a[3,0:types] = 1./3*(a[0,0:types]*(x[0]*beta[2,0:types] + x[1]*beta[1,0:types] + x[2]*beta[0,0:types]) + a[1,0:types]*(x[0]*beta[1,0:types] + x[1]*beta[0,0:types]) + a[2,0:types]*(x[0]*beta[0,0:types] - gamma))
a[0] = a[0] + a[1]*t + a[2]*t**2 + a[3]*t**3
a_steady.append(a[0].tolist())
Tp_steady.append(Tp[0] - 273)
a_steady = sp.array(a_steady).transpose()
L_range = sp.array(L_range)
# Make a list of colors cycling through the rgbcmyk series.
colors = [colorConverter.to_rgba(c) for c in ('r','g','b','c','y','m','k')]
# plot the population values and the steady state temperature
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
line_daisys = collections.LineCollection([zip(L_range, daisy) for daisy in a_steady])
line_daisys.set_array(L_range)
line_daisys.set_color(colors)
ax.add_collection(line_daisys, autolim=True)
ax.set_title('Population, Planetary Temperature vs Luminosity', fontsize=16)
ax.set_xlabel('Solar Luminosity (\%)')
ax.set_ylabel('Temperature ($^\circ$C)')
# ax = plt.twinx()
# line_temperature = plt.plot(L_range, Tp_steady, color='green')
# ax.set_ylabel('Population Coverage (\%)')
# ax.yaxis.tick_right()
# plt.legend((line_temperature, line_daisys), ('Tp', 'a'))
plt.show()
return a_steady, Tp_steady, L_range, colors
if __name__ == '__main__':
a, Tp, L, colors = main()
-------------------------------------------------------------------------
This SF.net email is sponsored by the 2008 JavaOne(SM) Conference
Don't miss this year's exciting event. There's still time to save $100.
Use priority code J8TL2D2.
http://ad.doubleclick.net/clk;198757673;13503038;p?http://java.sun.com/javaone
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
