Python Getting Last evalues of a loop-CodePudding

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

CodePudding user response：

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