No OneTemporary
Actions

Size

12 KB

Referenced Files

None

Subscribers

None

View Options

	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')

File Metadata

Mime Type: text/x-diff
Expires: Tue, Sep 30, 6:13 AM (11 h, 59 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6531308
Default Alt Text: (12 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions