#include < iomanip >
        #include"MiMeS.hpp"
        
        using numeric = long double;//make life easier if you want to change to double
        
        int main(){
            mimes::util::Timer _timer_;//use this to time it!
        
            // use chi_PATH to interpolate the axion mass.
            mimes::AxionMass < numeric > axionMass(chi_PATH,0,mimes::Cosmo < numeric >::mP);
        
            /*set $\maT^2$ for $T\geq T_{\rm max}$*/
            numeric TMax=axionMass.getTMax(), chiMax=axionMass.getChiMax();    
        
            axionMass.set_ma2_MAX(
                [&chiMax;,&TMax;](numeric T, numeric fa){ return chiMax/fa/fa*std::pow(T/TMax,-8.16);}
            );  
        
            /*set $\maT^2$ for $T\leq T_{\rm min}$*/
            numeric TMin=axionMass.getTMin(), chiMin=axionMass.getChiMin();    
        
            axionMass.set_ma2_MIN( 
                [&chiMin;,&TMin;](numeric T, numeric fa){ return chiMin/fa/fa;}
            );		
        
            /*this path contains the cosmology*/
            std::string inputFile = std::string(rootDir)+
                std::string("/UserSpace/InputExamples/MatterInput.dat");
        
            /*declare an instance of Axion*/
            mimes::Axion < numeric, 1, RODASPR2 < numeric > > ax(0.1, 1e16, 500, 1e-4, 1e3, 10, 1e-2, 
                            inputFile, &axionMass;, 1e-2, 1e-8, 1e-2, 1e-10, 1e-10, 0.85, 1.5, 0.85,
                            int(1e7) );
            /*solve the EOM!*/
            ax.solveAxion();
        
            std::cout<&ltstd;::setprecision(5)
            <<"theta_i="<&ltax.theta;_i<&ltstd;::setw(25)<<"f_a="<&ltax.fa;<<" GeV\n"<<"theta_osc~="<&ltax.theta;_osc
            <&ltstd;::setw(20)<<"T_osc~="<&ltax.T;_osc<<"GeV \n"<<"Omega h^2="<&ltax.relic;<<"\n";
            
            return 0;
        }
    

This script is simply run as python example.py. The output should look like:

 
                theta_i=0.1                     f_a=1e+16 GeV
                theta_osc~=0.045618             T_osc~=0.059194GeV 
                Omega h^2=0.1224
                0.0389824
            

        #you need these in order to print the time it took to run main in stderr
        from time import time; 
        from sys import stderr 
        
        #add the relative path for MiMeS/src (we assume here that this script is one birectory below the root of MiMeS)
        from sys import path as sysPath; sysPath.append('../src')   
        
        from interfacePy.AxionMass import AxionMass #import the AxionMass class
        from interfacePy.Axion import Axion #import the Axion class
        from interfacePy.Cosmo import mP #import the Planck mass from Cosmo
        
        def main():
        
            # AxionMass instance
            axionMass = AxionMass(r'../src/data/chi.dat',0,mP)
        
            # define $\maT^2$ for $T\leq T_{\rm min}$
            TMin, chiMin=axionMass.getTMin(), axionMass.getChiMin()
            
            axionMass.set_ma2_MIN( lambda T,fa: chiMin/fa/fa )
        
            # define $\maT^2$ for $T\geq T_{\rm max}$
            TMax, chiMax=axionMass.getTMax(), axionMass.getChiMax()
        
            axionMass.set_ma2_MAX( lambda T,fa: chiMax/fa/fa*pow(TMax/T,8.16))
        
            #in python it is more convenient to use relative paths
            inputFile="../UserSpace/InputExamples/MatterInput.dat"  
        
            ax = Axion(0.1, 1e16, 500, 1e-4, 1e3, 10, 1e-2, inputFile, axionMass, 
                1e-2, 1e-8, 1e-2, 1e-10, 1e-10, 0.85, 1.5, 0.85, int(1e7))
                
            ax.solveAxion()
        
            print(" theta_i=",ax.theta_i,"\t\t\t\t","f_a=",ax.fa," GeV\n","theta_osc~=",
                    ax.theta_osc,"\t","T_osc~=",ax.T_osc,"GeV \n","Omega h^2=",ax.relic)
        
            #once we are done we should run the destructor
            del ax,axionMass
        
        if __name__ == '__main__':
            _=time()
            main()
            print(round(time()-_,3),file=stderr)
        

This script is simply run as python example.py. The output should look like:

 
                        theta_i= 0.1                            f_a= 1e+16  GeV
                        theta_osc~= 0.04561752103101018         T_osc~= 0.05919445773033808 GeV 
                        Omega h^2= 0.12239928594726415
                    0.059