.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/end_to_end/plot_gj1214b_kempton23.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_end_to_end_plot_gj1214b_kempton23.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_end_to_end_plot_gj1214b_kempton23.py:


Observe a phase curve of a spotted star.
========================================

This example demonstrates stellar contamination of a phase curve.

A phase curve with a long enough baseline can be contaminated by
stellar variability. We take the phase curve of GJ1214 b, analyzed by :cite:t:`kempton+23`
using JWST MIRI-LRS, as an example.

.. GENERATED FROM PYTHON SOURCE LINES 11-26

.. code-block:: Python

    from pathlib import Path
    import numpy as np
    import matplotlib.pyplot as plt
    from astropy import units as u
    from cartopy import crs as ccrs
    import libpypsg
    from libpypsg.globes import PyGCM

    from VSPEC import ObservationModel,PhaseAnalyzer
    from VSPEC import params
    from VSPEC.config import MSH

    SEED = 1214
    libpypsg.docker.set_url_and_run()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Saved settings to /home/runner/.libpypsg/settings.json
    Reloading settings...




.. GENERATED FROM PYTHON SOURCE LINES 27-29

Create the configurations
-------------------------

.. GENERATED FROM PYTHON SOURCE LINES 29-180

.. code-block:: Python


    # Instrument
    inst = params.InstrumentParameters.miri_lrs()

    # Observation
    observation = params.ObservationParameters(
        observation_time=41.0*u.hr,
        integration_time=15*u.min
    )

    # PSG
    psg_params = params.psgParameters(
        gcm_binning=200,
        phase_binning=1,
        nmax=0,
        lmax=0,
        use_continuum_stellar=True,
        continuum=['Rayleigh', 'Refraction', 'CIA_all'],
        use_molecular_signatures=True
    )

    # Star and Planet
    star_teff = 3250*u.K
    star_rad = 0.215*u.R_sun
    planet_rad = 2.742*u.R_earth
    orbit_rad = 0.01490*u.AU
    orbit_period = 1.58040433*u.day
    planet_rot_period = orbit_period
    star_rot_period = 120 * u.day
    planet_mass = 8.17*u.M_earth
    star_mass = 0.178*u.M_sun
    inclination = 88.7*u.deg

    start_time_before_eclipse = 2*u.hr
    angle_before_eclipse = (2*np.pi*u.rad * start_time_before_eclipse/orbit_period).to(u.deg)
    initial_phase = 0*u.deg - angle_before_eclipse

    planet_params = params.PlanetParameters(
        name='GJ1214b',
        radius=planet_rad,
        gravity=params.GravityParameters('kg',planet_mass),
        semimajor_axis=orbit_rad,
        orbit_period=orbit_period,
        rotation_period=planet_rot_period,
        eccentricity=0,
        obliquity=0*u.deg,
        obliquity_direction=0*u.deg,
        init_phase=initial_phase,
        init_substellar_lon=0*u.deg
    )

    system_params = params.SystemParameters(
        distance=14.6427*u.pc,
        inclination=inclination,
        phase_of_periastron=0*u.deg
    )

    star_dict = {
        'teff': star_teff,
        'radius': star_rad
    }
    planet_dict = {'semimajor_axis': orbit_rad}

    gcm_dict = {
        'nlayer': 30,
        'nlon': 30,
        'nlat': 15,
        'epsilon': 6,
        'albedo': 0.3,
        'emissivity': 1.0,
        'lat_redistribution': 0.1,
        'gamma': 1.4,
        'psurf': 1*u.bar,
        'ptop': 1e-5*u.bar,
        'wind': {'U': '0 m/s','V':'0 m/s'},
        'molecules':{'CO2':0.99}
    }

    gcm = params.gcmParameters.from_dict({
        'star':star_dict,
        'planet':planet_dict,
        'gcm':{'vspec':gcm_dict,'mean_molec_weight':44}
    })

    star_kwargs = dict(
        psg_star_template='M',
        teff=star_teff,
        mass=star_mass,
        radius=star_rad,
        period=star_rot_period,
        misalignment=0*u.deg,
        misalignment_dir=0*u.deg,
        ld=params.LimbDarkeningParameters.proxima(),
        faculae=params.FaculaParameters.none(),
        flares=params.FlareParameters.none(),
        granulation=params.GranulationParameters.none(),
        grid_params=(500,1000),
    )

    quiet_star = params.StarParameters(
        spots=params.SpotParameters.none(),
        **star_kwargs
    )
    spotted_star = params.StarParameters(
        spots=params.SpotParameters(
            distribution='iso',
            initial_coverage=0.2,
            area_mean=300*MSH,
            area_logsigma=0.2,
            teff_umbra=2700*u.K,
            teff_penumbra=2700*u.K,
            equillibrium_coverage=0.2,
            burn_in=0*u.s,
            growth_rate=0.0/u.day,
            decay_rate=0*MSH/u.day,
            initial_area=10*MSH
        ),
        **star_kwargs
    )

    # Set parameters for simulation
    header_kwargs = dict(
        seed = SEED,
        verbose = 0,
        spec_grid = params.VSPECGridParameters(
            max_teff=3400*u.K,min_teff=2300*u.K,
            impl_bin='rust',impl_interp='scipy',fail_on_missing=False
        )
    )
    internal_params_kwargs = dict(
        planet=planet_params,
        system=system_params,
        obs=observation,
        gcm=gcm,
        psg=psg_params,
        inst=inst
    )

    params_quiet = params.InternalParameters(
        header=params.Header(data_path=Path('.vspec/gj1214_quiet'),**header_kwargs),
        star = quiet_star,
        **internal_params_kwargs
    )

    params_spotted = params.InternalParameters(
        header=params.Header(data_path=Path('.vspec/gj1214_spotted'),**header_kwargs),
        star = spotted_star,
        **internal_params_kwargs
    )









