diff --git a/UserSpace/Python/AnharmonicFactor.py b/UserSpace/Python/AnharmonicFactor.py index 2a71e35..6231fa8 100755 --- a/UserSpace/Python/AnharmonicFactor.py +++ b/UserSpace/Python/AnharmonicFactor.py @@ -1,36 +1,36 @@ from numpy import linspace, pi from sys import path as sysPath -from os import path as osPath -sysPath.append(osPath.join(osPath.dirname(__file__), '../../src')) +sysPath.append('../../src') + #load the module from interfacePy.AnharmonicFactor import anharmonicFactor #load the anharmonicFactor function form the AnharmonicFactor module for theta_maX in linspace(0,pi,50): print(theta_maX,anharmonicFactor(theta_maX)) # change to True to plot data alongsode the interpolation (if there are large differences there is something wrong) if False: import matplotlib.pyplot as plt fig=plt.figure(figsize=(9,4)) fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.99, right=0.9,wspace=0.0,hspace=0.25) sub = fig.add_subplot(1,1,1) Y=[] T=linspace(0,pi,50) for t in T: Y.append(anharmonicFactor(t)) sub.plot(T,Y,linestyle='-',linewidth=2,alpha=1,c='xkcd:black') sub.set_xlabel(r'$\theta_{\rm max}$') sub.xaxis.set_label_coords(0.5, -0.1) sub.set_ylabel(r'$ \dfrac{2\sqrt{2}}{\pi \ \theta_{\rm max}^2} \ \int_{-\theta_{\rm max}}^{\theta_{\rm max}} d\theta \ \sqrt{ \cos(\theta) - \cos(\theta_{\rm max}) }$') # sub.set_ylabel(r'$ f\left( \theta_{\rm max} \right)$') sub.yaxis.set_label_coords(-0.1,0.5) fig.savefig('anharmonic_factor_examplePlot.pdf',bbox_inches='tight') \ No newline at end of file diff --git a/UserSpace/Python/Axion.py b/UserSpace/Python/Axion.py index 11ad702..ebcc044 100755 --- a/UserSpace/Python/Axion.py +++ b/UserSpace/Python/Axion.py @@ -1,197 +1,193 @@ from time import time from sys import stderr from sys import path as sysPath -from os import path as osPath -sysPath.append(osPath.join(osPath.dirname(__file__), '../../src')) +sysPath.append('../../src') #load the module from interfacePy.Axion import Axion from interfacePy.AxionMass import AxionMass -from interfacePy.Cosmo import Cosmo +from interfacePy.Cosmo import mP -Hubble=Cosmo().Hubble -mP=Cosmo().mP - theta_i, fa=0.94435, 1e12 umax=500 TSTOP=1e-4 ratio_ini=1e3 N_convergence_max, convergence_lim=5, 1e-2 #this is fine, but you can experiment a bit. #radiation dominated example inputFile="../InputExamples/RDinput.dat" # Matter Domination example. # the NSC parameters (using the notation of 2012.07202) are: # T_end=1e-2 (GeV), c=3, T_ini=1e12 (GeV), and r=1e-1 # inputFile="../InputExamples/MatterInput.dat" # Kination Domination example. # the NSC parameters (using the notation of 2012.07202) are: # T_end=0, c=6, T_ini=1e3 (GeV), and r=1e10 # inputFile="../InputExamples/KinationInput.dat" # options for the solver # These variables are optional. Yoou can use the Axion class without them. initial_step_size=1e-2; #initial step the solver takes. minimum_step_size=1e-8; #This limits the sepsize to an upper limit. maximum_step_size=1e-2; #This limits the sepsize to a lower limit. absolute_tolerance=1e-8; #absolute tolerance of the RK solver relative_tolerance=1e-8; #relative tolerance of the RK solver beta=0.9; #controls how agreesive the adaptation is. Generally, it should be around but less than 1. #The stepsize does not increase more than fac_max, and less than fac_min. #This ensures a better stability. Ideally, fac_max=inf and fac_min=0, but in reality one must #tweak them in order to avoid instabilities. fac_max=1.2; fac_min=0.8; maximum_No_steps=int(1e7); #maximum steps the solver can take Quits if this number is reached even if integration is not finished. _=time() # AxionMass instance axionMass = AxionMass(r'../../src/data/chi.dat',0,mP) #------------------------------------------------------------------------------# # this is the axion mass squared beyond the interpolation limits for the current data # if yo don't specify it, the axion mass is taken to be constant beyond these limits TMax=axionMass.getTMax() chiMax=axionMass.getChiMax() TMin=axionMass.getTMin() chiMin=axionMass.getChiMin() axionMass.set_ma2_MAX( lambda T,fa: chiMax/fa/fa*pow(T/TMax,-8.16) ) axionMass.set_ma2_MIN( lambda T,fa: chiMin/fa/fa ) #------------------------------------------------------------------------------# # def ma2(T,fa): # TQCD=150*1e-3 # ma20=3.1575e-05/fa/fa # if T<=TQCD: # return ma20 # return ma20*pow((TQCD/T),8.16) # axionMass = AxionMass(ma2) # Axion instance ax=Axion(theta_i, fa, umax, TSTOP, ratio_ini, N_convergence_max, convergence_lim, inputFile, axionMass, initial_step_size,minimum_step_size, maximum_step_size, absolute_tolerance, relative_tolerance, beta, fac_max, fac_min, maximum_No_steps) # solve the EOM (this only gives you the relic, T_osc, theta_osc, and a_osc) ax.solveAxion() print(ax.theta_i, ax.fa, ax.theta_osc, ax.T_osc ,ax.relic) print(round(time()-_,3),file=stderr) # change to True in order to mkae the plots if False: import matplotlib.pyplot as plt from numpy import array as np_array from numpy import abs as np_abs ax.getPeaks()#this gives you the peaks of the oscillation ax.getPoints()#this gives you all the points of integration ax.getErrors()#this gives you local errors of integration #########-----\theta-----######### fig=plt.figure(figsize=(9,4)) fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.9, right=0.9,wspace=0.0,hspace=0.25) sub = fig.add_subplot(1,1,1) #this plot shows the peaks of the oscillation sub.plot(ax.T_peak,ax.theta_peak,linestyle=':',marker='+',color='xkcd:blue',linewidth=2) #this plot shows all the points sub.plot(ax.T,ax.theta,linestyle='-',linewidth=2,alpha=1,c='xkcd:black') sub.set_yscale('linear') sub.set_xscale('linear') sub.set_xlabel(r'$T ~[{\rm GeV}]$') sub.xaxis.set_label_coords(0.5, -0.1) sub.set_ylabel(r'$\theta$') sub.yaxis.set_label_coords(-0.1,0.5) sub.axhline(ax.theta_osc,linestyle=':',color='xkcd:red',linewidth=1.5) sub.axvline(ax.T_osc,linestyle='--',color='xkcd:gray',linewidth=1.5) fig.savefig('theta-T_examplePlot.pdf',bbox_inches='tight') #########-----\zeta-----######### fig=plt.figure(figsize=(9,4)) fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.9, right=0.9,wspace=0.0,hspace=0.25) sub = fig.add_subplot(1,1,1) sub.plot(ax.T,ax.zeta,linestyle='-',linewidth=2,alpha=1,c='xkcd:black') sub.plot(ax.T_peak,ax.zeta_peak,linestyle=':',marker='+',color='xkcd:blue',linewidth=2) sub.set_yscale('linear') sub.set_xscale('linear') sub.set_xlabel(r'$T ~[{\rm GeV}]$') sub.xaxis.set_label_coords(0.5, -0.1) sub.set_ylabel(r'$\zeta $') sub.yaxis.set_label_coords(-0.1,0.5) sub.axvline(ax.T_osc,linestyle='--',color='xkcd:gray',linewidth=1.5) fig.savefig('zeta-T_examplePlot.pdf',bbox_inches='tight') #########-----errors-----######### fig=plt.figure(figsize=(9,4)) fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.9, right=0.9,wspace=0.0,hspace=0.25) sub = fig.add_subplot(1,1,1) sub.plot(ax.T,np_abs(ax.dtheta/ax.theta),linestyle='-',linewidth=2,alpha=1,c='xkcd:black',label=r'$\dfrac{\delta \theta}{\theta}$') sub.plot(ax.T,np_abs(ax.dzeta/ax.zeta),linestyle='-',linewidth=2,alpha=1,c='xkcd:red',label=r'$\dfrac{\delta \zeta}{\zeta}$') sub.set_yscale('log') sub.set_xscale('linear') sub.set_xlabel(r'$T ~[{\rm GeV}]$') sub.xaxis.set_label_coords(0.5, -0.1) sub.set_ylabel(r'local errors') sub.yaxis.set_label_coords(-0.1,0.5) sub.legend(bbox_to_anchor=(1, 0.0),borderaxespad=0., borderpad=0.05,ncol=1,loc='lower right',fontsize=14,framealpha=0) sub.axvline(ax.T_osc,linestyle='--',color='xkcd:gray',linewidth=1.5) fig.savefig('errors_examplePlot.pdf',bbox_inches='tight') #########-----rho_a-----######### fig=plt.figure(figsize=(9,4)) fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.9, right=0.9,wspace=0.0,hspace=0.25) sub = fig.add_subplot(1,1,1) sub.plot(ax.T,ax.rho_axion,linestyle='-',linewidth=2,alpha=1,c='xkcd:black') sub.plot(ax.T_peak,ax.rho_axion_peak,linestyle=':',linewidth=2,alpha=1,c='xkcd:blue') sub.set_yscale('log') sub.set_xscale('linear') sub.set_xlabel(r'$T ~[{\rm GeV}]$') sub.xaxis.set_label_coords(0.5, -0.1) sub.set_ylabel(r'$\rho_{a}(T) ~[{\rm GeV}^{4}]$') sub.yaxis.set_label_coords(-0.1,0.5) sub.axvline(ax.T_osc,linestyle='--',color='xkcd:gray',linewidth=1.5) fig.savefig('rho_a-T_examplePlot.pdf',bbox_inches='tight') #don't forget the destructor del ax \ No newline at end of file diff --git a/UserSpace/Python/AxionMass.py b/UserSpace/Python/AxionMass.py index 140efdf..14f6bb9 100755 --- a/UserSpace/Python/AxionMass.py +++ b/UserSpace/Python/AxionMass.py @@ -1,39 +1,40 @@ from numpy import logspace, pi from sys import path as sysPath -from os import path as osPath -sysPath.append(osPath.join(osPath.dirname(__file__), '../../src')) +from sys import path as sysPath +sysPath.append('../../src') + #load the module from interfacePy.AxionMass import AxionMass axionMass = AxionMass(r'../../src/data/chi.dat',0,1e5) # this is the axion mass squared beyond the interpolation limits for the current data TMax=axionMass.getTMax() chiMax=axionMass.getChiMax() TMin=axionMass.getTMin() chiMin=axionMass.getChiMin() axionMass.set_ma2_MAX( lambda T,fa: chiMax/fa/fa*pow(T/TMax,-8.16) ) axionMass.set_ma2_MIN( lambda T,fa: chiMin/fa/fa ) def ma2(T,fa): TQCD=150*1e-3; ma20=3.1575e-05/fa/fa; if T<=TQCD: return ma20; return ma20*pow((TQCD/T),8.16) axionMassFunc = AxionMass(ma2) print(axionMassFunc.ma2(3,1)) print(axionMass.ma2(3,1)) print(axionMass.ma2(0,1)) print(axionMassFunc.ma2(0,1)) del axionMass del axionMassFunc \ No newline at end of file diff --git a/UserSpace/Python/Cosmo.py b/UserSpace/Python/Cosmo.py index 2b87730..67217c1 100755 --- a/UserSpace/Python/Cosmo.py +++ b/UserSpace/Python/Cosmo.py @@ -1,94 +1,93 @@ from numpy import logspace from sys import path as sysPath -from os import path as osPath -sysPath.append(osPath.join(osPath.dirname(__file__), '../../src')) +sysPath.append('../../src') #load the module from interfacePy import Cosmo cosmo=Cosmo('../../src/data/eos2020.dat',0,1e5) for T in logspace(-5,5,50): print( 'T=',T,'GeV\t', 'H=',cosmo.Hubble(T),'GeV\t', 'h_eff=',cosmo.heff(T),'\t', 'g_eff=',cosmo.geff(T),'\t', 's=',cosmo.s(T),'GeV^3\t', ) if False: import matplotlib.pyplot as plt #########-----g_eff and h_eff-----######### fig=plt.figure(figsize=(9,4)) fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.95, right=0.9,wspace=0.0,hspace=0.0) fig.suptitle('') sub = fig.add_subplot(1,1,1) T=logspace(-5,5,500) gt=[cosmo.geff(i) for i in T] ht=[cosmo.heff(i) for i in T] sub.plot(T,gt,linestyle='--',c='xkcd:red',label=r"$g_{\rm eff} (T)$") sub.plot(T,ht,linestyle=':',c='xkcd:black',label=r"$h_{\rm eff} (T)$") sub.set_xlabel(r'$T ~ [{\rm GeV}]$') sub.set_ylabel(r'rel. dof') sub.legend(bbox_to_anchor=(1, 0.0),borderaxespad=0., borderpad=0.05,ncol=1,loc='lower right',fontsize=14,framealpha=0) sub.set_yscale('log') sub.set_xscale('log') fig.savefig('rdofs-T_examplePlot.pdf',bbox_inches='tight') #########-----dg_effdT and dh_effdT-----######### fig=plt.figure(figsize=(9,4)) fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.95, right=0.9,wspace=0.0,hspace=0.0) fig.suptitle('') sub = fig.add_subplot(1,1,1) T=logspace(-5,5,500) dg=[cosmo.dgeffdT (i) for i in T] dh=[cosmo.dheffdT(i) for i in T] sub.plot(T,dg,linestyle='--',c='xkcd:red',label=r"$\dfrac{d g_{\rm eff}}{dT} (T)$") sub.plot(T,dh,linestyle=':',c='xkcd:black',label=r"$\dfrac{d h_{\rm eff}}{dT} (T)$") sub.set_xlabel(r'$T ~ [{\rm GeV}]$') sub.legend(bbox_to_anchor=(1, 0.5),borderaxespad=0., borderpad=0.05,ncol=1,loc='lower right',fontsize=14,framealpha=0) sub.set_yscale('symlog') sub.set_xscale('log') fig.savefig('drdofsdT-T_examplePlot.pdf',bbox_inches='tight') #########-----dh-----######### fig=plt.figure(figsize=(9,4)) fig.subplots_adjust(bottom=0.15, left=0.15, top = 0.95, right=0.9,wspace=0.0,hspace=0.0) fig.suptitle('') sub = fig.add_subplot(1,1,1) T=logspace(-5,5,500) dht=[cosmo.dh(i) for i in T] sub.plot(T,dht,linestyle='-',c='xkcd:black') sub.set_xlabel(r'$T ~ [{\rm GeV}]$') sub.set_ylabel(r'$\delta_h = 1 + \dfrac{1}{3} \dfrac{d \log h_{\rm eff} }{d \log T}$') sub.set_yscale('linear') sub.set_xscale('log') fig.savefig('dh-T_examplePlot.pdf',bbox_inches='tight')