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

.. only:: html

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

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

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

.. _sphx_glr_nesc_test_cases_nesc_case03.py:


====================================================
Case 3: Tumbling brick with dynamic damping, no drag
====================================================

==============  ===============
Verifies        Inertial coupling [dynamic damping model]
Gravitation     J2
Geodesy         WGS-84 rotating
Atmosphere      US 1976 STD
Winds           still air
Vehicle         Dragless rotating brick with aero damping
Notes           Drag coefficient set to zero
==============  ===============

.. GENERATED FROM PYTHON SOURCE LINES 16-47

.. code-block:: Python


    from simupy.block_diagram import BlockDiagram
    import simupy_flight
    import numpy as np

    from nesc_testcase_helper import plot_nesc_comparisons, int_opts, benchmark
    from nesc_testcase_helper import ft_per_m, kg_per_slug

    planet = simupy_flight.Planet(
        gravity=simupy_flight.earth_J2_gravity,
        winds=simupy_flight.get_constant_winds(),
        atmosphere=simupy_flight.atmosphere_1976,
        planetodetics=simupy_flight.Planetodetic(
            a=simupy_flight.earth_equitorial_radius,
            omega_p=simupy_flight.earth_rotation_rate,
            f=simupy_flight.earth_f,
        ),
    )

    Ixx = 0.001894220 * kg_per_slug / (ft_per_m**2)  # slug-ft2
    Iyy = 0.006211019 * kg_per_slug / (ft_per_m**2)  # slug-ft2
    Izz = 0.007194665 * kg_per_slug / (ft_per_m**2)  # slug-ft2
    Ixy = 0.0 * kg_per_slug / (ft_per_m**2)  # slug-ft2
    Iyz = 0.0 * kg_per_slug / (ft_per_m**2)  # slug-ft2
    Izx = 0.0 * kg_per_slug / (ft_per_m**2)  # slug-ft2
    m = 0.155404754 * kg_per_slug  # slug

    x = 0.0
    y = 0.0
    z = 0.0








.. GENERATED FROM PYTHON SOURCE LINES 48-49

Configure the constant dynamic damping

.. GENERATED FROM PYTHON SOURCE LINES 49-103

.. code-block:: Python


    S_A = 0.22222 / (ft_per_m**2)
    b_l = 1 / (3 * ft_per_m)
    c_l = 2 / (3 * ft_per_m)
    a_l = b_l
    aero_model = simupy_flight.get_constant_aero(Cp_b=-1.0, Cq_b=-1.0, Cr_b=-1.0)
    vehicle = simupy_flight.Vehicle(
        base_aero_coeffs=aero_model,
        m=m,
        I_xx=Ixx,
        I_yy=Iyy,
        I_zz=Izz,
        I_xy=Ixy,
        I_yz=Iyz,
        I_xz=Izx,
        x_com=x,
        y_com=y,
        z_com=z,
        x_mrc=x,
        y_mrc=y,
        z_mrc=z,
        S_A=S_A,
        a_l=a_l,
        b_l=b_l,
        c_l=c_l,
        d_l=0.0,
    )

    BD = BlockDiagram(planet, vehicle)
    BD.connect(planet, vehicle, inputs=np.arange(planet.dim_output))
    BD.connect(vehicle, planet, inputs=np.arange(vehicle.dim_output))

    lat_ic = 0.0 * np.pi / 180
    long_ic = 0.0 * np.pi / 180
    h_ic = 30_000 / ft_per_m
    V_N_ic = 0.0
    V_E_ic = 0.0
    V_D_ic = 0.0
    psi_ic = 0.0 * np.pi / 180
    theta_ic = 0.0 * np.pi / 180
    phi_ic = 0.0 * np.pi / 180
    omega_X_ic = 10.0 * np.pi / 180
    omega_Y_ic = 20.0 * np.pi / 180
    omega_Z_ic = 30.0 * np.pi / 180

    planet.initial_condition = planet.ic_from_planetodetic(
        long_ic, lat_ic, h_ic, V_N_ic, V_E_ic, V_D_ic, psi_ic, theta_ic, phi_ic
    )
    planet.initial_condition[-3:] = omega_X_ic, omega_Y_ic, omega_Z_ic

    with benchmark() as b:
        res = BD.simulate(30, integrator_options=int_opts)

    plot_nesc_comparisons(res, "03")



.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /nesc_test_cases/images/sphx_glr_nesc_case03_001.png
         :alt: nesc case03
         :srcset: /nesc_test_cases/images/sphx_glr_nesc_case03_001.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /nesc_test_cases/images/sphx_glr_nesc_case03_002.png
         :alt: nesc case03
         :srcset: /nesc_test_cases/images/sphx_glr_nesc_case03_002.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /nesc_test_cases/images/sphx_glr_nesc_case03_003.png
         :alt: nesc case03
         :srcset: /nesc_test_cases/images/sphx_glr_nesc_case03_003.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /nesc_test_cases/images/sphx_glr_nesc_case03_004.png
         :alt: nesc case03
         :srcset: /nesc_test_cases/images/sphx_glr_nesc_case03_004.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /nesc_test_cases/images/sphx_glr_nesc_case03_005.png
         :alt: nesc case03
         :srcset: /nesc_test_cases/images/sphx_glr_nesc_case03_005.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /nesc_test_cases/images/sphx_glr_nesc_case03_006.png
         :alt: nesc case03
         :srcset: /nesc_test_cases/images/sphx_glr_nesc_case03_006.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /nesc_test_cases/images/sphx_glr_nesc_case03_007.png
         :alt: nesc case03
         :srcset: /nesc_test_cases/images/sphx_glr_nesc_case03_007.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /nesc_test_cases/images/sphx_glr_nesc_case03_008.png
         :alt: nesc case03
         :srcset: /nesc_test_cases/images/sphx_glr_nesc_case03_008.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /nesc_test_cases/images/sphx_glr_nesc_case03_009.png
         :alt: nesc case03
         :srcset: /nesc_test_cases/images/sphx_glr_nesc_case03_009.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /nesc_test_cases/images/sphx_glr_nesc_case03_010.png
         :alt: nesc case03
         :srcset: /nesc_test_cases/images/sphx_glr_nesc_case03_010.png
         :class: sphx-glr-multi-img


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

 .. code-block:: none

    /home/runner/work/simupy-flight/simupy-flight/simupy_flight/kinematics.py:198: RuntimeWarning: invalid value encountered in scalar divide
      x60 = x6*(2*x26*x54 + x58*(-x13 + x14) + x59*(x33 + x55))/x53
    /home/runner/work/simupy-flight/simupy-flight/simupy_flight/dynamics.py:22: RuntimeWarning: divide by zero encountered in scalar divide
      x15 = numpy.select([numpy.greater(V_T, 0.0)], [(1/2)/V_T], default=0.0)
    /home/runner/work/simupy-flight/simupy-flight/simupy_flight/kinematics.py:198: RuntimeWarning: divide by zero encountered in scalar divide
      x60 = x6*(2*x26*x54 + x58*(-x13 + x14) + x59*(x33 + x55))/x53
    time to simulate: 1.686 s





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

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


.. _sphx_glr_download_nesc_test_cases_nesc_case03.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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