# Improving the script of a plotter for an ODE integrator to be compacted, readable, understandable but encapsulated

This question is the aftermath of this other in which I'm asking for recommendations improving the space of combinations generated from a set of parameters, bounds, and conditions (calculations involving time-expensive math).

In this case, assuming the latter is improved, I want to make this other part of the code more readable and aligned with Python principles. This is mainly a plotting task and the (physical) "conditions imposed" to store information are not really relevant, as the meaning is non programatic.

The code will assume you ran:

import sys
import warnings
warnings.simplefilter("ignore")

import time
import random

import numpy as np
from numpy import diff
import scipy.integrate as integrate
import scipy

import matplotlib.pyplot as plt
from IPython.display import clear_output

global kjup
kjup = 0.0009543

plt.rcParams.update({'font.size': 26, 'font.family' : 'Bitstream Vera Sans',
'font.weight' : 'normal'})

plt.rcParams['lines.linewidth'] = 3


and then

def Table_plus_Analytical(M,m1,m2,               #mean, lower ERROR, upper ERROR
alpha,                 #lowest VAL, highest VAL
ta1,ta2,
te1,te2,
running_time=5/60, #in minutes
mass_array_generation='arange', #random normal, random uniform, arange
params_array_generation='arange', #random uniform, arange
f1_res = -1.190,f2_res = 0.428,
q = 1,p = 3
):

start_time = time.time()

def Calculate_N_elements_per_array(M,m1,m2,alpha,ta1,ta2,te1,te2,running_time):

#parameters in a list; counting every parameter that will vary

iterated_parameters_count, parameters= 0, [M,m1,m2,alpha,ta1,ta2,te1,te2]

for par in parameters:
if isinstance(par, list)==True and len(par)>1:
iterated_parameters_count +=1
dim = float(iterated_parameters_count)

#estimating a linear response
aprox_slope = .01 #seconds/million of values

#Then time in seconds = aprox_slope * values (in millions)

running_time = 60*running_time #from minutes to seconds

num_of_values_per_arrange = (running_time/aprox_slope)**(1/dim)
N = num_of_values_per_arrange
#(N) = (t/m)**(1/dim) if Ntotal = (N)

print(str(int(N))+' elements in each parameter will be generated')
return int(N)

def Make_arrage(N, args, mass_array_generation, params_array_generation):

masses = args[:3]
elseargs = args[3:]

if mass_array_generation == 'random normal':
############### Convert masses to np.arrays if desirable (when mass is list)
mtemp = []

for mass in masses:
if isinstance(mass,int) or (isinstance(mass,list) and len(mass)==1):
mass = mass

if (isinstance(mass,list) and len(mass)>1):

mu = mass[0]
lo = mu - mass[1]
hi = mu + mass[2]

std = abs(hi-lo)/2
mass = np.random.normal(mu, std, size=N)

mtemp.append(mass)
masses = mtemp
###############

if mass_array_generation == 'random uniform':
############### Convert masses to np.arrays if desirable (when mass is list)
mtemp = []

for mass in masses:
if isinstance(mass,int) or (isinstance(mass,list) and len(mass)==1):
mass = mass

if (isinstance(mass,list) and len(mass)>1):

mu = mass[0]
lo = mu - mass[1]
hi = mu + mass[2]

std = abs(hi-lo)/2
mass = np.random.uniform(lo, hi, size=N)

mtemp.append(mass)
masses = mtemp

if mass_array_generation == 'arange':
############### Convert masses to np.arrays if desirable (when mass is list)
mtemp = []

for mass in masses:
if isinstance(mass,int) or (isinstance(mass,list) and len(mass)==1):
mass = mass

if (isinstance(mass,list) and len(mass)>1):

mu = mass[0]
lo = mu - mass[1]
hi = mu + mass[2]

std = abs(hi-lo)/2
mass = np.arange(lo, hi, (hi-lo)/N)

mtemp.append(mass)
masses = mtemp
###############

if params_array_generation == 'random uniform':
elseargs = [np.random.uniform(arg[0],arg[1],size=N) for arg in elseargs]

if params_array_generation == 'arange':
elseargs = [np.arange(arg[0], arg[1], abs(arg[0]-arg[1])/N) for arg in elseargs]

#Now masses are np.arrays and floats joint in a list.

array_args = masses+elseargs
return array_args

#################

q=1 #!#!

def Calculate_fs(q, alpha_array):

def db_1_2(psi, j, alph): #alpha is replaced with "alph" because it's already a variable for alpha limits
return (np.cos(j*psi))/((1-2*alph*np.cos(psi) + alph**2)**(1./2))

def db_3_2(psi, j, alph):
return (np.cos(j*psi))/((1-2*alph*np.cos(psi) + alph**2)**(3./2))

def b_1_2(j, alph):
return (1./np.pi) * integrate.quad(db_1_2, 0., 2*np.pi, args=(j, alph))[0]

def b_3_2(j, alph):
return (1./np.pi) * integrate.quad(db_3_2, 0., 2*np.pi, args=(j, alph))[0]

def Calc_f1(q, alph):
f1_result = (-1./2) * (2*(q+1)*b_1_2(q+1.,alph)
+ (alph/2) * ( b_3_2(q+2., alph)  +  b_3_2(q, alph)) - (alph**2) * b_3_2(q+1., alph))
return f1_result

def Calc_f2(q, alph):
f2_result = (1./2) * ((2.*q+1.) * b_1_2(q,alph)
+ (alph/2.) *(b_3_2(q+1., alph)  + b_3_2(q-1., alph)) - (alph**2) * b_3_2(q, alph))
if q!= 1.:
return f2_result
else:
return f2_result - 2*alph

f1_array=[]
f2_array=[]

for alpha_val in alpha_array:
f1_array.append(Calc_f1(q, alpha_val))
f2_array.append(Calc_f2(q, alpha_val))

return np.array(f1_array), np.array(f2_array)

def Calc_Delta_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2):

def Calc_A_factor(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2):

fact_1 = (3/(q**2. * te1))
fact_2 = (  (q/(q+1)) * (m1/alpha/m2) + 1 )
fact_3 = (alpha*m2/Ms)**2.
fact_4 = (  (q+1)  *  f1**2.  +  (q**2.  /  (q+1)  )  *  (m1/alpha/m2)  *  f2**2.  *  (te1  /  te2)  )

A = fact_1 * fact_2 * fact_3 * fact_4

return A

def Calc_B_factor(ta1,ta2):

fact_1 = 3/(2* ta1)
fact_2 = (1 - (ta1/ta2))

B = fact_1*fact_2
return B

A = Calc_A_factor(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
B = Calc_B_factor(ta1,ta2)

Delta_eq = ( -1 * A / B)**(1./2.)
return Delta_eq

def Calc_e1_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2):
Delta = Calc_Delta_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
e1= alpha*(m2/Ms) * (abs(f1)/(Delta *q))
return e1

def Calc_e2_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2):
Delta = Calc_Delta_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
e2 = (m1/Ms) * (abs(f2)/((q+1)*Delta))
return e2

def Calc_dta(ta1,ta2):
return ta1/ta2

def Calc_taitei(ta,te):
return ta/te

#################

def Calculate_Filter_Possibilites(data, f1_res=f1_res, f2_res=f2_res, q=q, p=p):

plug_in = np.stack(np.meshgrid(*data), axis=-1).reshape(-1, len(data)) #Cartesian product

"""
calculated_Delta_eq=Calc_Delta_eq(Ms,m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
calculated_e1eq = Calc_e1_eq(m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
calculated_e2eq = Calc_e2_eq(m1,m2,alpha,ta1,ta2,te1,te2,q,f1,f2)
"""
len_plugin = len(plug_in)

old_len= len(plug_in)

calculated_f1, calculated_f2 = Calculate_fs(q,plug_in[:,3])

q =q *np.ones(len_plugin)
p =p *np.ones(len_plugin)

plug_in = np.column_stack([plug_in,q[:,None], calculated_f1[:,None],calculated_f2[:,None],p[:,None]]) #append extra columns

calculated_Delta_eq=Calc_Delta_eq(plug_in[:,0],plug_in[:,1]*kjup,plug_in[:,2]*kjup,plug_in[:,3],plug_in[:,4],plug_in[:,5],plug_in[:,6],plug_in[:,7],plug_in[:,8],abs(plug_in[:,9]),abs(plug_in[:,10]))
calculated_e1eq = Calc_e1_eq(plug_in[:,0],plug_in[:,1]*kjup,plug_in[:,2]*kjup,plug_in[:,3],plug_in[:,4],plug_in[:,5],plug_in[:,6],plug_in[:,7],plug_in[:,8],abs(plug_in[:,9]),abs(plug_in[:,10]))
calculated_e2eq = Calc_e2_eq(plug_in[:,0],plug_in[:,1]*kjup,plug_in[:,2]*kjup,plug_in[:,3],plug_in[:,4],plug_in[:,5],plug_in[:,6],plug_in[:,7],plug_in[:,8],abs(plug_in[:,9]),abs(plug_in[:,10]))

