I have a definition and it repeats with a while loop. I want to pull the penultimate value of m and ar from the operations that occur in the definition in each loop.
To do this, I opened an out file and tried to print the second-to-last element in each while loop.
But when I look at this file as output, there is only 1 line. How can I fix this?
import numpy as np
import matplotlib.pyplot as plt
def main(rho_c):
f = open("radial.out","w")
### Constants
pi = 3.1415926535897
gamma =5./3.
m_sun = 2.998e33
K = 1e10
P_c = K*rho_c**gamma
km = 1e5 # cm
dr = 1e4 # cm
r = 0;
Rho = np.zeros(2)
Rho[0] = rho_c
Rho[1] = rho_c
ar = np.zeros(2) # Tum r degerleri array olarak buna atilacak
ar[0] = r
ar[1] = dr
m = np.zeros(2)
m_k1 = 0
m_k2 = 0
m_k3 = 0
m_k4 = 0
P = np.zeros(2)
P_k1 = 0
P_k2 = 0
P_k3 = 0
P_k4 = 0
m[0] = 0
P[0] = P_c
P[1] = P[0]
m[1] = 4*pi*dr**3*rho_c/3
def der_m(r,m,ro):
return 4.*np.pi*(r**2)*ro
def ro(P):
if P==P_c:
return rho_c
return (P/K)**(3./5.)
def der_P(r,P,ro,m):
G = 6.67430e-8 # Gravitational constant cm^3/(g*s)^
c = 2.998e+10 # cm/s
return - (G*m*ro/(r**2))*(1.+(4.*np.pi*(r**3)*P/(m*c**2) ) )* \
(1. + P/(ro*c**2) )/(1-(2.*G*m/(r*c**2)))
i = 1
r = r+i*dr
while P[-1] > 0:
# For m
""" r , m , ro """
m_k1 = der_m( r , m[i] , ro(P[i]) )
m_k2 = der_m( r+0.5*dr , m[i]+dr*m_k1*0.5 , ro(P[i]) )
m_k3 = der_m( r+0.5*dr , m[i]+dr*m_k2*0.5 , ro(P[i]) )
m_k4 = der_m( r+1.*dr , m[i]+dr*m_k3 , ro(P[i]) )
m = np.append(m, m[i] + (dr/6.)*( m_k1+2.*m_k2+2.*m_k3+m_k4 ) )
# For P
P_k1 = der_P(r , P[i] ,ro(P[i]) , m[i])
P_k2 = der_P(r+0.5*dr , P[i]+dr*P_k1*0.5 ,ro(P[i]+dr*P_k1*0.5) , m[i]+dr*m_k1*0.5)
P_k3 = der_P(r+0.5*dr , P[i]+dr*P_k2*0.5 ,ro(P[i]+dr*P_k2*0.5) , m[i]+dr*m_k1*0.5)
P_k4 = der_P(r+1.*dr , P[i]+dr*P_k3 ,ro(P[i]+dr*P_k3) , m[i])
P = np.append(P, P[i] + (dr/6.)*( P_k1+2.*P_k2+2.*P_k3+P_k4 ) )
r = r+dr
ar = np.append(ar,r)
Rho = np.append(Rho, ro(P[i]) )
i = i+1
Hm = m[i]/np.abs(der_m(r,m[i],ro(P[i])))
Hp = P[i]/np.abs(der_P(r,P[i],ro(P[i]),m[i]))
H = (Hm*Hp)/(Hm+Hp)
dr = H
data = ar[-2]/km, m[-2]/m_sun
f.write(str(data)+'\n')
np.savetxt("radial.out", data, delimiter=" ")
plt.plot(ar/km,m/m_sun)
plt.xlim([0, 13])
plt.ylim([0,1.035])
plt.xlabel(r'$R$ (km)')
plt.ylabel(r'$M/M_\odot$')
#plt.savefig('plot_mass.pdf')
f.close()
rho_c = 2e15 # Density centeral kg/me3
while rho_c < 2e16:
#logspace
main(rho_c)
rho_c = rho_c*1.2
>Solution :
You open a new out file every time main() is called.
import numpy as np
import matplotlib.pyplot as plt
f = open("radial.out","w") # outside main()
def main(rho_c):
### Constants
pi = 3.1415926535897
gamma =5./3.
m_sun = 2.998e33
K = 1e10
P_c = K*rho_c**gamma
km = 1e5 # cm
dr = 1e4 # cm
r = 0;
Rho = np.zeros(2)
Rho[0] = rho_c
Rho[1] = rho_c
ar = np.zeros(2) # Tum r degerleri array olarak buna atilacak
ar[0] = r
ar[1] = dr
m = np.zeros(2)
m_k1 = 0
m_k2 = 0
m_k3 = 0
m_k4 = 0
P = np.zeros(2)
P_k1 = 0
P_k2 = 0
P_k3 = 0
P_k4 = 0
m[0] = 0
P[0] = P_c
P[1] = P[0]
m[1] = 4*pi*dr**3*rho_c/3
def der_m(r,m,ro):
return 4.*np.pi*(r**2)*ro
def ro(P):
if P==P_c:
return rho_c
return (P/K)**(3./5.)
def der_P(r,P,ro,m):
G = 6.67430e-8 # Gravitational constant cm^3/(g*s)^
c = 2.998e+10 # cm/s
return - (G*m*ro/(r**2))*(1.+(4.*np.pi*(r**3)*P/(m*c**2) ) )* \
(1. + P/(ro*c**2) )/(1-(2.*G*m/(r*c**2)))
i = 1
r = r+i*dr
while P[-1] > 0:
# For m
""" r , m , ro """
m_k1 = der_m( r , m[i] , ro(P[i]) )
m_k2 = der_m( r+0.5*dr , m[i]+dr*m_k1*0.5 , ro(P[i]) )
m_k3 = der_m( r+0.5*dr , m[i]+dr*m_k2*0.5 , ro(P[i]) )
m_k4 = der_m( r+1.*dr , m[i]+dr*m_k3 , ro(P[i]) )
m = np.append(m, m[i] + (dr/6.)*( m_k1+2.*m_k2+2.*m_k3+m_k4 ) )
# For P
P_k1 = der_P(r , P[i] ,ro(P[i]) , m[i])
P_k2 = der_P(r+0.5*dr , P[i]+dr*P_k1*0.5 ,ro(P[i]+dr*P_k1*0.5) , m[i]+dr*m_k1*0.5)
P_k3 = der_P(r+0.5*dr , P[i]+dr*P_k2*0.5 ,ro(P[i]+dr*P_k2*0.5) , m[i]+dr*m_k1*0.5)
P_k4 = der_P(r+1.*dr , P[i]+dr*P_k3 ,ro(P[i]+dr*P_k3) , m[i])
P = np.append(P, P[i] + (dr/6.)*( P_k1+2.*P_k2+2.*P_k3+P_k4 ) )
r = r+dr
ar = np.append(ar,r)
Rho = np.append(Rho, ro(P[i]) )
i = i+1
Hm = m[i]/np.abs(der_m(r,m[i],ro(P[i])))
Hp = P[i]/np.abs(der_P(r,P[i],ro(P[i]),m[i]))
H = (Hm*Hp)/(Hm+Hp)
dr = H
data = ar[-2]/km, m[-2]/m_sun
f.write(str(data)+'\n')
np.savetxt("radial.out", data, delimiter=" ")
plt.plot(ar/km,m/m_sun)
plt.xlim([0, 13])
plt.ylim([0,1.035])
plt.xlabel(r'$R$ (km)')
plt.ylabel(r'$M/M_\odot$')
#plt.savefig('plot_mass.pdf')
rho_c = 2e15 # Density centeral kg/me3
while rho_c < 2e16:
#logspace
main(rho_c)
rho_c = rho_c*1.2
f.close() # outside main and after while