.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/04-topography/4.4-plot_ray_intensity.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_04-topography_4.4-plot_ray_intensity.py>`
        to download the full example code.

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

.. _sphx_glr_auto_examples_04-topography_4.4-plot_ray_intensity.py:


Plot the ray intensity map for a crater
=======================================

.. rubric:: By David Minton

This example demonstrates how to create a ray intensity map for a crater using the "basicmoon" morphology model.

The ray intensity model simulates the spatial modulation of ejecta distribution in the form of rays, a feature observed around many planetary impact craters. These rays are modeled using a layered pattern of radially aligned Gaussian contributions, whose intensity decays with radial distance from the crater rim. The number and angular spread of rays vary stochastically to produce a realistic, azimuthally varying intensity field.

The model uses the ray intensity function compiled from Rust code for performance. This function computes the contribution of each ray at a given point based on its radial distance and angular bearing from the crater center. The output intensity is plotted in a 2D space normalized to the crater radius and visualized on a logarithmic color scale to emphasize the structure of the rays.

Note, that this example is rather more complex than it really needs to be. In practice, the Morphology object is always associated with a Surface object, which contains its own set of methods for visualizing the surface morphology. For this example, we are bypassing the Surface object functionality entirely and generating our own grid for illustration purposes. Example 1.1-visualize_one_crater.py demonstrates a more "typical" approach.

.. GENERATED FROM PYTHON SOURCE LINES 16-68



.. image-sg:: /auto_examples/04-topography/images/sphx_glr_4.4-plot_ray_intensity_001.png
   :alt: Ray Intensity Map
   :srcset: /auto_examples/04-topography/images/sphx_glr_4.4-plot_ray_intensity_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    Creating a new grid
    Generating a mesh with icosphere level 4.
    Done






|

.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np
    from matplotlib.colors import LogNorm

    from cratermaker import Crater, Morphology

    simdir = "simdata-4_4"

    # Note, that for these examples we pass ask_overwrite=False and reset=True to the Simulation constructor. This will suppress
    # prompts that ask the user if they want to overwrite existing files, which would would prevent these examples from running on their
    # own when building the documentation pages. Alternatively, calling cm.cleanup(simdir) will remove all pre-existing output files.

    crater = Crater.maker(radius=10.0e3)
    # Because we are not explicitly passing a Surface object, the Morphology constructor will generate a default surface. We pass the "simdir" and "gridlevel" arguments to control the Surface generation, even though we don't make use of it directly here.
    morphology = Morphology.maker(simdir=simdir, gridlevel=4, ask_overwrite=False, reset=True)
    rc = crater.radius

    grid_size = 1000
    rmax = 20
    morphology.ejecta_truncation = rmax
    extent = rmax * rc
    x = np.linspace(-extent, extent, grid_size)
    y = np.linspace(-extent, extent, grid_size)
    xx, yy = np.meshgrid(x, y)

    # Convert Cartesian to polar coordinates
    r = np.sqrt(xx**2 + yy**2)
    theta = np.degrees(np.arctan2(yy, xx)) + 180.0

    # Compute ray intensity
    r_flat = r.ravel()
    theta_flat = theta.ravel()
    intensity_flat = morphology.ray_intensity(crater, r_flat, theta_flat)
    intensity = intensity_flat.reshape(r.shape)
    # Plot the intensity colormap
    fig, ax = plt.subplots(figsize=(6, 5))
    c = ax.imshow(
        intensity,
        norm=LogNorm(vmin=1e-4, vmax=intensity.max()),
        origin="lower",
        extent=[-rmax, rmax, -rmax, rmax],
        cmap="inferno",
    )
    ax.set_xlabel("x (r_c)")
    ax.set_ylabel("y (r_c)")
    ax.set_title("Ray Intensity Map")
    fig.colorbar(c, ax=ax, label="Intensity")

    plt.tight_layout()
    plt.show()
    print("Done")


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

   **Total running time of the script:** (0 minutes 4.722 seconds)


.. _sphx_glr_download_auto_examples_04-topography_4.4-plot_ray_intensity.py:

.. only:: html

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

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

      :download:`Download Jupyter notebook: 4.4-plot_ray_intensity.ipynb <4.4-plot_ray_intensity.ipynb>`

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

      :download:`Download Python source code: 4.4-plot_ray_intensity.py <4.4-plot_ray_intensity.py>`

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

      :download:`Download zipped: 4.4-plot_ray_intensity.zip <4.4-plot_ray_intensity.zip>`


.. only:: html

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

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