calculated_dta = Calc_dta(plug_in[:,4],plug_in[:,5])

#,calculated_Delta_eq[:,None], calculated_e1eq[:,None], calculated_e2eq[:,None], calculated_dta[:,None], calc_mult[:,None]

calculated_ta1te1 = Calc_taitei(plug_in[:,4],plug_in[:,6])
calculated_ta2te2 = Calc_taitei(plug_in[:,5],plug_in[:,7])

plug_in = np.column_stack([plug_in, calculated_Delta_eq[:,None], calculated_e1eq[:,None], calculated_e2eq[:,None], calculated_dta[:,None], calculated_ta1te1[:,None], calculated_ta2te2[:,None]])

#Free on Delta, constrained corrected on eccentricity of Planet 1 (that is constraining Delta and e2 implicitly)
plug_in = plug_in[(0.034 < calculated_Delta_eq) & (calculated_Delta_eq<0.037)& (calculated_dta >= 1)] #drop rows you don't need
#& (0.01 < calculated_e1eq) & (calculated_e1eq < 0.07) & (0.04 < calculated_e2eq) & (calculated_e2eq < 0.17)
new_len = len(plug_in)

print(str((new_len/old_len)*100)[:4]+'% of combinations satisfy the conditions; from '+str(old_len)+' possibilities, now there are '+str(new_len)+' possibilities.')
print('Conditions:\n(0.03 < calculated_Delta_eq) & (calculated_Delta_eq<0.13)& (0.042/5 < calculated_e1eq) & (calculated_e1eq < 0.042/2) & (calculated_dta >= 1)')

return plug_in

####################

N = Calculate_N_elements_per_array(M,m1,m2,alpha,ta1,ta2,te1,te2,running_time)
if N<1:
N=1

array_args = Make_arrage(N, [M,m1,m2,alpha,ta1,ta2,te1,te2], mass_array_generation, params_array_generation)

Table_filtered_possibilities = Calculate_Filter_Possibilites(array_args)

return Table_filtered_possibilities


for fun, you can check:


result_2 = Table_plus_Analytical([0.823,0.041,0.041], m1 = [0.978,0.06,0.06],m2 = [0.369,0.1,0.08],
alpha = [0.64,0.68],
ta1 = [1e3,3e5],ta2 = [1e2,1e5],
te1 = [1e1,3e3],te2 = [1e0,1e3],
running_time=11,
mass_array_generation='random normal', #random normal, random uniform, arange
params_array_generation='random uniform', #random uniform, arange
f1_res = 0,f2_res = 0,
q = 1,p = 3
)

with np.printoptions(precision=4):
print(result_2)



and then

Nbins=20
height=100

plt.figure(figsize=(10,7))
plt.locator_params(axis='x', nbins=4)
plt.hist(result_2.T[12],bins=Nbins)
plt.scatter(0.035,0,alpha=0, label=r'$N_{bins}=$'+str(Nbins))
plt.legend()

plt.vlines(0.0359,0,height, color='blue', linestyle='--')

plt.ylabel(r'count')

plt.xlabel(r'Analytical eq. $\Delta$')

plt.grid(True, alpha=0.1, linewidth=3, color='blue', linestyle='--')
plt.tight_layout()
plt.show()

plt.figure(figsize=(10,7))
#plt.xlim([0.02,0.07])
#plt.ylim([0,1e3])

plt.ylabel(r'count')

plt.hist(result_2.T[13],bins=Nbins, color='blue', alpha=0.3)
plt.hist(result_2.T[14],bins=Nbins, color='orange', alpha=0.3)
plt.scatter(0,0,alpha=0)
plt.legend([r'$N_{bins}=$'+str(Nbins),r'$e_{1eq}$',r'$e_{2eq}$'])
plt.xlabel(r'Analytical eq. eccentricities $e_{ieq}$')

