Using custom functions

Custom functions allow the user to modify the model used in CONAN to fit lightcurves, RVs, or indeed any type of dataset.

There are two types of custom functions that can be implemented in CONAN:

Model-modifying function:

Here the custom function is designed to modify CONAN's built-in LC or RV model (see here). For example, including a Rossiter-McLaughlin (RM) anomaly signal to the RV model requires the user to define the custom_RM_model() and an operator function op_func() that defines how to combine the output of the RM model with the inbuilt RV model. An example is given below:

def custom_func(phase_angle, A, B, C,extra_args={}):     # A,B,C are the custom parameters and no extra arguments
    return A*np.sin(phase_angle) + B*np.cos(2*phase_angle) + C

def op_func(transit_model, custom_model):   # operation function to combine the custom model with the transit model
    return transit_model + custom_model

# A is fixed, B has uniform prior, C has gaussian prior
custom_lcfunc = lc_obj.add_custom_LC_function(func      = custom_func, 
                                              x         = "phase_angle",
                                              func_args = dict(A=0.1, B=(0.2,0,0.01), C=(0.3,0.01), 
                                              op_func   = op_func)

# this custom function has now been registered and will be used in the light-curve model.

Model-replacing function:

Here the custom function completely replaces CONAN's built-in LC or RV model. For example, to model transits of non-spherical planets, the user will define the model and specify to replace_model. This functionality is useful as it allows the user to fit any type of dataset with their own model. An example of using catwoman to fit non-spherical planet shape due to limb asymmetry is shown below:

def catwoman_func(t, rp2, phi, extra_args=dict(ld_law="quadratic"), 
                 LC_pars=dict(Duration=0,rho_star=0,RpRs=0,Impact_para=0,T_0=0,
                              Period=0,Eccentricity=0,omega=90,q1=0,q2=0,D_occ=0,
                              Fn=0,ph_off=0,A_ev=0,f1_ev=0,A_db=0) ):
    import catwoman
    import numpy as np
    import astropy.constants as c
    import astropy.units as u

    # create a catwoman model. CONAN transit pars in LDpars can be passed to the function
    params  = catwoman.TransitParams()
    params.t0  = LC_pars["T_0"]
    params.per = LC_pars["Period"]
    params.rp  = LC_pars["RpRs"]
    #convert stellar density to a/R*
    rho = LC_pars["rho_star"]
    G  = (c.G.to(u.cm**3/(u.g*u.second**2))).value
    aR = ( rho*G*(P*(u.day.to(u.second)))**2 / (3*np.pi)) **(1/3.)
    params.a   =  aR
    params.inc = np.arccos(LC_pars["Impact_para"]/aR)
    params.ecc = LC_pars["Eccentricity"]
    params.w   = LC_pars["omega"]
    params.limb_dark = extra_args["ld_law"]
    #convert from kipping parameterisation to quadratic
    u1 = 2*np.sqrt(LC_pars["q1"])*LC_pars["q2"]
    u2 = np.sqrt(LC_pars["q1"])*(1-2*LC_pars["q2"])
    params.u = [u1, u2]

    params.phi = phi                        #angle of rotation of top semi-circle (in degrees)
    params.rp2 = rp2                        #bottom semi-circle radius (in units of stellar radii)

    model = catwoman.TransitModel(params,t)         #initalises model
    return  model.light_curve(params)                #calculates light curve

lc_obj.add_custom_LC_function(func            = catwoman_func, 
                              x               = "time",
                              func_args       = dict(rp2=(0.1,0.01),phi=(-90,0,90)),  # the only 2 non-standard parameters
                              extra_args      = dict(ld_law="quadratic") 
                              op_func         = None,
                              replace_LCmodel = True)