.. GENERATED FROM PYTHON SOURCE LINES 181-186

Map the planetary surface
-------------------------

Before we run ``VSPEC``, let's look at the planet.


.. GENERATED FROM PYTHON SOURCE LINES 186-208

.. code-block:: Python


    gcm_data:PyGCM = gcm.get_gcm()

    tsurf = gcm_data.tsurf.dat.to_value(u.K)

    fig = plt.figure()
    proj = ccrs.Robinson(central_longitude=0)
    ax = fig.add_subplot(projection=proj)

    lats = gcm_data.lat_start + np.arange(gcm_data.shape[2]+1) * gcm_data.dlat.to_value(u.deg)
    lons = gcm_data.lon_start + np.arange(gcm_data.shape[1]+1) * gcm_data.dlon.to_value(u.deg)

    im = ax.pcolor(lons,lats,tsurf.T,cmap='gist_heat',transform=ccrs.PlateCarree())
    gl = ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,
        color='grey', alpha=0.8, linestyle='--')
    gl.top_xlabels = False
    gl.right_ylabels = False


    _=fig.colorbar(im,ax=ax,label='Surface Temperature (K)')





.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_001.png
   :alt: plot gj1214b kempton23
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 209-212

Run the spotless model
----------------------


.. GENERATED FROM PYTHON SOURCE LINES 212-217