plt.vlines(0.042,0,height, color='blue', linestyle='--')
plt.vlines(0.0622,0,height, color='orange', linestyle='--')
#plt.yticks([30,60,90])
plt.tight_layout()
plt.grid(True, alpha=0.1, linewidth=3, color='blue', linestyle='--')



and you'll see what this first fragment does:

My first question: how can I display the number of bins in a better way?

Moving now to more interesting lines of code: this will generate plots from a differential equations set that is also solved. The script is

global frsizes #figures_resize_factor
frsizes=1.9

def plot_deltas(N_times, t, n1, n2, svd_deltaeq, tolerance = 0.1):

abort=0

p2_p1=np.zeros((N_times,1))
for j in range(N_times):
p2_p1[j]=n1[j]/n2[j]

Delta_array= (p2_p1)/2 -1
minval=np.min(Delta_array)
maxval=np.max(Delta_array)

"""
if np.isnan(minval)==True or np.isinf(minval)==True or np.isinf(maxval)==True or np.mean(abs(Delta_array))>1:
print(r'the abs. mean value of $\Delta$ exceeds 1'+' or NaN or +/- inf vals. were encountered\n')
abort=1
return np.nan, np.nan, abort
"""
dot_Delta = diff(Delta_array.reshape(-1))/diff(t)

plt.figure(figsize=(17*frsizes,3.5*frsizes))

plt.plot(t, Delta_array)
plt.plot(t[:-1], dot_Delta*10**2.)
plt.hlines(0.0359,t[0],t[-1], color='green', linestyles='dashed')
plt.hlines(svd_deltaeq,t[0],t[-1], color='blue', linestyles='dashed')
plt.legend([r'$\Delta = \Delta(t) \; | \; \Delta(t_f) =$'+str(round(Delta_array[-1][0],4)),  r'$\frac{d\Delta}{dt} (t) \; | \; \frac{d\Delta}{dt} (t_f)\rightarrow 0$',  r'$\Delta_{obs}^{(Trif.21)} = 0.0359^{+0.0011}_{-0.0015}$',  r'$\Delta_{eq_{theo}}=$'+str(round(svd_deltaeq,4))], markerscale =2., fontsize = 20, loc = 'upper right')

plt.title(r'Deviation $\Delta$ from $q=1$ resonance')
plt.ylabel(r'$\Delta$ and d$\Delta$/dt $\times 10^2$')
plt.xlabel('Time $t$ [yr]')
plt.xlim(t[0],t[-1])

plt.grid(True, ls = '--', c = 'b', linewidth = .3)

def Candidate_conditions(meandotDel, Delta_at_t, theo_deltaeq):
if abs(meandotDel) < 1e-6 and 0.034 < Delta_at_t< 0.037:
return True
else:
return False

if Candidate_conditions(np.mean(dot_Delta[-100]), Delta_array[-1][0], svd_deltaeq) == True:
pass
#plt.scatter(int(2e4),0.06,color='green',marker='$\mathrm{candidate\; with\; }\dot{\Delta}\, \sim\, 0$', s=int(2e5))

plt.tight_layout()
plt.show()

print(r'$\Delta_{t_f}$'+' is near to: '+str(Delta_array[-1][0]))
print('mean(dDelta/dt)[-100] is near to: '+str(np.mean(dot_Delta[-100]))+'\n')

return Delta_array[-1][0], dot_Delta, abort

def plot_e(N_times, e1, e2, t, svd_e1eq, svd_e2eq):

plt.figure(figsize=(17*frsizes,7*frsizes))

plt.subplot(2, 1, 1)
plt.title(r'Eccentricities $e_1$, $e_2$ variation in time')

plt.plot(t, e1)
plt.plot(t, e2)

plt.ylabel(r'$e_i$')
plt.xlabel(r'Time $t$ [yr]')
plt.legend([r'$e_{1eq}=$ '+str(e1[-1])[:7]+'; '+r'$e_{1eq_{theo}}=$ '+str(svd_e1eq)[:7]+'; '+r'$e_{1obs}^{\mathrm{(Trif.21)}} = 0.0420^{+0.0255}_{-0.0075}$',r'$e_{2eq}=$ '+str(e2[-1])[:7]+'; '+r'$e_{2eq_{theo}}=$ '+str(svd_e2eq)[:7]+'; '+r'$e_{2obs}^{\mathrm{(Trif.21)}} = 0.0622^{+0.0452}_{-0.0211}$'], markerscale =2., fontsize = 20)
plt.grid(True, ls = '--', c = 'b', linewidth = .3)

