# Tidal-Love-Number-Project
README FOR ADMIXED COMPACT OBJECT STABILITY CURVE SOLVER~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

This document is intended to explain how to use this code. All inputs are in natural units (MeV units of energy for all quantities). This code produces stability curves for compact object
containing 2 admixed perfect fluids (to be denoted fluid 1 and fluid 2). An example of the uses of this is to investigate neutron stars admixed with dark matter.

Sections:
0: Define an EOS
1: Definition of radius
2: Make a stability curve
3: Using find_density_for_mass_ratio
4: Computing properties of a single NS

0: Define an EOS========================================================================================================================================================================

EOS are functions which take in a density and output a pressure, these can be custom definined as a typical function in python, or there are a couple premade options.

.The class FermiEOS can be used to create an equation of state of a relativistic Fermi gas. A FermiEOS object is initialised with arguments (in order) particle mass, 
mediator mass, and the gauge coupling. Then the object can be called as any other EOS.

.The class ArrayEOS can be used to create and equation of state from a set of densities and their corresponding pressures. An ArrayEOS object is initialised with one argument, a
numpy array: numpy.array([densities,pressures]) where densities is an ordered list of the densities and pressures is a list of the corresponding pressures. Its very important that 
the list of densities is ordered as otherwise the function will not work, this is also the case for other arrays listed as "ordered" in this document.

1: Definition of radius=================================================================================================================================================================

Radius definitions are functions that take in inputs, in order, numpy.array([radii,density1]),numpy.array([radii,density2]) and numpy.array([radii,mass_enc]) where radii are an 
ordered list of radii, density1 and density2 are the corresponding fluid 1 and 2 densities respectively, and mass_enc is the corresponding total enclosed mass. These
functions return a tuple (fluid 1 total radius, fluid 2 total radius, radius). These can be custom defined or there is one premade option.

.The class simple_rad_def creates a radius definition in which a fixed proportion of the entire mass of the compact object is within the radius. The object can be initialised with one
argument, this fixed proportion. 

2: Make a stability curve===============================================================================================================================================================

To make a stability curve, you can use the function mass_radius_lovenum.
This function outputs a tuple (Rs, Rs1, Rs2, Ms, Ms1, Ms2) each of which is an array of corresponding radius, entire fluid 1 radius, entire fluid 2 radius, mass, entire 
fluid 1 mass, entire fluid 2 mass. If a keyword argument lovenum = True then there will be an extra element in the tuple, Lambdas, which is the tidal deformability.
Unlike the input quantities, these quantities are in units km and solar masses.

This function with its minimum arguments is def mass_radius_lovenum(EOS1, EOS2, rhorange):
.EOS1 is the equation of state of fluid 1.
.EOS2 is the equation of state of fluid 2.
.rhorange is set of (total) central densities at which points will be plotted. This can be sent as e.g. an numpy array or a standard array.

There are also keyword arguments (kwargs):
.ratio - desired ratio between the masses of fluid 2 and fluid 1.
.rspan - range of r values solutions to be solved over. The upper bound should to be suffieiently large that the terminating condition is reached. A 2 element array [r_min,r_max] 
where r_min is the starting  radius for the solution and r_max is the largest possible end point. There is a terminating condition of stopping at 1e-8 times the central density which
makes it worth choosing a large r_max (e.g. r_max = 1e30) and a small r_min e.g. r_min = 1e5.
.lovenum - if true outputs tidal deformabilities as a final element to the output array which is an array of the corresponding tidal deformabilities.
.tol - this is "tolerance", the proportion of the ratio that is allowed as an error. Typically 1e-3 is good enough to obtain smooth looking stability curves.
.RadDef - the definition of the radius used.
.METHOD - the method used to solve the differential equation, scipy.solve_ivp is used so see the possible arguments from the documentation of that function.
.RunTol - the tolerance on the solving of the differential equation, its important to choose a sufficiently small tolerance else the lines will not be smooth.

3: Using find_density_for_mass_ratio====================================================================================================================================================

This function allows you to choose a ratio between the fluid 1 and 2 masses and a central density then outputs properties of the compact object sastifying the input properties.
This function has outputs a dictionary containing elements with the following keys:

"rhocentralDM" - the dark matter central density
"DM","MM","M" - the total dark matter mass, total nuclear matter mass and mass respectively.
"RDM","RM","R" - the total dark matter radius, total nuclear matter radius and radius repectively.
"ratio" - the total dark matter mass/total nuclear matter mass.

