def main(): T_night = 500.0 T_ss = 1000.0 gamma = 13.0 T_rad_eq_formula(T_night, T_ss, gamma) def radians(degrees): from math import pi return degrees * pi / 180.0 def T_rad_eq_formula(T_night, T_ss, gamma): from numarray import arange, cos, where # Define longitude array and locate where alpha1 = 90-gamma/2 is in it two = 2.0; four = two*two; fourth = 1.0/four alpha1, alpha2 = [90 - gamma/two, 90 + gamma/two] lon0 = arange(0.0, 90.25, 0.25) lon1 = arange(alpha1, alpha2+0.25, 0.25) lon0_radians = radians(lon0) lon1_radians = radians(lon1) i1 = where(lon0 == alpha1) # Define T0 function T0 = T_night**four + (T_ss**four - T_night**four)*cos(lon0_radians) T0 = T0**fourth T0_alpha1 = T0[i1] # Define T1 function T1 = (T_night - T0_alpha1)/gamma * (lon1 - alpha2) + T_night # Plot T_eq on the layer from pylab import plot, xlabel, ylabel, xlim, ylim, show, legend, gca, setp plot(lon0, T0, 'k-', lon1, T1, 'k--', linewidth = 2) xlabel(r'$\alpha$', fontsize = 20) ylabel(r'$\rm{T}_{\rm{eq}}(\alpha)$', fontsize = 20) leg = legend([r'$\rm{T}_{\rm{eq}}^0$', r'$\rm{T}_{\rm{eq}}^1$'], handlelen = 0.125, shadow = True) ltext = leg.get_texts() setp(ltext, fontsize = 16) show() if __name__ == '__main__': main()