"""
clip = 200

plt.subplot(2, 1, 2)

plt.plot(t[-clip:], e1[-clip:])
plt.plot(t[-clip:], e2[-clip:])

plt.ylabel(r'$e_i$')
plt.xlabel(r'Time $t$ [yr]')
plt.legend([r'$e_{1eq}=$ '+str(e1[-1])[:7],r'$e_{2eq}=$ '+str(e2[-1])[:7]])
plt.grid(True, ls = '--', c = 'b', linewidth = .3)
"""

plt.show()
return [e1[-1], e2[-1]]

def plot_smas(t, n1, n2):

plt.figure(figsize=(17*frsizes,3.5*frsizes))

plt.plot(t, (2*np.pi/np.array(n1))**(2/3))
plt.plot(t, (2*np.pi/np.array(n2))**(2/3))

plt.title(r'$\propto$ semi-major axes $a_i$ in time')
plt.ylabel(r'$\left(2\pi \: n_i^{-1}\right)^{2/3}$')
plt.xlabel('Time $t$ [yr]')
plt.legend([r'$n_1$',r'$n_2$'],markerscale =2., fontsize = 15)
plt.grid(True, ls = '--', c = 'b', linewidth = .3)

plt.show()

return [0,0]

def plot_Pandes(t, e1, e2, n1, n2, N_times):

plt.figure(figsize=(17*frsizes,3.5*frsizes))

p2_p1=np.zeros((N_times,1))
for j in range(N_times):
p2_p1[j]=n1[j]/n2[j]

plt.plot(e1, p2_p1)
plt.plot(e2, p2_p1)

plt.title(r'Eccentricity change in $P_2/P_1$')
plt.ylabel(r'$\left(\frac{n_1}{n_2}\right)$')
plt.xlabel('Eccentricities $e_i$')
plt.legend([r'$e_1$',r'$e_2$'], markerscale =2., fontsize = 15)
plt.grid(True, ls = '--', c = 'b', linewidth = .3)
plt.ylim([2.04-0.3, 2.34])

plt.show()

plt.figure(figsize=(17*frsizes,3.5*frsizes))

plt.plot(t, p2_p1)

plt.title(r'Period ratios $P_2/P_1$')
plt.ylabel(r'$\left(\frac{n_1}{n_2}\right)$')
plt.xlabel('Time $t$ [yr]')
plt.grid(True, ls = '--', c = 'b', linewidth = .3)

plt.show()

def plot_phi_t(t, phi1, phi2):

plt.figure(figsize=(17*frsizes,10.5*frsizes))

plt.title(r'Resonance angles $\phi_1$, $\phi_2$ evolution in time')

plt.subplot(3, 1, 1)
plt.semilogx(t, phi2)
plt.semilogx(t, phi1)

plt.ylabel(r'$\phi_i$  [°]')
#plt.xlabel(r'Time $log(t)$ [yr]')
plt.legend([r'$\phi_1$',r'$\phi_2$'], markerscale =2., fontsize = 15)
plt.grid(True, ls = '--', c = 'b', linewidth = .3)

plt.subplot(3, 1, 2)

plt.plot(t, np.cos(np.array(phi1)))
plt.plot(t, np.cos(np.array(phi2)))

plt.ylabel(r'$\cos{\phi_i}$')
plt.xlabel(r'Time $t$ [yr]')
plt.legend([r'$\phi_1$',r'$\phi_2$'], markerscale =2., fontsize = 15)
plt.grid(True, ls = '--', c = 'b', linewidth = .3)

"""
plt.subplot(3, 1, 3)

plt.semilogx(t, [i % 180 for i in phi2])
plt.semilogx(t, [i % 180 for i in phi1])

plt.ylabel(r'wrap. $0° < \phi_i < 180°$')
#plt.xlabel(r'Time $log(t)$ [yr]')
plt.legend([r'$\phi_1$',r'$\phi_2$'])
plt.grid(True, ls = '--', c = 'b', linewidth = .3)

"""
plt.show()

#Values of --- Planets b and c are extracted from https://arxiv.org
#Ms in Solar masses
#m1 and m2 in Jup. masses
#timescales in years