If a keyword argument lovenum = True then there will be an extra element in the dictionary with key "k2", which is the tidal love number.

This function with its minimum arguments is find_density_for_mass_ratio(EOS,EOSDM, ratio, tol,centraldensityT,rho2guessmin,rho2guessmax,rspan,lovenum,RunTol):
.EOS1 is the equation of state of fluid 1.
.EOS2 is the equation of state of fluid 2.
.ratio - desired ratio between the masses of fluid 2 and fluid 1.
.tol - this is "tolerance", the proportion of the ratio that is allowed as an error. Typically 1e-3 is good enough to obtain smooth looking stability curves.
.centraldensityT - the desired total central density
.rho2guessmin - the lower bound estimate for the fluid 2 density, if in doubt choose 0
.rho2guessmax - the upper bound estimate for the fluid 2 density, if in doubt choose rhocentral
.rspan - range of r values solutions to be solved over. The upper bound should to be suffieiently large that the terminating condition is reached. A 2 element array [r_min,r_max]
.lovenum - if true outputs tidal deformabilities as a final element to the output array which is an array of the corresponding tidal deformabilities.
.RunTol - the tolerance on the solving of the differential equation, its important to choose a sufficiently small tolerance else the lines will not be smooth.

There are also keyword arguments (kwargs):
.cumulative - If cumulative = False nothing happens. You can alternatively set cumulative to an array of central densities and it will try to estimate the next value and use them as 
upper and lower bounds, this is what mass_radius_lovenum does with previous runs. It does this by making a 2nd order taylor series with the last 3 data points (so the array must be 
longer than 3 for it to have any effect).
.RadDef - the definition of the radius used.
.METHOD - the method used to solve the differential equation, scipy.solve_ivp is used so see the possible arguments from the documentation of that function.

4: Computing properties of a single compact object======================================================================================================================================

First create an object of the Compact_Object class.
This object takes in 1 argument:
.state_eq - the EOS of the fluid 1.
and 1 keyword argument (kwarg):
.EOS2in - the EOS of fluid 2.

This object will have some defaults you can modify:
self.RadDef = raddefault
self.method = 'Radau'
self.tol = 0.000001
method is the method used to solve the ODE, RadDef is the definition of the radius and tol is the tolerance on the solving of the ODE. These can be easily modified.

next use self.profiles_data(self,rhocentral,rspan,rhocentralDM = 0.0)
this has the arguments:
.centraldensity1 - the fluid 1 central density
.rspan - range of r values solutions to be solved over. The upper bound should to be suffieiently large that the terminating condition is reached. A 2 element array [r_min,r_max]
and keyword argument:
.centraldensity2 - the fluid 2 central density

This will solve the differential equation and save the data in the class object:
.self.density1,self.density2,self.densityT are of form np.array([rs,densities]) where rs is an array of the radii corresponding to the fluid 1, fluid 2 and total densities respectively.
.self.pressure1,self.pressure2,self.pressureT are of form np.array([rs,pressures]) where rs is an array of the radii corresponding to the fluid 1, fluid 2 and total pressures respectively.
.self.mass1,self.mass2,self.massT are of form np.array([rs,masses]) where rs is an array of the radii corresponding to the fluid 1, fluid 2 and total masses respectively.
.self.M1,self.M2,self.M are the entire fluid 1, entire fluid 2 and total masses respectively
.self.R1,self.R2,self.R are the entire fluid 1, entire fluid 2 and total radii respectively
.self.centraldensity1,self.centraldensity2 are the fluid 1 and fluid 2 central densities

To calculate the tidal love number first use self.solvefory(), this takes no arguments but solves the output of the ODE solver (scipy.solve_ivp) in variable self.ysolution.
Next use self.find_love_number(), this again takes no arguments.
This saves new quantity to the Compact_Object class object:
.self.k2, the tidal love number.

To convert this to a tidal deformability the function k2_to_deformability which takes 3 inputs:
.Rs an array of the radii of the compact objects
.Ms an array of the masses of the compact objects
.k2s an array of the tidal love numbers of the compact objects
it then outputs and array of the corresponding tidal deformabilities. The reason this takes in arrays rather than individual values is as it is used to calculate all tidal deformabilities
at once for compact_objects of interest in mass_radius_lovenum