# V & V Script for NASA Inflow Measurement Reports


```python
"""V & V Script for NASA Inflow Measurement Reports"""
import csdl_alpha as csdl
import numpy as np
from BladeAD.utils.define_rotor_analysis import run_rotor_analysis
from BladeAD.core.airfoil.xfoil.two_d_airfoil_model import TwoDMLAirfoilModel
from BladeAD import NASA_INFLOW_DATA_FOLDER
import pandas as pd
from BladeAD.utils.verify_dynamic_inflow_models_plotting import plot_vnv_data


# start CSDL recorder
recorder = csdl.Recorder(inline=True)
recorder.start()

# Discretization & rotor geometry definition
num_nodes = 1
num_radial = 50
num_azimuthal = 100
num_blades = 4

radius = csdl.Variable(name="radius", value=0.860552)
chord_profile = csdl.Variable(
    name="chord_profile", 
    shape=(num_radial, ), 
    value=0.06604
)

# Note that the twist profile is linearly varying by 8 degrees
# and is such that the twist is zero at 3/4 of the blade radius
root_twist = np.deg2rad(8 - 2.470588235)
tip_twist = root_twist - np.deg2rad(8)
twist_profile = csdl.Variable(
    name="twist_profile", 
    shape=(num_radial, ), 
    value=np.linspace(root_twist, tip_twist, num_radial)
)

# The rotor blade uses a NACA 0012 airfoil
airfoil_model = TwoDMLAirfoilModel(
    airfoil_name="naca_0012",
)

# Thrust origin is the same for all test cases
thrust_origin=csdl.Variable(value=np.array([0. ,0., 0.]))

# --------------------NASA test case 1--------------------
alpha = -3.00 # rotor disk tilt angle, i.e., angle of attack
tilt_angle = np.deg2rad(90 + alpha)
tx = np.cos(tilt_angle)
tz = np.sin(tilt_angle)
thrust_vector=csdl.Variable(value=np.array([tx, 0, -tz]))
mesh_velocity_np = np.zeros((num_nodes, 3))
mesh_velocity_np[:, 0] = 28.50 # Free stream velocity
mesh_velocity = csdl.Variable(name="mesh_velocity", shape=(num_nodes, 3), value=mesh_velocity_np)
rpm = csdl.Variable(name="rpm", shape=(num_nodes,), value=2113)

# Read in the NASA data
df_1 = pd.read_csv(f"{NASA_INFLOW_DATA_FOLDER}/nasa_report_rotor_inflow_data_mu_015.csv")

# ----- Run Pitt-Peters model -----
pitt_peters_ouputs = run_rotor_analysis(
    model="pitt_peters",
    radius=radius,
    chord_profile=chord_profile,
    twist_profile=twist_profile,
    rpm=rpm,
    mesh_velocity=mesh_velocity,
    thrust_vector=thrust_vector,
    thrust_origin=thrust_origin,
    num_nodes=num_nodes,
    num_radial=num_radial,
    num_azimuthal=num_azimuthal,
    num_blades=num_blades,
    theta_0=np.deg2rad(9.37),
    theta_1_c=np.deg2rad(1.11),
    theta_1_s=np.deg2rad(-3.23),
    xi_0=np.deg2rad(-0.95),
    airfoil_model=airfoil_model,
)

inflow_pp = pitt_peters_ouputs.inflow
C_T_pp = pitt_peters_ouputs.thrust_coefficient.value

# Plot the inflow data for pitt-peters model
plot_vnv_data(
    nasa_data_df=df_1,
    nasa_C_T=0.0064,
    blade_ad_data=[inflow_pp.value],
    num_radial=num_radial,
    num_azimuthal=num_azimuthal,
    model_name=["Pitt--Peters model"],
    model_C_T=[C_T_pp[0]],
)

# ----- Run Peters-He model for different number of flow states-----
def run_peters_he(Q, M):
    return run_rotor_analysis(
        model="peters_he",
        radius=radius,
        chord_profile=chord_profile,
        twist_profile=twist_profile,
        rpm=rpm,
        mesh_velocity=mesh_velocity,
        thrust_vector=thrust_vector,
        thrust_origin=thrust_origin,
        num_nodes=num_nodes,
        num_radial=num_radial,
        num_azimuthal=num_azimuthal,
        num_blades=num_blades,
        theta_0=np.deg2rad(9.37),
        theta_1_c=np.deg2rad(1.11),
        theta_1_s=np.deg2rad(-3.23),
        xi_0=np.deg2rad(-0.95),
        airfoil_model=airfoil_model,
        Q=Q,
        M=M,
    )

peters_he_outputs = [run_peters_he(Q, M) for Q, M in [(2, 1), (4, 3), (6, 5)]]

# Extract inflow and thrust coefficient values
inflow_ph = [output.inflow.value for output in peters_he_outputs]
C_T_ph = [output.thrust_coefficient.value[0] for output in peters_he_outputs]

# Plot Peters-He model results
plot_vnv_data(
    nasa_data_df=df_1,
    nasa_C_T=0.0064,
    blade_ad_data=inflow_ph,
    num_radial=num_radial,
    num_azimuthal=num_azimuthal,
    model_name=[
        "Peters--He w/ 6 flow states",
        "Peters--He w/ 15 flow states",
        "Peters--He w/ 28 flow states",
    ],
    model_C_T=C_T_ph,
)

# --------------------NASA test case 2--------------------
alpha = -3.04 # rotor disk tilt angle, i.e., angle of attack
tilt_angle = np.deg2rad(90 + alpha)
tx = np.cos(tilt_angle)
tz = np.sin(tilt_angle)
thrust_vector=csdl.Variable(value=np.array([tx, 0, -tz]))
mesh_velocity_np = np.zeros((num_nodes, 3))
mesh_velocity_np[:, 0] = 43.8605044 # Free stream velocity
mesh_velocity = csdl.Variable(name="mesh_velocity", shape=(num_nodes, 3), value=mesh_velocity_np)
rpm = csdl.Variable(name="rpm", shape=(num_nodes,), value=2113)

# Read in the NASA data
df_2 = pd.read_csv(f"{NASA_INFLOW_DATA_FOLDER}/nasa_report_rotor_inflow_data_mu_023.csv")

# ----- Run Pitt-Peters model -----
pitt_peters_ouputs = run_rotor_analysis(
    model="pitt_peters",
    radius=radius,
    chord_profile=chord_profile,
    twist_profile=twist_profile,
    rpm=rpm,
    mesh_velocity=mesh_velocity,
    thrust_vector=thrust_vector,
    thrust_origin=thrust_origin,
    num_nodes=num_nodes,
    num_radial=num_radial,
    num_azimuthal=num_azimuthal,
    num_blades=num_blades,
    theta_0=np.deg2rad(8.16),
    theta_1_c=np.deg2rad(1.52),
    theta_1_s=np.deg2rad(-4.13),
    xi_0=np.deg2rad(-0.9),
    airfoil_model=airfoil_model,
)

inflow_pp = pitt_peters_ouputs.inflow
C_T_pp = pitt_peters_ouputs.thrust_coefficient.value

# Plot the inflow data for pitt-peters model
plot_vnv_data(
    nasa_data_df=df_2,
    nasa_C_T=0.0064,
    blade_ad_data=[inflow_pp.value],
    num_radial=num_radial,
    num_azimuthal=num_azimuthal,
    model_name=["Pitt--Peters model"],
    model_C_T=[C_T_pp[0]],
)

# ----- Run Peters-He model for different number of flow states-----
def run_peters_he(Q, M):
    return run_rotor_analysis(
        model="peters_he",
        radius=radius,
        chord_profile=chord_profile,
        twist_profile=twist_profile,
        rpm=rpm,
        mesh_velocity=mesh_velocity,
        thrust_vector=thrust_vector,
        thrust_origin=thrust_origin,
        num_nodes=num_nodes,
        num_radial=num_radial,
        num_azimuthal=num_azimuthal,
        num_blades=num_blades,
        theta_0=np.deg2rad(8.16),
        theta_1_c=np.deg2rad(1.52),
        theta_1_s=np.deg2rad(-4.13),
        xi_0=np.deg2rad(-0.9),
        airfoil_model=airfoil_model,
        Q=Q,
        M=M,
    )

peters_he_outputs = [run_peters_he(Q, M) for Q, M in [(2, 2), (4, 4), (6, 6)]]

# Extract inflow and thrust coefficient values
inflow_ph = [output.inflow.value for output in peters_he_outputs]
C_T_ph = [output.thrust_coefficient.value[0] for output in peters_he_outputs]

# Plot Peters-He model results
plot_vnv_data(
    nasa_data_df=df_2,
    nasa_C_T=0.0064,
    blade_ad_data=inflow_ph,
    num_radial=num_radial,
    num_azimuthal=num_azimuthal,
    model_name=[
        "Peters--He w/ 6 flow states",
        "Peters--He w/ 15 flow states",
        "Peters--He w/ 28 flow states",
    ],
    model_C_T=C_T_ph,
)


# --------------------NASA test case 3--------------------
alpha = -5.7 # rotor disk tilt angle, i.e., angle of attack
tilt_angle = np.deg2rad(90 + alpha)
tx = np.cos(tilt_angle)
tz = np.sin(tilt_angle)
thrust_vector=csdl.Variable(value=np.array([tx, 0, -tz]))
mesh_velocity_np = np.zeros((num_nodes, 3))
mesh_velocity_np[:, 0] = 66.75 # Free stream velocity
mesh_velocity = csdl.Variable(name="mesh_velocity", shape=(num_nodes, 3), value=mesh_velocity_np)
rpm = csdl.Variable(name="rpm", shape=(num_nodes,), value=2113)

# Read in the NASA data
df_3 = pd.read_csv(f"{NASA_INFLOW_DATA_FOLDER}/nasa_report_rotor_inflow_data_mu_35.csv")

# ----- Run Pitt-Peters model -----
pitt_peters_ouputs = run_rotor_analysis(
    model="pitt_peters",
    radius=radius,
    chord_profile=chord_profile,
    twist_profile=twist_profile,
    rpm=rpm,
    mesh_velocity=mesh_velocity,
    thrust_vector=thrust_vector,
    thrust_origin=thrust_origin,
    num_nodes=num_nodes,
    num_radial=num_radial,
    num_azimuthal=num_azimuthal,
    num_blades=num_blades,
    theta_0=np.deg2rad(9.2),
    theta_1_c=np.deg2rad(0.3),
    theta_1_s=np.deg2rad(-6.80),
    xi_0=np.deg2rad(-1.3),
    airfoil_model=airfoil_model,
)

inflow_pp = pitt_peters_ouputs.inflow
C_T_pp = pitt_peters_ouputs.thrust_coefficient.value

# Plot the inflow data for pitt-peters model
plot_vnv_data(
    nasa_data_df=df_3,
    nasa_C_T=0.0064,
    blade_ad_data=[inflow_pp.value],
    num_radial=num_radial,
    num_azimuthal=num_azimuthal,
    model_name=["Pitt--Peters model"],
    model_C_T=[C_T_pp[0]],
)

# ----- Run Peters-He model for different number of flow states-----
def run_peters_he(Q, M):
    return run_rotor_analysis(
        model="peters_he",
        radius=radius,
        chord_profile=chord_profile,
        twist_profile=twist_profile,
        rpm=rpm,
        mesh_velocity=mesh_velocity,
        thrust_vector=thrust_vector,
        thrust_origin=thrust_origin,
        num_nodes=num_nodes,
        num_radial=num_radial,
        num_azimuthal=num_azimuthal,
        num_blades=num_blades,
        theta_0=np.deg2rad(9.2),
        theta_1_c=np.deg2rad(0.3),
        theta_1_s=np.deg2rad(-6.80),
        xi_0=np.deg2rad(-1.3),
        airfoil_model=airfoil_model,
        Q=Q,
        M=M,
    )

peters_he_outputs = [run_peters_he(Q, M) for Q, M in [(2, 2), (4, 4), (6, 6)]]

# Extract inflow and thrust coefficient values
inflow_ph = [output.inflow.value for output in peters_he_outputs]
C_T_ph = [output.thrust_coefficient.value[0] for output in peters_he_outputs]

# Plot Peters-He model results
plot_vnv_data(
    nasa_data_df=df_3,
    nasa_C_T=0.0064,
    blade_ad_data=inflow_ph,
    num_radial=num_radial,
    num_azimuthal=num_azimuthal,
    model_name=[
        "Peters--He w/ 6 flow states",
        "Peters--He w/ 15 flow states",
        "Peters--He w/ 28 flow states",
    ],
    model_C_T=C_T_ph,
)
```