.. code-block:: Python


    model_quiet = ObservationModel(params=params_quiet)
    model_quiet.build_planet()
    model_quiet.build_spectra()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Creating interpolators:
    thermal
    thermal, combined
    thermal, combined, stellar
    thermal, combined, stellar, photon noise
    thermal, combined, stellar, photon noise, detector noise
    thermal, combined, stellar, photon noise, detector noise, telescope noise
    thermal, combined, stellar, photon noise, detector noise, telescope noise, background noise
    thermal, combined, stellar, photon noise, detector noise, telescope noise, background noise, transit
    Finished!

    Loading Spectra:   0%|          | 0/12 [00:00<?, ?it/s]
    Loading Spectra:   8%|▊         | 1/12 [00:00<00:03,  3.35it/s]
    Loading Spectra:  17%|█▋        | 2/12 [00:00<00:02,  3.35it/s]
    Loading Spectra:  25%|██▌       | 3/12 [00:00<00:02,  3.35it/s]
    Loading Spectra:  33%|███▎      | 4/12 [00:01<00:02,  3.35it/s]
    Loading Spectra:  42%|████▏     | 5/12 [00:01<00:02,  3.35it/s]
    Loading Spectra:  50%|█████     | 6/12 [00:01<00:01,  3.36it/s]
    Loading Spectra:  58%|█████▊    | 7/12 [00:02<00:01,  3.36it/s]
    Loading Spectra:  67%|██████▋   | 8/12 [00:02<00:01,  3.37it/s]
    Loading Spectra:  75%|███████▌  | 9/12 [00:02<00:00,  3.36it/s]
    Loading Spectra:  83%|████████▎ | 10/12 [00:02<00:00,  3.36it/s]
    Loading Spectra:  92%|█████████▏| 11/12 [00:03<00:00,  3.36it/s]
    Loading Spectra: 100%|██████████| 12/12 [00:03<00:00,  3.36it/s]
    Loading Spectra: 100%|██████████| 12/12 [00:03<00:00,  3.36it/s]




.. GENERATED FROM PYTHON SOURCE LINES 218-220

Plot the lightcurve
-------------------

.. GENERATED FROM PYTHON SOURCE LINES 220-272

