2017-01-14 53 views
-2

我想把這個Mathematica代碼改寫成Python,但是我很困惑,那麼請幫幫我!非常感謝。也許整合的結果與陣列不同,但是,我不知道!python plot整合結果,如何做像Mathematica那樣做

IMAGE: mathematica code which I want to rewrite in python

我用那麼Python代碼困擾....

# -*- coding: utf-8 -*- 
import numpy as np 
import matplotlib.pyplot as plt 
import matplotlib as mpl 
import math 
import scipy as sp 
from pylab import * 
#from scipy import * 
from scipy.integrate import quad, dblquad, tplquad 

plt.figure('God Bless : fig1') 
plt.title(r'fig1_') 
plt.xlabel(r'z') 
plt.ylabel('$L_0$(erg $s^{-1}$)') 
#plt.axis([0,10,-2,2]) 


Limit = 1.e48 
c = 2.997e10 
Om = 0.27 
H0 = 70e5/(1.e6*3.86e18) 
z = np.arange(0,10, 0.01) 
def dLz(z): 
    return 1./(1.-Om + Om*(1.+z)**3.)**(1./2) 
val, abserr = quad(dLz, 0, 10) 
print ("integral value =", val, ", absolute error =", abserr)  
dL  = val 
Fmin = 2.0e-8 
Llimit = 4.0*np.pi*dL**2*Fmin 
plt.plot(z, Limit) 


plt.savefig('fig_1.eps', dpi=300) 
plt.show() 
+0

也許DEF DLZ(Z):返回1/math.sqrt(1 +嗡* Z *(3 + Z *(3 + Z))) – Bill

回答

0

我認爲這就是它!哈利路亞!

#----------- DETERMIN zimax ------------- 

Limit = 1.e48 
c = 2.997e10 
Om = 0.27 
H0 = 70e5/(1.e6*3.86e18) 
Fmin = 2.0e-8 
z = np.arange(0,10, 0.1) 
def dLz(z): 
    return (c/H0)*(1.+z)/(1.-Om + Om*(1.+z)**3.)**(1./2) 
x_lower = 0 
vals = [] 
Llimits = [] 
for x_upper in z : 
    val, abserr = quad(dLz, x_lower, x_upper) 
    vals.append(val) 
    print ("integral value =", val, ", absolute error =", abserr) 
    dL  = val 
    Llimit = 4.0*np.pi*dL**2*Fmin 
    Llimits.append(Llimit) # add to array 
plt.plot(z, Llimits, 'b--') 

# ------------------------------------------ 
相關問題