#IN MJUPS

def tester(svd_Ms, svd_m1, svd_m2, svd_alpha, svd_ta1, svd_ta2, svd_te1, svd_te2, svd_deltaeq, svd_e1eq, svd_e2eq, test_a1=3, test_e10 = 0.1, test_e20 = 0.1, test_phi10=10, test_phi20=0, plot_detailed = 0, test_tmax= int(3e4), test_Ntimes = int(1e4)):
warnings.simplefilter("ignore")

# we fix f1 and f2 to speed up the code, and only evolve it close to the commensurability
q = 1. ## Near q+1:q mean motion resonance
p = 3. ## Conservation of angular momentum

alpha_res=(q/(q+1))**(2/3)

def db_1_2(psi, j, alpha):
return (np.cos(j*psi))/((1-2*alpha*np.cos(psi) + alpha**2)**(1./2))

def db_3_2(psi, j, alpha):
return (np.cos(j*psi))/((1-2*alpha*np.cos(psi) + alpha**2)**(3./2))

def b_1_2(j, alpha):
return (1./np.pi) * integrate.quad(db_1_2, 0., 2*np.pi, args=(j, alpha))[0]

def b_3_2(j, alpha):
return (1./np.pi) * integrate.quad(db_3_2, 0., 2*np.pi, args=(j, alpha))[0]

def f1(q, alpha):
f1 = (-1./2) * (2*(q+1)*b_1_2(q+1.,alpha)
+ (alpha/2) * ( b_3_2(q+2., alpha)  +  b_3_2(q, alpha)) - (alpha**2) * b_3_2(q+1., alpha))
return f1

def f2(q, alpha):
f2 = (1./2) * ((2.*q+1.) * b_1_2(q,alpha)
+ (alpha/2.) *(b_3_2(q+1., alpha)  + b_3_2(q-1., alpha)) - (alpha**2) * b_3_2(q, alpha))
if q!= 1.:
return f2
else:
return f2 - 2*alpha

a1 = test_a1
a2 = a1 / alpha

f1_res=f1(q,alpha_res)
f2_res=f2(q,alpha_res)

def equations(y, t, q, alpha, m1, m2, ms, p, ta1, ta2, te1, te2):
n1, n2, e1, e2, phi1, phi2 = y

dn1dt = -3 * q * (n1 ** 2) * alpha * (m2 / ms) * (e1 * f1_res * np.sin(phi1)
+ e2 * f2_res * np.sin(phi2)) + 3 * (n1 / (2 * ta1)) + (p * n1 * (e1 ** 2)) / (te1)

dn2dt = 3 * (q + 1) * (n2 ** 2) * (m1 / ms) * (e1 * f1_res * np.sin(phi1)
+ e2 * f2_res * np.sin(phi2))+ 3 * (n2 / (2 * ta2)) + (p * n2 * (e2 ** 2)) / (te2)

de1dt = -n1 * alpha * (m2 / ms) * f1_res * np.sin(phi1) - (e1 / te1)
de2dt = -n2 * (m1 / ms) * f2_res * np.sin(phi2) - (e2 / te2)

dphi1dt = (q + 1) * n2 - q * n1 - n1 * alpha * (m2 / ms) * (1 / e1) * f1_res * np.cos(phi1)
dphi2dt = (q + 1) * n2 - q * n1 - n2 * (m1 / ms) * (1 / e2) * f2_res * np.cos(phi2)

dydt = [dn1dt, dn2dt, de1dt, de2dt, dphi1dt, dphi2dt]
return dydt

kjup = 0.0009543 #1 Jupiter mass in Solar Mass.
G = 39.476926 #AU^3/M⊙*yr^2

#ms = 0.823 ## Solar Mass Star (Trifonov et al. 2021 Table 1.)
#m1 = 0.978*kjup #Inner planet (b) from Trifonov et al. 2021 Table 4.
#m2 = 0.369*kjup #Outer planet (c) from Trifonov et al. 2021 Table 4.

ms = svd_Ms
m1 = svd_m1*kjup
m2 = svd_m2*kjup

P1 = np.sqrt(((4*np.pi**2)/(ms * G) )*a1**3 ) ## Years
P2 = np.sqrt(((4*np.pi**2)/(ms * G) )*a2**3 ) ## Years
n1_0 = (2*np.pi)/P1
n2_0 = (2*np.pi)/P2

