2017-01-14 6 views
-2

このMathematicaコードをPythonに書き直したいのですが、困惑しました。どうもありがとう。おそらく、統合の結果は配列と異なりますが、私は知りません!Pythonプロットの積分結果、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):return 1/math.sqrt(1 + Om * 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--') 

# ------------------------------------------ 
関連する問題