Space Shuttle Reentry Trajectory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Space Shuttle Reentry Trajectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Pkg\n",
    "Pkg.activate(joinpath(@__DIR__))\n",
    "Pkg.instantiate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Interpolations\n",
    "using KNITRO\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "using SparseArrays\n",
    "using SparseDiffTools\n",
    "using SparsityDetection\n",
    "\n",
    "using BenchmarkTools\n",
    "using ForwardDiff\n",
    "using ReverseDiff\n",
    "using Random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Global variables\n",
    "const w  = 203000  # weight (lb)\n",
    "const g₀ = 32.174  # acceleration (ft/sec^2)\n",
    "const m  = w / g₀  # mass (slug)\n",
    "\n",
    "# Aerodynamic and atmospheric forces on the vehicle\n",
    "const ρ₀ =  0.002378\n",
    "const hᵣ =  23800.0\n",
    "const Rₑ =  20902900.0\n",
    "const μ  =  0.14076539e17\n",
    "const S  =  2690.0\n",
    "const a₀ = -0.20704\n",
    "const a₁ =  0.029244\n",
    "const b₀ =  0.07854\n",
    "const b₁ = -0.61592e-2\n",
    "const b₂ =  0.621408e-3\n",
    "const c₀ =  1.0672181\n",
    "const c₁ = -0.19213774e-1\n",
    "const c₂ =  0.21286289e-3\n",
    "const c₃ = -0.10117249e-5\n",
    "\n",
    "c_L(a) = a₀ + a₁ * rad2deg(a)\n",
    "c_D(a) = b₀ + b₁ * rad2deg(a) + b₂ * rad2deg(a)^2\n",
    "\n",
    "D(h, v, a) = 0.5 * c_D(a) * S * ρ(h) * v^2\n",
    "L(h, v, a) = 0.5 * c_L(a) * S * ρ(h) * v^2\n",
    "\n",
    "r(h) = Rₑ + h\n",
    "g(h) = μ / r(h)^2\n",
    "ρ(h) = ρ₀ * exp(-h / hᵣ)\n",
    "\n",
    "# Aerodynamic heating on the vehicle wing leading edge\n",
    "qₐ(a) = c₀ + c₁ * rad2deg(a) + c₂ * rad2deg(a)^2 + c₃ * rad2deg(a)^3\n",
    "qᵣ(h, v) = 17700 * √ρ(h) * (0.0001 * v)^3.07\n",
    "q(h, v, a) = qₐ(a) * qᵣ(h, v)\n",
    "\n",
    "# Motion of the vehicle as a differential-algebraic system of equations (DAEs)\n",
    "δh(v, γ) = v * sin(γ)\n",
    "δϕ(h, θ, v, γ, ψ) = (v / r(h)) * cos(γ) * sin(ψ) / cos(θ)\n",
    "δθ(h, v, γ, ψ) = (v / r(h)) * cos(γ) * cos(ψ)\n",
    "δv(h, v, γ, α) = -(D(h, v, α) / m) - g(h) * sin(γ)\n",
    "δγ(h, v, γ, α, β) = (L(h, v, α) / (m * v)) * cos(β) + cos(γ) * ((v / r(h)) - (g(h) / v))\n",
    "δψ(h, θ, v, γ, ψ, α, β) = (1 / (m * v * cos(γ))) * L(h, v, α) * sin(β) + (v / (r(h) * cos(θ))) * cos(γ) * sin(ψ) * sin(θ)\n",
    "\n",
    "# Initial conditions\n",
    "hₛ = 2.6          # altitude (ft) / 1e5\n",
    "ϕₛ = deg2rad(0)   # longitude (rad)\n",
    "θₛ = deg2rad(0)   # latitude (rad)\n",
    "vₛ = 2.56         # velocity (ft/sec) / 1e4\n",
    "γₛ = deg2rad(-1)  # flight path angle (rad)\n",
    "ψₛ = deg2rad(90)  # azimuth (rad)\n",
    "αₛ = deg2rad(0)   # angle of attack (rad)\n",
    "βₛ = deg2rad(0)   # bank angle (rad)\n",
    "tₛ = 1.00         # time step (sec)\n",
    "\n",
    "# Final conditions, the so-called Terminal Area Energy Management (TAEM)\n",
    "hₜ = 0.8          # altitude (ft) / 1e5\n",
    "vₜ = 0.25         # velocity (ft/sec) / 1e4\n",
    "γₜ = deg2rad(-5)  # flight path angle (rad)\n",
    "\n",
    "xₛ = [hₛ, ϕₛ, θₛ, vₛ, γₛ, ψₛ, αₛ, βₛ, tₛ]\n",
    "xₜ = [hₜ, ϕₛ, θₛ, vₜ, γₜ, ψₛ, αₛ, βₛ, tₛ]\n",
    "\n",
    "@show xₛ xₜ;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Direct Transcription Formulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 2008  # number of segments\n",
    "M = N + 1  # number of mesh points\n",
    "\n",
    "# Linear initial guess between the boundary conditions\n",
    "interp_linear = LinearInterpolation([1, M], [xₛ, xₜ])\n",
    "initial_guess = [value for x = 1:M for value in interp_linear(x)]\n",
    "\n",
    "nₓ = length(xₛ)  # mesh points dimension\n",
    "n = nₓ * M  # number of decision variables\n",
    "@assert n == length(initial_guess)\n",
    "\n",
    "ind_h = 1:nₓ:n  # indices of the `h` decision variables\n",
    "ind_ϕ = 2:nₓ:n  # indices of the `ϕ` decision variables\n",
    "ind_θ = 3:nₓ:n  # indices of the `θ` decision variables\n",
    "ind_v = 4:nₓ:n  # indices of the `v` decision variables\n",
    "ind_γ = 5:nₓ:n  # indices of the `γ` decision variables\n",
    "ind_ψ = 6:nₓ:n  # indices of the `ψ` decision variables\n",
    "ind_α = 7:nₓ:n  # indices of the `α` decision variables\n",
    "ind_β = 8:nₓ:n  # indices of the `β` decision variables\n",
    "ind_t = 9:nₓ:n  # indices of the `t` decision variables\n",
    "@assert length(ind_h) == length(ind_v) == length(ind_α) == M\n",
    "\n",
    "m_dyn = 6 * N\n",
    "ind_con_dyn = 1:m_dyn  # indices of the dynamics' constraints (defects)\n",
    "@assert length(ind_con_dyn) == 6 * (M - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forward dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EulerStep(x, δx, h) = x + h * δx\n",
    "\n",
    "function dynamics_defects!(dx, x)\n",
    "    hᵢ, ϕᵢ, θᵢ, vᵢ, γᵢ, ψᵢ, αᵢ, βᵢ, tᵢ, hⱼ, ϕⱼ, θⱼ, vⱼ, γⱼ, ψⱼ = x\n",
    "\n",
    "    # Unit correction due to decision variable scaling\n",
    "    hᵢ *= 1e5; hⱼ *= 1e5; vᵢ *= 1e4; vⱼ *= 1e4\n",
    "\n",
    "    dx[1] = EulerStep(hᵢ, δh(        vᵢ, γᵢ            ), tᵢ) - hⱼ\n",
    "    dx[2] = EulerStep(ϕᵢ, δϕ(hᵢ, θᵢ, vᵢ, γᵢ, ψᵢ        ), tᵢ) - ϕⱼ\n",
    "    dx[3] = EulerStep(θᵢ, δθ(hᵢ,     vᵢ, γᵢ, ψᵢ        ), tᵢ) - θⱼ\n",
    "    dx[4] = EulerStep(vᵢ, δv(hᵢ,     vᵢ, γᵢ,     αᵢ    ), tᵢ) - vⱼ\n",
    "    dx[5] = EulerStep(γᵢ, δγ(hᵢ,     vᵢ, γᵢ,     αᵢ, βᵢ), tᵢ) - γⱼ\n",
    "    dx[6] = EulerStep(ψᵢ, δψ(hᵢ, θᵢ, vᵢ, γᵢ, ψᵢ, αᵢ, βᵢ), tᵢ) - ψⱼ\n",
    "\n",
    "    return dx\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "struct JacobianData{T1 <: SparseMatrixCSC,T2 <: ForwardColorJacCache}\n",
    "    jac_dyn::T1\n",
    "    jac_cache::T2\n",
    "    length_jac::Int\n",
    "\n",
    "    function JacobianData(; input = rand(15), output = zeros(6))\n",
    "        sparsity_pattern = jacobian_sparsity(dynamics_defects!, output, input)\n",
    "        jac_dyn = convert.(Float64, sparse(sparsity_pattern))\n",
    "\n",
    "        jac_cache = ForwardColorJacCache(dynamics_defects!, input, dx = output,\n",
    "                                         colorvec = matrix_colors(jac_dyn),\n",
    "                                         sparsity = sparsity_pattern)\n",
    "\n",
    "        length_jac = nnz(jac_dyn)\n",
    "\n",
    "        T1 = typeof(jac_dyn)\n",
    "        T2 = typeof(jac_cache)\n",
    "\n",
    "        new{T1,T2}(jac_dyn, jac_cache, length_jac)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input = rand(15)\n",
    "jac_data = JacobianData(input=input)\n",
    "forwarddiff_color_jacobian!(jac_data.jac_dyn, dynamics_defects!, input, jac_data.jac_cache)\n",
    "\n",
    "jacIndexConsCB = Cint.(hcat([rowvals(jac_data.jac_dyn) .+ i*6 for i = 0:N-1]...) .- 1)\n",
    "jacIndexVarsCB = Cint.(hcat([vcat([fill(j + i*nₓ, length(nzrange(jac_data.jac_dyn, j))) for j = 1:15]...) for i = 0:N-1]...) .- 1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hessian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function hessian(f!, x, lambda)\n",
    "    function extended_gradient(x::AbstractArray{T}) where {T}\n",
    "        dx = zeros(T, 6)\n",
    "        ForwardDiff.jacobian(f!, dx, x)' * lambda\n",
    "    end\n",
    "\n",
    "    ForwardDiff.jacobian(extended_gradient, x)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hessIndexVars1CB = nothing\n",
    "hessIndexVars2CB = nothing\n",
    "\n",
    "hess_cis = nothing\n",
    "\n",
    "let\n",
    "    global hessIndexVars1CB, hessIndexVars2CB, hess_cis\n",
    "\n",
    "    num_samples = 50\n",
    "    samples = [hessian(dynamics_defects!, rand(15), rand(6))\n",
    "               for _ in 1:num_samples]  # Workaround for sparsity detection\n",
    "\n",
    "    hess = samples |> sum |> triu |> sparse\n",
    "    @assert nnz(hess) == 31\n",
    "\n",
    "    ind_hess_vars1 = rowvals(hess)\n",
    "    ind_hess_vars2 = vcat([fill(j, length(nzrange(hess, j)))\n",
    "                           for j = 1:size(hess, 2)]...)\n",
    "\n",
    "    println(\"ind_hess_vars1 = \", ind_hess_vars1)\n",
    "    println(\"ind_hess_vars2 = \", ind_hess_vars2)\n",
    "\n",
    "    hess_cis = CartesianIndex.(ind_hess_vars1, ind_hess_vars2)\n",
    "\n",
    "    hessIndexVars1CB = Cint.(hcat([ind_hess_vars1 .+ i * nₓ for i = 0:M - 2]...) .- 1)\n",
    "    hessIndexVars2CB = Cint.(hcat([ind_hess_vars2 .+ i * nₓ for i = 0:M - 2]...) .- 1)\n",
    "end\n",
    "\n",
    "@assert length(hessIndexVars1CB) == N * 31\n",
    "@assert length(hessIndexVars1CB) == length(hessIndexVars2CB)\n",
    "@assert eltype(hessIndexVars1CB) == eltype(hessIndexVars2CB) == Cint\n",
    "@assert length(hess_cis) == 31"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Callbacks to evaluate the constraints and their derivatives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function cb_eval_fc_con_dyn(kc, cb, evalRequest, evalResult, userParams)\n",
    "    x = evalRequest.x\n",
    "\n",
    "    for i = 0:M - 2\n",
    "        ind_xᵢ = (1:nₓ) .+ i * nₓ\n",
    "        ind_xⱼ = ind_xᵢ .+ nₓ\n",
    "        ind_con = (1:6) .+ i * 6\n",
    "\n",
    "        xᵢ = x[ind_xᵢ][1:9]\n",
    "        xⱼ = x[ind_xⱼ][1:6]\n",
    "\n",
    "        @views dynamics_defects!(evalResult.c[ind_con], [xᵢ; xⱼ])\n",
    "    end\n",
    "\n",
    "    return 0\n",
    "end\n",
    "\n",
    "function cb_eval_ga_con_dyn(kc, cb, evalRequest, evalResult, userParams)\n",
    "    x = evalRequest.x\n",
    "    jac_data = userParams\n",
    "\n",
    "    for i = 0:M - 2\n",
    "        ind_xᵢ = (1:nₓ) .+ i * nₓ\n",
    "        ind_xⱼ = ind_xᵢ .+ nₓ\n",
    "        ind_con = (1:6) .+ i * 6\n",
    "\n",
    "        xᵢ = x[ind_xᵢ][1:9]\n",
    "        xⱼ = x[ind_xⱼ][1:6]\n",
    "\n",
    "        forwarddiff_color_jacobian!(jac_data.jac_dyn, dynamics_defects!, [xᵢ; xⱼ], jac_data.jac_cache)\n",
    "\n",
    "        aux_jac = i * jac_data.length_jac\n",
    "        ind_jac = aux_jac + 1 : aux_jac + jac_data.length_jac\n",
    "        evalResult.jac[ind_jac] = nonzeros(jac_data.jac_dyn)\n",
    "    end\n",
    "\n",
    "    return 0\n",
    "end\n",
    "\n",
    "function cb_eval_h_con_dyn(kc, cb, evalRequest, evalResult, userParams)\n",
    "    x = evalRequest.x\n",
    "    lambda = evalRequest.lambda\n",
    "    sigma = evalRequest.sigma\n",
    "\n",
    "    for i = 0:M - 2\n",
    "        ind_xᵢ = (1:nₓ) .+ i * nₓ\n",
    "        ind_xⱼ = ind_xᵢ .+ nₓ\n",
    "\n",
    "        xᵢ = x[ind_xᵢ][1:9]\n",
    "        xⱼ = x[ind_xⱼ][1:6]\n",
    "\n",
    "        ind_lambda = (1:6) .+ i * 6  # .+ n  # `n`: number of decision variables\n",
    "\n",
    "        hess = hessian(dynamics_defects!, [xᵢ; xⱼ], lambda[ind_lambda])\n",
    "\n",
    "        ind_hess = (1:31) .+ i * 31\n",
    "\n",
    "        evalResult.hess[ind_hess] = hess[hess_cis]\n",
    "    end\n",
    "\n",
    "    return 0\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lm = KNITRO.LMcontext()\n",
    "kc = KNITRO.KN_new_lm(lm)\n",
    "\n",
    "KNITRO.KN_add_vars(kc, n)\n",
    "\n",
    "!isempty(initial_guess) && KNITRO.KN_set_var_primal_init_values_all(kc, initial_guess)\n",
    "\n",
    "KNITRO.KN_set_var_lobnds(kc, M, collect(Cint, ind_h .- 1), zeros(M))  # `h` bounds\n",
    "KNITRO.KN_set_var_lobnds(kc, M, collect(Cint, ind_v .- 1), fill(1 / 1e4, M))  # `v` bounds\n",
    "\n",
    "KNITRO.KN_set_var_lobnds(kc, M, collect(Cint, ind_θ .- 1), fill(deg2rad(-89), M))  # `θ` lower bounds\n",
    "KNITRO.KN_set_var_upbnds(kc, M, collect(Cint, ind_θ .- 1), fill(deg2rad( 89), M))  # `θ` upper bounds\n",
    "\n",
    "KNITRO.KN_set_var_lobnds(kc, M, collect(Cint, ind_γ .- 1), fill(deg2rad(-89), M))  # `γ` lower bounds\n",
    "KNITRO.KN_set_var_upbnds(kc, M, collect(Cint, ind_γ .- 1), fill(deg2rad( 89), M))  # `γ` upper bounds\n",
    "\n",
    "KNITRO.KN_set_var_lobnds(kc, M, collect(Cint, ind_α .- 1), fill(deg2rad(-90), M))  # `α` lower bounds\n",
    "KNITRO.KN_set_var_upbnds(kc, M, collect(Cint, ind_α .- 1), fill(deg2rad( 90), M))  # `α` upper bounds\n",
    "\n",
    "KNITRO.KN_set_var_lobnds(kc, M, collect(Cint, ind_β .- 1), fill(deg2rad(-89), M))  # `β` lower bounds\n",
    "KNITRO.KN_set_var_upbnds(kc, M, collect(Cint, ind_β .- 1), fill(deg2rad(  1), M))  # `β` upper bounds\n",
    "\n",
    "# KNITRO.KN_set_var_lobnds(kc, collect(Cint, ind_t .- 1), fill(0.50, M))  # `t` lower bounds\n",
    "KNITRO.KN_set_var_fxbnds(kc, M, collect(Cint, ind_t .- 1), fill(1.00, M))  # Fix time steps\n",
    "\n",
    "# Fix initial and final conditions\n",
    "ind_fixed_vars = [ind_h[1]:ind_ψ[1] ; [ind_h[end], ind_v[end], ind_γ[end]]]\n",
    "val_fixed_vars = [hₛ, ϕₛ, θₛ, vₛ, γₛ, ψₛ, hₜ, vₜ, γₜ]\n",
    "KNITRO.KN_set_var_fxbnds(kc, length(val_fixed_vars), collect(Cint, ind_fixed_vars .- 1), val_fixed_vars)\n",
    "\n",
    "KNITRO.KN_add_cons(kc, m_dyn)\n",
    "\n",
    "KNITRO.KN_set_con_eqbnds(kc, m_dyn, collect(Cint, ind_con_dyn .- 1), zeros(m_dyn))  # defects\n",
    "\n",
    "# This callback does not evaluate the objective function. As such,\n",
    "# we pass `false` as the second argument to `KN_add_eval_callback()`.\n",
    "cb_dyn = KNITRO.KN_add_eval_callback(kc, false, collect(Cint, ind_con_dyn .- 1), cb_eval_fc_con_dyn)\n",
    "\n",
    "# Similarly to above, this callback does not evaluate the gradient of the objective.\n",
    "# As such, we pass `nV = 0`, and `objGradIndexVars = C_NULL` in the `KN_set_cb_grad()` call.\n",
    "KNITRO.KN_set_cb_grad(kc, cb_dyn, cb_eval_ga_con_dyn,\n",
    "                      nV=0, objGradIndexVars=C_NULL,\n",
    "                      jacIndexCons=jacIndexConsCB,\n",
    "                      jacIndexVars=jacIndexVarsCB)\n",
    "\n",
    "nnzH = length(hessIndexVars1CB)\n",
    "KNITRO.KN_set_cb_hess(kc, cb_dyn, nnzH, cb_eval_h_con_dyn,\n",
    "                      hessIndexVars1=hessIndexVars1CB,\n",
    "                      hessIndexVars2=hessIndexVars2CB)\n",
    "\n",
    "KNITRO.KN_set_cb_user_params(kc, cb_dyn, jac_data)\n",
    "\n",
    "# Add objective\n",
    "objIndex, objCoef = ind_θ[end], 1.0\n",
    "KNITRO.KN_add_obj_linear_struct(kc, objIndex - 1, objCoef)\n",
    "KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MAXIMIZE)\n",
    "\n",
    "KNITRO.KN_set_newpt_callback(kc, (kc, x, duals, userParams) -> (print(\"\"); return 0))\n",
    "\n",
    "KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_ALG,        KNITRO.KN_ALG_BAR_DIRECT)\n",
    "KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_BAR_MURULE, KNITRO.KN_BAR_MURULE_QUALITY)\n",
    "KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_LINSOLVER,  KNITRO.KN_LINSOLVER_MA57)\n",
    "KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_MAXTIMECPU, 20.0)  # Default: 1.0e8\n",
    "\n",
    "# KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_HESSOPT, KNITRO.KN_HESSOPT_LBFGS)\n",
    "# KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_LMSIZE, 5)  # limited-memory pairs stored\n",
    "\n",
    "# KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, KNITRO.KN_OUTLEV_ITER_VERBOSE)\n",
    "\n",
    "# # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
    "# # Perform a derivative check.\n",
    "# KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_DERIVCHECK, KNITRO.KN_DERIVCHECK_SECOND)\n",
    "# KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_DERIVCHECK_TERMINATE, KNITRO.KN_DERIVCHECK_STOPALWAYS)\n",
    "# KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_DERIVCHECK_TOL, 1.0e-5)  # Default value: 1.0e-6\n",
    "# KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_DERIVCHECK_TYPE, KNITRO.KN_DERIVCHECK_CENTRAL)\n",
    "\n",
    "KNITRO.KN_solve(kc)\n",
    "nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)\n",
    "cpu_time = KNITRO.KN_get_solve_time_real(kc)\n",
    "\n",
    "KNITRO.KN_free(kc)\n",
    "KNITRO.KN_release_license(lm)\n",
    "\n",
    "print(\"Final latitude θ = $(round(rad2deg(x[ind_θ[end]]), digits = 2))°\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = cumsum([0 ; x[ind_t][1:end-1]])\n",
    "\n",
    "plt_altitude = plot(ts, x[ind_h], legend = nothing, title = \"Altitude (100,000 ft)\")\n",
    "plt_longitude = plot(ts, rad2deg.(x[ind_ϕ]), legend = nothing, title = \"Longitude (deg)\")\n",
    "plt_latitude = plot(ts, rad2deg.(x[ind_θ]), legend = nothing, title = \"Latitude (deg)\")\n",
    "plt_velocity = plot(ts, x[ind_v], legend = nothing, title = \"Velocity (1000 ft/sec)\")\n",
    "plt_flight_path = plot(ts, rad2deg.(x[ind_γ]), legend = nothing, title = \"Flight Path (deg)\")\n",
    "plt_azimuth = plot(ts, rad2deg.(x[ind_ψ]), legend = nothing, title = \"Azimuth (deg)\")\n",
    "\n",
    "plt = plot(plt_altitude,  plt_velocity, plt_longitude, plt_flight_path, plt_latitude,\n",
    "           plt_azimuth, layout = grid(3, 2), linewidth = 2, size = (700, 700))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt_attack_angle = plot(ts[1:end-1], rad2deg.(x[ind_α][1:end-1]), legend = nothing, title = \"Angle of Attack (deg)\")\n",
    "plt_bank_angle = plot(ts[1:end-1], rad2deg.(x[ind_β][1:end-1]), legend = nothing, title = \"Bank Angle (deg)\")\n",
    "plt_heating = plot(ts, q.(x[ind_h]*1e5, x[ind_v]*1e4, x[ind_α]), legend = nothing, title = \"Heating (BTU/ft/ft/sec)\") \n",
    "\n",
    "plt = plot(plt_attack_angle, plt_bank_angle, plt_heating, layout = grid(3, 1), linewidth = 2, size = (700, 700))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt = plot(rad2deg.(x[ind_ϕ]), rad2deg.(x[ind_θ]), x[ind_h], linewidth = 2,\n",
    "           legend = nothing, title = \"Space Shuttle Reentry Trajectory\",\n",
    "           xlabel = \"Longitude (deg)\", ylabel = \"Latitude (deg)\", zlabel = \"Altitude (100,000 ft)\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.2",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}