ta1 = svd_ta1
ta2 = svd_ta2
te1 = svd_te1
te2 = svd_te2        #te2 =  3.e2*(q/(q+1))**(2/3) ## P1/P2 * tdamp

e1_0 = test_e10
e2_0 = test_e20

phi1_0 = test_phi10
phi2_0 = test_phi20

tmax = test_tmax #e3-e4
N_times = test_Ntimes #e3-1.e4 times

y0 = [n1_0, n2_0, e1_0, e2_0, phi1_0, phi2_0]
t = np.linspace(0, tmax, N_times)

sol = integrate.odeint(equations, y0, t, args=(q, alpha, m1, m2, ms, p, ta1, ta2, te1, te2))

n1 = [i[0] for i in sol]
n2 = [i[1] for i in sol]
e1 = [i[2] for i in sol]
e2 = [i[3] for i in sol]
phi1 = [i[4] for i in sol]
phi2 = [i[5] for i in sol]

#plt.rcParams.update({'font.size': 26})
#plt.rcParams['lines.linewidth'] = 3

Delta_final, dot_Delta, abort = plot_deltas(N_times, t, n1, n2, svd_deltaeq)

if np.isnan(np.min(dot_Delta)) == True:
mean_dot_Delta = np.nan
if np.isnan(np.min(dot_Delta)) == False:
mean_dot_Delta = np.mean(dot_Delta[-100:])

e1_final, e2_final, a1_final, a2_final = [0,0,0,0]

if plot_detailed == 1 and abort == 0:

e1_final, e2_final = plot_e(N_times, e1, e2, t, svd_e1eq, svd_e2eq)
#plot_Pandes(t, e1, e2, n1, n2, N_times)
plot_phi_t(t, phi1, phi2)
#a1_final, a2_final = plot_smas(t, n1, n2)

return np.array([Delta_final, mean_dot_Delta, svd_deltaeq, e1_final, svd_e1eq, e2_final, svd_e2eq, a1_final, a2_final])

#Delta final, mean(dot(Delta)), theo_Delta,
#e1 final, theo_e1,
#e2 final, theo_e2,
#a1_final, a2_final



and this code is also filtering solutions aborting (or not) a calculation if it satisfies certain conditions. In this last part I use both:

def combinate_shuffle_test(
N_results_evol,
M = [0.823,0.041,0.041], m1 = [0.978,0.06,0.06],m2 = [0.369,0.1,0.08],
alpha = [0.60,0.65],
ta1 = [1e3,1e6],ta2 = [1e2,1e6],
te1 = [1e1,1e4],te2 = [1e0,1e4],
running_time=2/6,
mass_array_generation='random normal', #random normal, random uniform, arange
params_array_generation='random uniform', #random uniform, arange
f1_res = -1.190,f2_res = 1.688,
q = 1,p = 3
):

result = Table_plus_Analytical(M,m1,m2,
alpha,
ta1,ta2,
te1,te2,
running_time,
mass_array_generation,
params_array_generation,
f1_res,f2_res,q,p)

"""
with np.printoptions(precision=4):
print(result)
"""

print('\n\nQWE\n\nResults:')

random.shuffle(result)

N_result=result[:N_results_evol]
c = 1
for combination in N_result:

print('\n\n test N°'+str(c))
svd_Ms, svd_m1, svd_m2, svd_alpha, svd_ta1, svd_ta2, svd_te1, svd_te2, q, f1, f2, p, svd_deltaeq, svd_e1eq, svd_e2eq, svd_dta, *args = combination

data = tester(svd_Ms, svd_m1, svd_m2, svd_alpha, svd_ta1, svd_ta2, svd_te1, svd_te2, svd_deltaeq, svd_e1eq, svd_e2eq, test_a1=2, test_e10 = 0.05, test_e20 = 0.05, test_phi10=10, test_phi20=0, plot_detailed = 1, test_tmax= int(3e4), test_Ntimes = int(1e3))
print('\n\n')
c+=1

#Filter!!!!!!!!!!!!!!!!!!!!!
filtered_indexes=[]
c=0
if abs(data[0]-0.0359)<0.01: #0.01 = 7% relative error
print('match encountered for the condition: abs(data[0]-0.0359)<0.0025 (~7%err)')
filtered_indexes.append(c)
#delta(f)
c+=1

