.. 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_spotted_transit.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_spotted_transit.py>`
        to download the full example code.

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

.. _sphx_glr_auto_examples_end_to_end_plot_spotted_transit.py:


Observe a transit of a spotted star.
========================================

This example demonstrates stellar contamination of a transit spectrum.

A transit can change the spot coverage of a star and produce a signal
that is difficult to distinguish from atmospheric absorption. We aim to simulate
data from :cite:t:`2023ApJ...948L..11M`.

.. GENERATED FROM PYTHON SOURCE LINES 11-28

.. code-block:: Python

    from pathlib import Path
    import numpy as np
    import matplotlib.pyplot as plt
    from astropy import units as u
    import libpypsg

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

    SEED = 10
    PRE_TRANSIT = 8
    IN_TRANSIT = 13
    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 29-34

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

:cite:t:`2023ApJ...948L..11M` observed super-Earth GJ 486b with
JWST NIRSpec/G395H

.. GENERATED FROM PYTHON SOURCE LINES 34-252

.. code-block:: Python


    # Instrument

    inst = params.InstrumentParameters(
        telescope=params.SingleDishParameters.jwst(),
        bandpass=params.BandpassParameters(
            wl_blue=2.87*u.um,
            wl_red=5.14*u.um,
            resolving_power=200,
            wavelength_unit=u.um,
            flux_unit=u.Unit('W m-2 um-1')
        ),
        detector=params.DetectorParameters(
            beam_width=5*u.arcsec,
            integration_time=0.5*u.s,
            ccd=params.ccdParameters(
                pixel_sampling=8,
                read_noise=16.8*u.electron,
                dark_current=0.005*u.electron/u.s,
                throughput=0.3,
                emissivity=0.1,
                temperature=50*u.K
            )
        )
    )

    # Observation

    observation = params.ObservationParameters(
        observation_time=3.53*u.hr,
        integration_time=8*u.min
    )

    # PSG

    psg_kwargs = dict(
        gcm_binning=200,
        phase_binning=1,
        nmax=0,
        lmax=0,
        use_continuum_stellar=True,
        continuum=['Rayleigh', 'Refraction'],
    )
    psg_params = params.psgParameters(
        use_molecular_signatures=True,
        **psg_kwargs
    )
    psg_no_atm = params.psgParameters(
        use_molecular_signatures=False,
        **psg_kwargs
    )

    # Star and Planet

    star_teff = 3343*u.K
    star_rad = 0.339*u.R_sun
    a_rstar = 11.229 # Moran+23 Table 1
    rp_rstar = 0.03709 # Moran+23 Table 1
    planet_rad = star_rad * rp_rstar
    orbit_rad = star_rad * a_rstar
    orbit_period = 1.467119*u.day
    planet_rot_period = orbit_period
    star_rot_period = np.inf*u.s # assume the star does not change.
    planet_mass = 2.82*u.M_earth
    star_mass = 0.323*u.M_sun
    inclination = 89.06*u.deg
    planet_norm_factor = {
        'rock_quiet': 1,
        'rock_spotted': 0.97,
        'water_quiet': 0.98,
    } # We will multiply the radius by these scalars rather than fit
      # later (because I have run this code before and know how much
      # each model is off by).

    print(f'Planet radius: {planet_rad.to_value(u.R_earth)} R_earth')
    print(f'Semimajor axis: {orbit_rad.to_value(u.AU)} AU')

    observation_angle = (2*np.pi*u.rad * observation.observation_time/orbit_period).to(u.deg)
    initial_phase = 180*u.deg - 0.5*observation_angle

    planet_kwargs = dict(
        name='GJ486b',
        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
    )

    quiet_rock_planet = params.PlanetParameters(
        radius=planet_rad*planet_norm_factor['rock_quiet'],
        **planet_kwargs
    )
    quiet_water_planet = params.PlanetParameters(
        radius=planet_rad*planet_norm_factor['water_quiet'],
        **planet_kwargs
    )
    spotted_rock_planet = params.PlanetParameters(
        radius=planet_rad*planet_norm_factor['rock_spotted'],
        **planet_kwargs
    )


    system_params = params.SystemParameters(
        distance=8.07*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': 40,
        'nlon': 60,
        'nlat': 30,
        'epsilon': 100,
        'albedo': 0.3,
        'emissivity': 1.0,
        'lat_redistribution': 0.0,
        'gamma': 1.0,
        'psurf': 1*u.bar,
        'ptop': 1e-10*u.bar,
        'wind': {'U': '0 m/s','V':'0 m/s'},
        'molecules':{'H2O':0.99}
    }

    # Create two sets of GCM Parameters

    h2o_atm = {'molecules':{'H2O':0.99}}
    gcm_h2o = params.gcmParameters.from_dict({
        'star':star_dict,
        'planet':planet_dict,
        'gcm':{'vspec':dict(gcm_dict,**h2o_atm),'mean_molec_weight':18}
    })
    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.lambertian(),
        faculae=params.FaculaParameters.none(),
        flares=params.FlareParameters.none(),
        granulation=params.GranulationParameters.none(),
        grid_params=100000,
    )
    quiet_star = params.StarParameters(
        spots=params.SpotParameters.none(),
        **star_kwargs
    )
    spotted_star = params.StarParameters(
        spots=params.SpotParameters(
            distribution='iso',
            initial_coverage=0.11,
            area_mean=500*MSH,
            area_logsigma=0.2,
            teff_umbra=2700*u.K,
            teff_penumbra=2700*u.K,
            equillibrium_coverage=0.0,
            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,
        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(
        system=system_params,
        obs=observation,
        gcm=gcm_h2o,
        inst=inst
    )

    # Make the three cases

    params_rock_quiet = params.InternalParameters(
        header=params.Header(data_path=Path('.vspec/rock_quiet'),**header_kwargs),
        star = quiet_star,
        psg=psg_no_atm,
        planet=quiet_rock_planet,
        **internal_params_kwargs
    )
    params_h2o_quiet = params.InternalParameters(
        header=params.Header(data_path=Path('.vspec/h2o_quiet'),**header_kwargs),
        star = quiet_star,
        psg=psg_params,
        planet=quiet_water_planet,
        **internal_params_kwargs
    )

    params_rock_spotted = params.InternalParameters(
        header=params.Header(data_path=Path('.vspec/rock_spotted'),**header_kwargs),
        star = spotted_star,
        psg=psg_no_atm,
        planet=spotted_rock_planet,
        **internal_params_kwargs
    )





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

 .. code-block:: none

    Planet radius: 1.371472837835719 R_earth
    Semimajor axis: 0.017702612840063636 AU




.. GENERATED FROM PYTHON SOURCE LINES 253-257

Run VSPEC for the simplest case
-------------------------------

We read in the config file and run the model.

.. GENERATED FROM PYTHON SOURCE LINES 257-262

.. code-block:: Python


    model_rock_quiet = ObservationModel(params_rock_quiet)
    model_rock_quiet.build_planet()
    model_rock_quiet.build_spectra()





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

 .. code-block:: none

    Starting at phase 161.9544290544939 deg, observe for 3.53 h in 26 steps
    Phases = [161.95 163.32 164.68 166.04 167.41 168.77 170.13 171.5  172.86 174.22
     175.59 176.95 178.31 179.68 181.04 182.4  183.77 185.13 186.49 187.86
     189.22 190.58 191.95 193.31 194.67 196.03 197.4 ] deg

    Build Planet:   0%|          | 0/27 [00:00<?, ?it/s]
    Build Planet:   4%|▎         | 1/27 [00:00<00:21,  1.20it/s]
    Build Planet:   7%|▋         | 2/27 [00:01<00:20,  1.20it/s]
    Build Planet:  11%|█         | 3/27 [00:02<00:19,  1.22it/s]
    Build Planet:  15%|█▍        | 4/27 [00:03<00:18,  1.22it/s]
    Build Planet:  19%|█▊        | 5/27 [00:04<00:18,  1.22it/s]
    Build Planet:  22%|██▏       | 6/27 [00:04<00:17,  1.23it/s]
    Build Planet:  26%|██▌       | 7/27 [00:05<00:16,  1.22it/s]
    Build Planet:  30%|██▉       | 8/27 [00:06<00:15,  1.23it/s]
    Build Planet:  33%|███▎      | 9/27 [00:07<00:14,  1.22it/s]
    Build Planet:  37%|███▋      | 10/27 [00:08<00:13,  1.22it/s]
    Build Planet:  41%|████      | 11/27 [00:09<00:13,  1.22it/s]
    Build Planet:  44%|████▍     | 12/27 [00:09<00:12,  1.22it/s]
    Build Planet:  48%|████▊     | 13/27 [00:10<00:11,  1.22it/s]
    Build Planet:  52%|█████▏    | 14/27 [00:11<00:10,  1.22it/s]
    Build Planet:  56%|█████▌    | 15/27 [00:12<00:09,  1.23it/s]
    Build Planet:  59%|█████▉    | 16/27 [00:13<00:09,  1.22it/s]
    Build Planet:  63%|██████▎   | 17/27 [00:13<00:08,  1.21it/s]
    Build Planet:  67%|██████▋   | 18/27 [00:14<00:07,  1.22it/s]
    Build Planet:  70%|███████   | 19/27 [00:15<00:06,  1.22it/s]
    Build Planet:  74%|███████▍  | 20/27 [00:16<00:05,  1.22it/s]
    Build Planet:  78%|███████▊  | 21/27 [00:17<00:04,  1.22it/s]
    Build Planet:  81%|████████▏ | 22/27 [00:18<00:04,  1.22it/s]
    Build Planet:  85%|████████▌ | 23/27 [00:18<00:03,  1.23it/s]
    Build Planet:  89%|████████▉ | 24/27 [00:19<00:02,  1.22it/s]
    Build Planet:  93%|█████████▎| 25/27 [00:20<00:01,  1.23it/s]
    Build Planet:  96%|█████████▋| 26/27 [00:21<00:00,  1.22it/s]
    Build Planet: 100%|██████████| 27/27 [00:22<00:00,  1.22it/s]
    Build Planet: 100%|██████████| 27/27 [00:22<00:00,  1.22it/s]
    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!

    Build Spectra:   0%|          | 0/26 [00:00<?, ?it/s]

    Loading Spectra:   0%|          | 0/12 [00:00<?, ?it/s]

    Loading Spectra:   8%|▊         | 1/12 [00:00<00:03,  3.39it/s]

    Loading Spectra:  17%|█▋        | 2/12 [00:00<00:02,  3.38it/s]

    Loading Spectra:  25%|██▌       | 3/12 [00:00<00:02,  3.38it/s]

    Loading Spectra:  33%|███▎      | 4/12 [00:01<00:02,  3.37it/s]

    Loading Spectra:  42%|████▏     | 5/12 [00:01<00:02,  3.37it/s]

    Loading Spectra:  50%|█████     | 6/12 [00:01<00:01,  3.38it/s]

    Loading Spectra:  58%|█████▊    | 7/12 [00:02<00:01,  3.38it/s]

    Loading Spectra:  67%|██████▋   | 8/12 [00:02<00:01,  3.39it/s]

    Loading Spectra:  75%|███████▌  | 9/12 [00:02<00:00,  3.39it/s]

    Loading Spectra:  83%|████████▎ | 10/12 [00:02<00:00,  3.39it/s]

    Loading Spectra:  92%|█████████▏| 11/12 [00:03<00:00,  3.38it/s]

    Loading Spectra: 100%|██████████| 12/12 [00:03<00:00,  3.38it/s]
    Loading Spectra: 100%|██████████| 12/12 [00:03<00:00,  3.38it/s]

    Build Spectra:   4%|▍         | 1/26 [00:03<01:33,  3.72s/it]
    Build Spectra:   8%|▊         | 2/26 [00:03<00:39,  1.63s/it]
    Build Spectra:  12%|█▏        | 3/26 [00:04<00:22,  1.04it/s]
    Build Spectra:  15%|█▌        | 4/26 [00:04<00:14,  1.53it/s]
    Build Spectra:  19%|█▉        | 5/26 [00:04<00:10,  2.08it/s]
    Build Spectra:  23%|██▎       | 6/26 [00:04<00:07,  2.67it/s]
    Build Spectra:  27%|██▋       | 7/26 [00:04<00:05,  3.25it/s]
    Build Spectra:  31%|███       | 8/26 [00:04<00:04,  3.79it/s]
    Build Spectra:  35%|███▍      | 9/26 [00:05<00:03,  4.27it/s]
    Build Spectra:  38%|███▊      | 10/26 [00:05<00:03,  4.65it/s]
    Build Spectra:  42%|████▏     | 11/26 [00:05<00:03,  4.36it/s]
    Build Spectra:  46%|████▌     | 12/26 [00:05<00:03,  4.20it/s]
    Build Spectra:  50%|█████     | 13/26 [00:06<00:03,  4.08it/s]
    Build Spectra:  54%|█████▍    | 14/26 [00:06<00:03,  4.00it/s]
    Build Spectra:  58%|█████▊    | 15/26 [00:06<00:02,  3.95it/s]
    Build Spectra:  62%|██████▏   | 16/26 [00:06<00:02,  3.93it/s]
    Build Spectra:  65%|██████▌   | 17/26 [00:07<00:02,  3.91it/s]
    Build Spectra:  69%|██████▉   | 18/26 [00:07<00:02,  3.96it/s]
    Build Spectra:  73%|███████▎  | 19/26 [00:07<00:01,  4.39it/s]
    Build Spectra:  77%|███████▋  | 20/26 [00:07<00:01,  4.75it/s]
    Build Spectra:  81%|████████  | 21/26 [00:07<00:00,  5.02it/s]
    Build Spectra:  85%|████████▍ | 22/26 [00:08<00:00,  5.24it/s]
    Build Spectra:  88%|████████▊ | 23/26 [00:08<00:00,  5.41it/s]
    Build Spectra:  92%|█████████▏| 24/26 [00:08<00:00,  5.49it/s]
    Build Spectra:  96%|█████████▌| 25/26 [00:08<00:00,  5.60it/s]
    Build Spectra: 100%|██████████| 26/26 [00:08<00:00,  5.68it/s]
    Build Spectra: 100%|██████████| 26/26 [00:08<00:00,  2.99it/s]




.. GENERATED FROM PYTHON SOURCE LINES 263-268

Load in the data
----------------

We can use VSPEC to read in the synthetic
data we just created.

.. GENERATED FROM PYTHON SOURCE LINES 268-271

.. code-block:: Python


    data_rock_quiet = PhaseAnalyzer(model_rock_quiet.directories['all_model'])





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

 .. code-block:: none

    /opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/VSPEC/analysis.py:116: RuntimeWarning: No Layer info, maybe globes or molecular signatures are off
      warnings.warn(




.. GENERATED FROM PYTHON SOURCE LINES 272-276

Plot the transit
----------------

Lets plot the lightcurve of each channel.

.. GENERATED FROM PYTHON SOURCE LINES 276-312

.. code-block:: Python


    def plot_transit(data:PhaseAnalyzer,title:str,color:str):
        time_from_mid_transit = (data.time-0.5*(data.time[-1]+data.time[0])).to(u.hour)

        fig,axes = plt.subplots(2,1,tight_layout=True)
        axes[0].scatter(
            time_from_mid_transit,
            data.lightcurve('total',(0,-1),normalize=0,noise=False),
            label = 'white light curve',c=color
        )
        axes[0].set_xlabel('Time past mid-transit (hour)')
        axes[0].set_ylabel('Transit depth (ppm)')
        axes[0].legend()
        axes[0].set_title(title)

        # standardize the epochs to use for analysis
        pre_transit = PRE_TRANSIT
        in_transit = IN_TRANSIT

        unocculted_spectrum = data.spectrum('total',pre_transit)
        occulted_spectrum = data.spectrum('total',in_transit)
        lost_to_transit = unocculted_spectrum-occulted_spectrum
        transit_depth = (lost_to_transit/unocculted_spectrum).to_value(u.dimensionless_unscaled)

        axes[1].plot(data.wavelength, 1e6*(transit_depth),c=color)
        axes[1].set_xlabel(f'Wavelength ({data.wavelength.unit})')
        axes[1].set_ylabel('Transit depth (ppm)')
        ylo,yhi = axes[1].get_ylim()
        if yhi-ylo < 0.5:
            mean = 0.5*(ylo+yhi)
            axes[1].set_ylim(mean-0.25,mean+0.25)

        return fig

    plot_transit(data_rock_quiet,'Spotless Star and Bare Rock','xkcd:lavender').show()




.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_spotted_transit_001.png
   :alt: Spotless Star and Bare Rock
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_spotted_transit_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 313-320

Run the other models
--------------------

Let's do the same analysis for the other cases.

Spotless Star, H2O Planet
+++++++++++++++++++++++++

.. GENERATED FROM PYTHON SOURCE LINES 320-329

.. code-block:: Python


    model_h2o_quiet = ObservationModel(params_h2o_quiet)
    model_h2o_quiet.build_planet()
    model_h2o_quiet.build_spectra()

    data_h2o_quiet = PhaseAnalyzer(model_h2o_quiet.directories['all_model'])

    plot_transit(data_h2o_quiet,'Spotless Star and 1 bar H2O Atmosphere','xkcd:azure').show()




.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_spotted_transit_002.png
   :alt: Spotless Star and 1 bar H2O Atmosphere
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_spotted_transit_002.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    Starting at phase 161.9544290544939 deg, observe for 3.53 h in 26 steps
    Phases = [161.95 163.32 164.68 166.04 167.41 168.77 170.13 171.5  172.86 174.22
     175.59 176.95 178.31 179.68 181.04 182.4  183.77 185.13 186.49 187.86
     189.22 190.58 191.95 193.31 194.67 196.03 197.4 ] deg

    Build Planet:   0%|          | 0/27 [00:00<?, ?it/s]
    Build Planet:   4%|▎         | 1/27 [00:01<00:28,  1.11s/it]
    Build Planet:   7%|▋         | 2/27 [00:02<00:28,  1.14s/it]
    Build Planet:  11%|█         | 3/27 [00:03<00:26,  1.11s/it]
    Build Planet:  15%|█▍        | 4/27 [00:04<00:25,  1.11s/it]
    Build Planet:  19%|█▊        | 5/27 [00:05<00:24,  1.10s/it]
    Build Planet:  22%|██▏       | 6/27 [00:06<00:22,  1.09s/it]
    Build Planet:  26%|██▌       | 7/27 [00:07<00:22,  1.12s/it]
    Build Planet:  30%|██▉       | 8/27 [00:08<00:21,  1.11s/it]
    Build Planet:  33%|███▎      | 9/27 [00:10<00:20,  1.14s/it]
    Build Planet:  37%|███▋      | 10/27 [00:11<00:19,  1.12s/it]
    Build Planet:  41%|████      | 11/27 [00:12<00:18,  1.15s/it]
    Build Planet:  44%|████▍     | 12/27 [00:13<00:17,  1.16s/it]
    Build Planet:  48%|████▊     | 13/27 [00:14<00:16,  1.18s/it]
    Build Planet:  52%|█████▏    | 14/27 [00:15<00:15,  1.18s/it]
    Build Planet:  56%|█████▌    | 15/27 [00:17<00:14,  1.19s/it]
    Build Planet:  59%|█████▉    | 16/27 [00:18<00:13,  1.19s/it]
    Build Planet:  63%|██████▎   | 17/27 [00:19<00:11,  1.19s/it]
    Build Planet:  67%|██████▋   | 18/27 [00:20<00:10,  1.18s/it]
    Build Planet:  70%|███████   | 19/27 [00:21<00:09,  1.16s/it]
    Build Planet:  74%|███████▍  | 20/27 [00:22<00:07,  1.14s/it]
    Build Planet:  78%|███████▊  | 21/27 [00:24<00:06,  1.13s/it]
    Build Planet:  81%|████████▏ | 22/27 [00:25<00:05,  1.12s/it]
    Build Planet:  85%|████████▌ | 23/27 [00:26<00:04,  1.11s/it]
    Build Planet:  89%|████████▉ | 24/27 [00:27<00:03,  1.11s/it]
    Build Planet:  93%|█████████▎| 25/27 [00:28<00:02,  1.11s/it]
    Build Planet:  96%|█████████▋| 26/27 [00:29<00:01,  1.10s/it]
    Build Planet: 100%|██████████| 27/27 [00:30<00:00,  1.11s/it]
    Build Planet: 100%|██████████| 27/27 [00:30<00:00,  1.13s/it]
    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!

    Build Spectra:   0%|          | 0/26 [00:00<?, ?it/s]

    Loading Spectra:   0%|          | 0/12 [00:00<?, ?it/s]

    Loading Spectra:   8%|▊         | 1/12 [00:00<00:03,  3.36it/s]

    Loading Spectra:  17%|█▋        | 2/12 [00:00<00:02,  3.35it/s]

    Loading Spectra:  25%|██▌       | 3/12 [00:00<00:02,  3.36it/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.35it/s]

    Loading Spectra:  58%|█████▊    | 7/12 [00:02<00:01,  3.35it/s]

    Loading Spectra:  67%|██████▋   | 8/12 [00:02<00:01,  3.35it/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.35it/s]
    Loading Spectra: 100%|██████████| 12/12 [00:03<00:00,  3.35it/s]

    Build Spectra:   4%|▍         | 1/26 [00:03<01:35,  3.81s/it]
    Build Spectra:   8%|▊         | 2/26 [00:04<00:40,  1.71s/it]
    Build Spectra:  12%|█▏        | 3/26 [00:04<00:23,  1.03s/it]
    Build Spectra:  15%|█▌        | 4/26 [00:04<00:15,  1.39it/s]
    Build Spectra:  19%|█▉        | 5/26 [00:04<00:11,  1.84it/s]
    Build Spectra:  23%|██▎       | 6/26 [00:04<00:08,  2.29it/s]
    Build Spectra:  27%|██▋       | 7/26 [00:05<00:07,  2.70it/s]
    Build Spectra:  31%|███       | 8/26 [00:05<00:05,  3.06it/s]
    Build Spectra:  35%|███▍      | 9/26 [00:05<00:05,  3.37it/s]
    Build Spectra:  38%|███▊      | 10/26 [00:05<00:04,  3.62it/s]
    Build Spectra:  42%|████▏     | 11/26 [00:06<00:04,  3.44it/s]
    Build Spectra:  46%|████▌     | 12/26 [00:06<00:04,  3.28it/s]
    Build Spectra:  50%|█████     | 13/26 [00:06<00:04,  3.20it/s]
    Build Spectra:  54%|█████▍    | 14/26 [00:07<00:03,  3.17it/s]
    Build Spectra:  58%|█████▊    | 15/26 [00:07<00:03,  3.16it/s]
    Build Spectra:  62%|██████▏   | 16/26 [00:07<00:03,  3.14it/s]
    Build Spectra:  65%|██████▌   | 17/26 [00:08<00:02,  3.12it/s]
    Build Spectra:  69%|██████▉   | 18/26 [00:08<00:02,  3.16it/s]
    Build Spectra:  73%|███████▎  | 19/26 [00:08<00:02,  3.42it/s]
    Build Spectra:  77%|███████▋  | 20/26 [00:08<00:01,  3.64it/s]
    Build Spectra:  81%|████████  | 21/26 [00:09<00:01,  3.80it/s]
    Build Spectra:  85%|████████▍ | 22/26 [00:09<00:01,  3.93it/s]
    Build Spectra:  88%|████████▊ | 23/26 [00:09<00:00,  4.02it/s]
    Build Spectra:  92%|█████████▏| 24/26 [00:09<00:00,  4.09it/s]
    Build Spectra:  96%|█████████▌| 25/26 [00:10<00:00,  4.13it/s]
    Build Spectra: 100%|██████████| 26/26 [00:10<00:00,  4.17it/s]
    Build Spectra: 100%|██████████| 26/26 [00:10<00:00,  2.51it/s]




.. GENERATED FROM PYTHON SOURCE LINES 330-332

Spotted Star, CO2 Planet
++++++++++++++++++++++++

.. GENERATED FROM PYTHON SOURCE LINES 332-342

.. code-block:: Python


    model_rock_spotted = ObservationModel(params_rock_spotted)
    model_rock_spotted.build_planet()
    model_rock_spotted.build_spectra()

    data_rock_spotted = PhaseAnalyzer(model_rock_spotted.directories['all_model'])

    plot_transit(data_rock_spotted,'Spotted Star and Bare Rock','xkcd:golden yellow').show()





.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_spotted_transit_003.png
   :alt: Spotted Star and Bare Rock
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_spotted_transit_003.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    Starting at phase 161.9544290544939 deg, observe for 3.53 h in 26 steps
    Phases = [161.95 163.32 164.68 166.04 167.41 168.77 170.13 171.5  172.86 174.22
     175.59 176.95 178.31 179.68 181.04 182.4  183.77 185.13 186.49 187.86
     189.22 190.58 191.95 193.31 194.67 196.03 197.4 ] deg

    Build Planet:   0%|          | 0/27 [00:00<?, ?it/s]
    Build Planet:   4%|▎         | 1/27 [00:00<00:21,  1.20it/s]
    Build Planet:   7%|▋         | 2/27 [00:01<00:20,  1.20it/s]
    Build Planet:  11%|█         | 3/27 [00:02<00:19,  1.22it/s]
    Build Planet:  15%|█▍        | 4/27 [00:03<00:18,  1.22it/s]
    Build Planet:  19%|█▊        | 5/27 [00:04<00:18,  1.22it/s]
    Build Planet:  22%|██▏       | 6/27 [00:04<00:17,  1.22it/s]
    Build Planet:  26%|██▌       | 7/27 [00:05<00:16,  1.22it/s]
    Build Planet:  30%|██▉       | 8/27 [00:06<00:15,  1.22it/s]
    Build Planet:  33%|███▎      | 9/27 [00:07<00:14,  1.22it/s]
    Build Planet:  37%|███▋      | 10/27 [00:08<00:13,  1.22it/s]
    Build Planet:  41%|████      | 11/27 [00:09<00:13,  1.22it/s]
    Build Planet:  44%|████▍     | 12/27 [00:09<00:12,  1.22it/s]
    Build Planet:  48%|████▊     | 13/27 [00:10<00:11,  1.21it/s]
    Build Planet:  52%|█████▏    | 14/27 [00:11<00:10,  1.21it/s]
    Build Planet:  56%|█████▌    | 15/27 [00:12<00:09,  1.22it/s]
    Build Planet:  59%|█████▉    | 16/27 [00:13<00:09,  1.22it/s]
    Build Planet:  63%|██████▎   | 17/27 [00:13<00:08,  1.22it/s]
    Build Planet:  67%|██████▋   | 18/27 [00:14<00:07,  1.22it/s]
    Build Planet:  70%|███████   | 19/27 [00:15<00:06,  1.21it/s]
    Build Planet:  74%|███████▍  | 20/27 [00:16<00:05,  1.21it/s]
    Build Planet:  78%|███████▊  | 21/27 [00:17<00:04,  1.22it/s]
    Build Planet:  81%|████████▏ | 22/27 [00:18<00:04,  1.22it/s]
    Build Planet:  85%|████████▌ | 23/27 [00:18<00:03,  1.22it/s]
    Build Planet:  89%|████████▉ | 24/27 [00:19<00:02,  1.22it/s]
    Build Planet:  93%|█████████▎| 25/27 [00:20<00:01,  1.22it/s]
    Build Planet:  96%|█████████▋| 26/27 [00:21<00:00,  1.22it/s]
    Build Planet: 100%|██████████| 27/27 [00:22<00:00,  1.22it/s]
    Build Planet: 100%|██████████| 27/27 [00:22<00:00,  1.22it/s]
    Generated 94 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!

    Build Spectra:   0%|          | 0/26 [00:00<?, ?it/s]

    Loading Spectra:   0%|          | 0/12 [00:00<?, ?it/s]

    Loading Spectra:   8%|▊         | 1/12 [00:00<00:03,  3.34it/s]

    Loading Spectra:  17%|█▋        | 2/12 [00:00<00:03,  3.32it/s]

    Loading Spectra:  25%|██▌       | 3/12 [00:00<00:02,  3.33it/s]

    Loading Spectra:  33%|███▎      | 4/12 [00:01<00:02,  3.34it/s]

    Loading Spectra:  42%|████▏     | 5/12 [00:01<00:02,  3.33it/s]

    Loading Spectra:  50%|█████     | 6/12 [00:01<00:01,  3.34it/s]

    Loading Spectra:  58%|█████▊    | 7/12 [00:02<00:01,  3.35it/s]

    Loading Spectra:  67%|██████▋   | 8/12 [00:02<00:01,  3.35it/s]

    Loading Spectra:  75%|███████▌  | 9/12 [00:02<00:00,  3.35it/s]

    Loading Spectra:  83%|████████▎ | 10/12 [00:02<00:00,  3.35it/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.35it/s]

    Build Spectra:   4%|▍         | 1/26 [00:03<01:38,  3.93s/it]
    Build Spectra:   8%|▊         | 2/26 [00:04<00:43,  1.82s/it]
    Build Spectra:  12%|█▏        | 3/26 [00:04<00:26,  1.15s/it]
    Build Spectra:  15%|█▌        | 4/26 [00:04<00:18,  1.21it/s]
    Build Spectra:  19%|█▉        | 5/26 [00:05<00:13,  1.53it/s]
    Build Spectra:  23%|██▎       | 6/26 [00:05<00:10,  1.82it/s]
    Build Spectra:  27%|██▋       | 7/26 [00:05<00:09,  2.08it/s]
    Build Spectra:  31%|███       | 8/26 [00:06<00:07,  2.29it/s]
    Build Spectra:  35%|███▍      | 9/26 [00:06<00:06,  2.45it/s]
    Build Spectra:  38%|███▊      | 10/26 [00:07<00:06,  2.58it/s]
    Build Spectra:  42%|████▏     | 11/26 [00:07<00:06,  2.48it/s]
    Build Spectra:  46%|████▌     | 12/26 [00:07<00:05,  2.41it/s]
    Build Spectra:  50%|█████     | 13/26 [00:08<00:05,  2.38it/s]
    Build Spectra:  54%|█████▍    | 14/26 [00:08<00:05,  2.37it/s]
    Build Spectra:  58%|█████▊    | 15/26 [00:09<00:04,  2.35it/s]
    Build Spectra:  62%|██████▏   | 16/26 [00:09<00:04,  2.33it/s]
    Build Spectra:  65%|██████▌   | 17/26 [00:10<00:03,  2.33it/s]
    Build Spectra:  69%|██████▉   | 18/26 [00:10<00:03,  2.35it/s]
    Build Spectra:  73%|███████▎  | 19/26 [00:10<00:02,  2.50it/s]
    Build Spectra:  77%|███████▋  | 20/26 [00:11<00:02,  2.61it/s]
    Build Spectra:  81%|████████  | 21/26 [00:11<00:01,  2.70it/s]
    Build Spectra:  85%|████████▍ | 22/26 [00:11<00:01,  2.75it/s]
    Build Spectra:  88%|████████▊ | 23/26 [00:12<00:01,  2.80it/s]
    Build Spectra:  92%|█████████▏| 24/26 [00:12<00:00,  2.84it/s]
    Build Spectra:  96%|█████████▌| 25/26 [00:12<00:00,  2.86it/s]
    Build Spectra: 100%|██████████| 26/26 [00:13<00:00,  2.88it/s]
    Build Spectra: 100%|██████████| 26/26 [00:13<00:00,  1.97it/s]
    /opt/hostedtoolcache/Python/3.11.10/x64/lib/python3.11/site-packages/VSPEC/analysis.py:116: RuntimeWarning: No Layer info, maybe globes or molecular signatures are off
      warnings.warn(




.. GENERATED FROM PYTHON SOURCE LINES 343-347

Plot the observed spectra
-------------------------

Let's compare the transits. We also load in the actual JWST data.

.. GENERATED FROM PYTHON SOURCE LINES 349-354

Get the data
++++++++++++

Reduced data from :cite:t:`2023ApJ...948L..11M` is publicly available.
However, you must download it from the figure caption of the online version.

.. GENERATED FROM PYTHON SOURCE LINES 354-370

.. code-block:: Python


    import pandas as pd

    try:
        filename = Path(__file__).parent / 'moran2023_fig3.txt'
    except NameError:
        filename = 'moran2023_fig3.txt'

    df = pd.read_fwf(filename,colspecs=[(0,8),(9,14),(15,20),(21,25),(26,28)],
        header=20,names=['Reduction','Wave','Width','Depth','e_Depth'])
    used_eureka = df['Reduction']=='Eureka'
    moran_x = df.loc[used_eureka,'Wave']
    moran_y = df.loc[used_eureka,'Depth']
    moran_yerr = df.loc[used_eureka,'e_Depth']
    moran_mean = np.mean(moran_y)








.. GENERATED FROM PYTHON SOURCE LINES 371-374

Make the figure
+++++++++++++++


.. GENERATED FROM PYTHON SOURCE LINES 374-401

.. code-block:: Python


    fig, ax = plt.subplots(1,1)

    for data,label,color in zip(
        [data_rock_quiet,data_h2o_quiet,data_rock_spotted],
        ['Rock', 'H2O', 'Rock+Spots'],
        ['xkcd:lavender','xkcd:azure','xkcd:golden yellow']
    ):
        pre_transit = PRE_TRANSIT
        in_transit = IN_TRANSIT
        unocculted_spectrum = data.spectrum('total',pre_transit)
        occulted_spectrum = data.spectrum('total',in_transit)
        lost_to_transit = unocculted_spectrum-occulted_spectrum
        transit_depth = (lost_to_transit/unocculted_spectrum).to_value(u.dimensionless_unscaled)*1e6
        depth_mean = np.mean(transit_depth)
        shift = moran_mean/depth_mean
        shift = 1
        ax.plot(data.wavelength,transit_depth*shift,label=label,color=color)

    ax.errorbar(moran_x,moran_y,yerr=moran_yerr,
        fmt='o',color='k',label='Moran+23 Eureka!')

    ax.set_xlabel('Wavelength (um)')
    ax.set_ylabel('Transit depth (ppm)')
    ax.legend()





.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_spotted_transit_004.png
   :alt: plot spotted transit
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_spotted_transit_004.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x7fd410499c50>



.. GENERATED FROM PYTHON SOURCE LINES 402-406

Make the figure again
+++++++++++++++++++++

This time without the bare rock

.. GENERATED FROM PYTHON SOURCE LINES 406-429

.. code-block:: Python


    fig, ax = plt.subplots(1,1,figsize=(5,3.5),tight_layout=True)

    for data,label,color in zip(
        [data_h2o_quiet,data_rock_spotted],
        ['H2O', 'Rock+Spots'],
        ['xkcd:azure','xkcd:golden yellow']
    ):
        pre_transit = PRE_TRANSIT
        in_transit = IN_TRANSIT
        unocculted_spectrum = data.spectrum('total',pre_transit)
        occulted_spectrum = data.spectrum('total',in_transit)
        lost_to_transit = unocculted_spectrum-occulted_spectrum
        transit_depth = (lost_to_transit/unocculted_spectrum).to_value(u.dimensionless_unscaled)
        ax.plot(data.wavelength,transit_depth*1e6,label=label,color=color)

    ax.errorbar(moran_x,moran_y,yerr=moran_yerr,
        fmt='o',color='k',label='Moran+23 Eureka!')

    ax.set_xlabel('Wavelength (um)')
    ax.set_ylabel('Transit depth (ppm)')
    ax.legend()




.. image-sg:: /auto_examples/end_to_end/images/sphx_glr_plot_spotted_transit_005.png
   :alt: plot spotted transit
   :srcset: /auto_examples/end_to_end/images/sphx_glr_plot_spotted_transit_005.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x7fd418530f90>




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

   **Total running time of the script:** (2 minutes 40.511 seconds)


.. _sphx_glr_download_auto_examples_end_to_end_plot_spotted_transit.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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