CONAN.utils.CA_08_PCmodel
=========================

.. py:function:: CONAN.utils.CA_08_PCmodel(phi, Fd, C1, D1, C2, D2)

   Calculate the phase curve of a planet using the Cowan & Agol 2008 model.
   See also https://iopscience.iop.org/article/10.3847/1538-3881/ae0e52/meta
   and https://arxiv.org/pdf/2305.06240

   The equation is given as F = Fd + C1*(cos(theta) - 1) + D1*sin(theta) + C2*(cos(2*theta) - 1) + D2*sin(2*theta)
   where theta = phi + pi, and phi is the phase angle of the planet (true anomaly+omega-pi/2) in radians.

   :param phi: phase angle (2*pi*phase for circular orbit) or true anomaly+omega-pi/2 in radians.
   :type phi: array-like
   :param Fd: Dayside flux/occultation depth
   :type Fd: float
   :param C1: first cosine coefficient
   :type C1: float
   :param D1: first sine coefficient
   :type D1: float
   :param C2: second cosine coefficient
   :type C2: float
   :param D2: second sine coefficient
   :type D2: float

   :returns: **F** -- planetary flux as a function of phase angle
   :rtype: array-like

   .. rubric:: Example

   >>> from utils import CA_08_PCmodel
   >>> import numpy as np
   >>> phases = np.linspace(-0.25, 0.75, 10000)
   >>> phi = 2 * np.pi * phases
   >>> Fn, Fd = 50e-6, 400e-6
   >>> C1 = (Fd - Fn)/2
   >>> D1, C2, D2 = -139e-6, 46e-6, -15e-6  #D1 somehow influences the phase shift. find an expression for it
   >>> res = CA_08_PCmodel(phi, Fd, C1, D1, C2, D2)

   >>> from CONAN.utils import cosine_atm_variation
   >>> pc = CA_08_PCmodel(phi, Fd, C1, D1, C2, D2)
   >>> cos = cosine_atm_variation(phi, Fd, Fn,delta).pc
   >>> plt.figure()
   >>> plt.plot(phases, pc)
   >>> plt.plot(phases, cos)
   >>> plt.axvline(0.5, color='k', ls='--')
   >>> # plt.axvline(0.5 - np.radians(delta) / (2 * np.pi))
   >>> plt.axhline(Fd-2*C1, color='k', ls='--')
   >>> plt.axvline(0)