filtered_configs, filtered_results_evolved= [], []

for index in filtered_indexes:
filtered_configs.append(np.array(N_result[index]))

return filtered_configs, filtered_results_evolved



The next cell will create a "master_results.npy" in which relevant data will be added; also will call combinate_shuffle_test and will save relevant results. If you run this

import os

#os.remove('master_results.npy')

master_configs = [[],[]]
master_evols = [[],[]]
#np.savetxt('master_results.txt', np.array([master_configs,master_evols]), fmt='%1.6f', newline='\n', delimiter=', ') #To save

while True:

clear_output(wait=True)

try:

master_configs, master_evols = call[0].tolist(), call[1].tolist() #To load
#os.remove
except:
pass

fil_cfgs, fil_res_evol = combinate_shuffle_test(200, running_time=11) # call = ...

#fil_cfgs = call[0]
#fil_res_evol = call[1]

master_configs += fil_cfgs
master_evols += fil_res_evol

#os.remove('master_results.npy') #
np.save('master_results.npy', np.array([master_configs,master_evols]), allow_pickle=True) #To save

if len(fil_cfgs)>0:

from inputimeout import inputimeout, TimeoutOccurred
try:
c = inputimeout(prompt=str(len(fil_cfgs))+' configs. appended; Press enter | 10 seconds timeout | You can load the configurations with np.load(master_results.npy) [0] and [1]\n', timeout=10)
except TimeoutOccurred:
c = 'Timeout!'
print(c)
print([result[0] for result in fil_res_evol])



you'll see:

• Combinations are generated and filtered. For instance:
6 elements in each parameter will be generated
2.80% of combinations satisfy the conditions; from 7776 possibilities, now there are 218 possibilities.
Conditions:
(0.03 < calculated_Delta_eq) & (calculated_Delta_eq<0.13)& (0.042/5 < calculated_e1eq) & (calculated_e1eq < 0.042/2) & (calculated_dta >= 1)

• Combinations are solved, filtered and (if) plotted. An example is:

• A fragment of combinations will be considered, shuffled and solved, filtered, (if) plotted, appending interesting ones.
• A loop through different shuffled-fragments of filtered-combinations, solving, plotting and saving some in master_results.npy.
• I can do fun reading masters, such as
[Df,mdD,teoD,e1f,teoe1,e2f,teoe2] = [[row[i] for row in np.load('master_results.npy',allow_pickle=True)[1][2:]] for i in [0,1,2,3,4,5,6]]
#Df
val=e1f
index_max = max(range(len(val)), key=val.__getitem__)
print(index_max)
plt.hist(val, bins=35)
plt.vlines(np.median(val),0,150, color='yellow',linestyle='--')
print(np.median(val))
print(np.max(val))


which is just a simple test.

This will also bring good information about differences on analytical and numerical calculations. See, e.g.,

[Df,mdD,teoD,e1f,teoe1,e2f,teoe2] = [[row[i] for row in np.load('master_results.npy',allow_pickle=True)[1][2:]] for i in [0,1,2,3,4,5,6]]
#3rd plot

plt.figure(figsize=(12,8))
plt.scatter(Df,teoD, s=40, alpha=0.6)
plt.xlabel(r"$\Delta(t_f)$ from Eqs. of Motion")
plt.ylabel(r"$\Delta_{\mathrm{eq}}^{\mathrm{(Ch)}}=\sqrt{-\mathcal{A}/\mathcal{B}}$")
plt.grid(True, color='blue', linewidth=3, linestyle='--', alpha=0.1)
plt.xlim([0.03, 0.04])
plt.ylim([0.03, 0.04])
xdum=np.linspace(0.03,0.04,100)
plt.plot(xdum,xdum,color='green', linestyle='--')
plt.show()


which outputs:

could you please bring suggestions to improve my lack of coding skills?

Best,

• Is this an expanded version of your other question? codereview.stackexchange.com/q/271017/232484 Dec 15, 2021 at 4:08
• @Teepeemm Thanks very much for this other comment. Is this an expanded question? Maybe, but also not actually. My other question is pursueing an improvement in the way I calculate and append coefficents and combinations of parameters. This other question assume that this is not a problem (because the time-management issue is assumed to be fixed) and is asking for general reviews and comments to the way I generate plots and write code.
– nuwe
Dec 15, 2021 at 12:26