.. code-block:: Python


    data_quiet = PhaseAnalyzer(model_quiet.directories['all_model'])
    flux_unit = u.Unit('W m-2 um-1')
    def get_star(data:PhaseAnalyzer):
        i_eclipse1 = np.argmin(data.lightcurve('total',(0,-1))[:data.n_images//4])
        i_eclipse2 = np.argmin(data.lightcurve('total',(0,-1))[3*data.n_images//4:]) + 3*data.n_images//4
        time = (data.time-data.time[0]).to_value(u.hr)
        star_spec1 = data.spectrum('total',i_eclipse1).to_value(flux_unit)
        star_spec2 = data.spectrum('total',i_eclipse2).to_value(flux_unit)
    
        def func(t:float):
            m = (star_spec2 - star_spec1)/(time[i_eclipse2]-time[i_eclipse1])
            x = t-time[i_eclipse1]
            b = star_spec1
            y = m * x + b
            return y

        return func


    def plot_lc(data:PhaseAnalyzer):
        fig,axes = plt.subplots(2,1,tight_layout=True)

        axes[0].scatter((data.time-data.time[0]).to(u.hr),
            data.lightcurve('total',(0,-1)),label='white light',s=5,c='k')
        axes[0].set_xlabel('Time since start of observation (hour)')
        axes[0].set_ylabel('Flux (W m-2 um-1)')
        axes[0].legend()
        first_four = data.time-data.time[0] <= 4*u.hour
        axins = axes[0].inset_axes([0.08, 0.15, 0.35, 0.5])
        axins.scatter((data.time-data.time[0]).to(u.hr)[first_four],
            data.lightcurve('total',(0,-1))[first_four],label='white light',s=5,c='k')
        axes[0].indicate_inset_zoom(axins)

        interp = get_star(data)
        t = (data.time-data.time[0]).to_value(u.hr)

        n_steps = 10
        colors = plt.get_cmap('viridis')
        indices = np.arange(start=0,stop=data.n_images,step=data.n_images//n_steps)

        for index in indices:
            star_spec = interp(t[index])
            pl_spec = data.spectrum('total',index).to_value(flux_unit) - star_spec
            axes[1].plot(data.wavelength,1e6*pl_spec/star_spec,c=colors(index/data.n_images))
        axes[1].set_xlabel('Wavelength (um)')
        axes[1].set_ylabel('Planet flux (ppm)')
        return fig

    plot_lc(data_quiet).show()





.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_002.png
   :alt: plot gj1214b kempton23
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 273-277

Plot the spectroscopic phase curve
----------------------------------

We can throw out the transit points while we're at it.

.. GENERATED FROM PYTHON SOURCE LINES 277-314

.. code-block:: Python



    def get_phase_map(data:PhaseAnalyzer):
        white_light_curve = data_quiet.lightcurve('total',(0,-1),normalize=0)
        points_to_use = white_light_curve > 0.5*(np.median(white_light_curve)+ np.min(white_light_curve))

        interp = get_star(data)
        ts = (data.time-data.time[0]).to_value(u.hr)

        # get the planet flux, except plance nan during transit
        star_im = np.array([interp(t) for t in ts]).T
        total_im = data.total.to_value(flux_unit)
        pl_im = np.where(
            points_to_use,
            total_im-star_im,
            np.nan
        )
        return pl_im,star_im
    

    def plot_phasecurve(data:PhaseAnalyzer):
        pl_im,star_im = get_phase_map(data)

        fig,ax = plt.subplots(1,1)
        im = ax.pcolormesh(
            (data.time-data.time[0]).to_value(u.hr),
            data.wavelength.to_value(u.um),
            pl_im/star_im*1e6,
            cmap='viridis'
        )
        fig.colorbar(im,ax=ax,label='Planet flux (ppm)')
        ax.set_xlabel('Time since start of observation (hour)')
        ax.set_ylabel('Wavelength (um)')
        return fig

    plot_phasecurve(data_quiet).show()




.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_003.png
   :alt: plot gj1214b kempton23
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 315-324

We can easily make out the phase curve because the star is static.

Run the spotted model
---------------------

Because we are using the same planet parameters, we won't rerun PSG
for this. Instead, we will just rerun the stellar part of the code.
In a way this is cheating but it will save time. Be careful because we
are overwriting our old data.

.. GENERATED FROM PYTHON SOURCE LINES 324-331

.. code-block:: Python


    model_spotted = ObservationModel(params_spotted)
    model_spotted.build_planet()
    model_spotted.build_spectra()

    data_spotted = PhaseAnalyzer(model_spotted.directories['all_model'])





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Generated 124 mature spots
    Creating interpolators:
    thermal
    thermal, combined
    thermal, combined, stellar
    thermal, combined, stellar, photon noise
    thermal, combined, stellar, photon noise, detector noise
    thermal, combined, stellar, photon noise, detector noise, telescope noise
    thermal, combined, stellar, photon noise, detector noise, telescope noise, background noise
    thermal, combined, stellar, photon noise, detector noise, telescope noise, background noise, transit
    Finished!

    Loading Spectra:   0%|          | 0/12 [00:00<?, ?it/s]
    Loading Spectra:   8%|▊         | 1/12 [00:00<00:03,  3.27it/s]
    Loading Spectra:  17%|█▋        | 2/12 [00:00<00:03,  3.28it/s]
    Loading Spectra:  25%|██▌       | 3/12 [00:00<00:02,  3.29it/s]
    Loading Spectra:  33%|███▎      | 4/12 [00:01<00:02,  3.30it/s]
    Loading Spectra:  42%|████▏     | 5/12 [00:01<00:02,  3.30it/s]
    Loading Spectra:  50%|█████     | 6/12 [00:01<00:01,  3.31it/s]
    Loading Spectra:  58%|█████▊    | 7/12 [00:02<00:01,  3.29it/s]
    Loading Spectra:  67%|██████▋   | 8/12 [00:02<00:01,  3.30it/s]
    Loading Spectra:  75%|███████▌  | 9/12 [00:02<00:00,  3.30it/s]
    Loading Spectra:  83%|████████▎ | 10/12 [00:03<00:00,  3.31it/s]
    Loading Spectra:  92%|█████████▏| 11/12 [00:03<00:00,  3.31it/s]
    Loading Spectra: 100%|██████████| 12/12 [00:03<00:00,  3.31it/s]
    Loading Spectra: 100%|██████████| 12/12 [00:03<00:00,  3.30it/s]




.. GENERATED FROM PYTHON SOURCE LINES 332-336

Plot the lightcurve, again
--------------------------

We redo our earlier analysis

.. GENERATED FROM PYTHON SOURCE LINES 336-339

.. code-block:: Python


    plot_lc(data_spotted).show()




.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_004.png
   :alt: plot gj1214b kempton23
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 340-343

And the phase curve
-------------------


.. GENERATED FROM PYTHON SOURCE LINES 343-346

.. code-block:: Python


    plot_phasecurve(data_spotted).show()




.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_005.png
   :alt: plot gj1214b kempton23
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 347-353

Compare phase curves
--------------------

:cite:t:`kempton+23` break the spectrum up into :math:`0.5 \mu m` bins
to analyze the phase curve. We will do the same to observe the effects of
stellar contamination.

.. GENERATED FROM PYTHON SOURCE LINES 353-388

.. code-block:: Python


    def get_lc(data:PhaseAnalyzer,w1:u.Quantity,w2:u.Quantity):
        """Get the lightcurve given a bandpass"""
        wl = data.wavelength
        i_left = int(np.argwhere(wl > w1)[0])
        try:
            i_right = int(np.argwhere(wl > w2)[0])
        except IndexError:
            i_right = -1
        interp = get_star(data)
        ts = (data.time-data.time[0]).to_value(u.hr)
        star_im = np.array([interp(t) for t in ts]).T
        total_im = data.total.to_value(flux_unit)
        pl_im = total_im-star_im
        lc = 1e6*pl_im[i_left:i_right,:]/star_im[i_left:i_right,:]
        return lc.mean(axis=0)

    bin_edges = np.arange(5.0,12.0,0.5)
    n_ax = len(bin_edges)
    fig,axes = plt.subplots(n_ax,1,figsize=(7,10),sharex=True)
    for edge,ax in zip(bin_edges,axes):
        w1,w2 = edge*u.um,(edge+0.5)*u.um
        quiet_lc = get_lc(data_quiet,w1,w2)
        spotted_lc = get_lc(data_spotted,w1,w2)
        time = (data_quiet.time - data_quiet.time[0]).to(u.hr)
        ax.plot(time,(quiet_lc),c='xkcd:azure',label='No Spots')
        ax.plot(time,(spotted_lc),c='xkcd:lavender',label='Spotted')
        ax.text(0.7,0.7,f'{w1:.1f} - {w2:.1f}',transform=ax.transAxes)
        ax.set_ylim(-100,700)
    fig.subplots_adjust(hspace=0,wspace=0)
    axes[0].legend()
    axes[-1].set_xlabel('Time (hour)')
    _ = axes[n_ax//2].set_ylabel('Planet Flux (ppm)')





.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_006.png
   :alt: plot gj1214b kempton23
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_006.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    /home/runner/work/VSPEC/VSPEC/examples/end_to_end/plot_gj1214b_kempton23.py:357: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
      i_left = int(np.argwhere(wl > w1)[0])
    /home/runner/work/VSPEC/VSPEC/examples/end_to_end/plot_gj1214b_kempton23.py:359: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
      i_right = int(np.argwhere(wl > w2)[0])




.. GENERATED FROM PYTHON SOURCE LINES 389-394

2D residuals
------------

Let's take a look at how much of the planet flux (from the spotted model)
is actaully due to spots.

.. GENERATED FROM PYTHON SOURCE LINES 394-408

.. code-block:: Python


    pl_quiet,star_quiet = get_phase_map(data_quiet)
    pl_spotted,_ = get_phase_map(data_spotted)

    contribution_from_spots = pl_spotted-pl_quiet
    contrast = (contribution_from_spots/star_quiet*1e6)
    t = (data_quiet.time - data_quiet.time[0]).to_value(u.hr)
    wl = data_quiet.wavelength.to_value(u.um)

    fig,ax = plt.subplots(1,1)
    im = ax.pcolormesh(t,wl,contrast)
    fig.colorbar(im,ax=ax,label='False planet flux (ppm)')
    ax.set_xlabel('Time since start of observation (hour)')
    _=ax.set_ylabel('Wavelength (um)')



.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_007.png
   :alt: plot gj1214b kempton23
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_gj1214b_kempton23_007.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (11 minutes 20.780 seconds)


.. _sphx_glr_download_auto_examples_end_to_end_plot_gj1214b_kempton23.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_gj1214b_kempton23.ipynb <plot_gj1214b_kempton23.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_gj1214b_kempton23.py <plot_gj1214b_kempton23.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_gj1214b_kempton23.zip <plot_gj1214b_kempton23.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_