#!/bin/bash
# Phase 3: External aerodynamics parametric yaw sweep
# Yaw angles: 0, 5, 10, 15, 20 degrees

set -e
source /usr/lib/openfoam/openfoam2406/etc/bashrc

BASEDIR="$(cd "$(dirname "$0")" && pwd)"
cd "$BASEDIR"

# Step 1: Generate mesh (only once)
echo "=== Running blockMesh ==="
blockMesh 2>&1 | tail -5

echo "=== Running snappyHexMesh ==="
snappyHexMesh -overwrite 2>&1 | tail -10

echo "=== Mesh complete ==="

# Check mesh
checkMesh 2>&1 | tail -20

# Save the base mesh and 0/ directory
cp -r constant/polyMesh constant/polyMesh_base
cp -r 0 0.base

# Yaw angles to sweep
ANGLES="0 5 10 15 20"

# Results file
RESULTS="$BASEDIR/results.csv"
echo "yaw_deg,Ux,Uy,Cd,Cs,Cl,Cmx,Cmy,Cmz" > "$RESULTS"

for ANGLE in $ANGLES; do
    echo ""
    echo "========================================="
    echo "=== Yaw angle: ${ANGLE} degrees ==="
    echo "========================================="

    # Compute velocity components
    # Yaw = rotation about Z axis
    # Ux = U * cos(yaw), Uy = U * sin(yaw)
    UX=$(python3 -c "import math; print(f'{10*math.cos(math.radians($ANGLE)):.6f}')")
    UY=$(python3 -c "import math; print(f'{10*math.sin(math.radians($ANGLE)):.6f}')")

    echo "Ux=$UX, Uy=$UY"

    # Clean previous run but keep mesh
    rm -rf [1-9]* processor*
    # Restore base 0/ directory
    rm -rf 0
    cp -r 0.base 0

    # Update 0/U with new velocity direction
    # Also update dragDir and liftDir (keep drag in freestream direction)
    # dragDir = (cos, sin, 0), sideForce perpendicular = (-sin, cos, 0)
    DRAG_X=$UX_NORM
    cat > 0/U << UEOF
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform ($UX $UY 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform ($UX $UY 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    sides
    {
        type            fixedValue;
        value           uniform ($UX $UY 0);
    }

    top
    {
        type            slip;
    }

    bottom
    {
        type            slip;
    }

    droneBody
    {
        type            noSlip;
    }
}
UEOF

    # Update controlDict with correct drag/side force directions for this yaw
    COSX=$(python3 -c "import math; print(f'{math.cos(math.radians($ANGLE)):.6f}')")
    SINX=$(python3 -c "import math; print(f'{math.sin(math.radians($ANGLE)):.6f}')")
    NSINX=$(python3 -c "import math; print(f'{-math.sin(math.radians($ANGLE)):.6f}')")

    cat > system/controlDict << CEOF
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      controlDict;
}

application     simpleFoam;

startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         500;
deltaT          1;

writeControl    timeStep;
writeInterval   500;
purgeWrite      1;

writeFormat     ascii;
writePrecision  8;
writeCompression off;
timeFormat      general;
timePrecision   6;

runTimeModifiable true;

functions
{
    forceCoeffs1
    {
        type            forceCoeffs;
        libs            (forces);
        writeControl    timeStep;
        writeInterval   1;
        log             true;

        patches         (droneBody);
        rho             rhoInf;
        rhoInf          1.225;
        liftDir         (0 0 1);
        dragDir         ($COSX $SINX 0);
        CofR            (0 0 0);
        pitchAxis       (0 1 0);
        magUInf         10;
        lRef            0.092;
        Aref            0.006648;
    }

    forces1
    {
        type            forces;
        libs            (forces);
        writeControl    timeStep;
        writeInterval   1;
        log             true;

        patches         (droneBody);
        rho             rhoInf;
        rhoInf          1.225;
        CofR            (0 0 0);
    }
}
CEOF

    echo "=== Running simpleFoam for yaw=$ANGLE ==="
    simpleFoam 2>&1 | tail -20

    # Extract final force coefficients
    # forceCoeffs output: Time Cd Cs Cl CmRoll CmPitch CmYaw Cd(f) Cd(r) Cs(f) Cs(r) Cl(f) Cl(r)
    COEFFS_FILE="postProcessing/forceCoeffs1/0/coefficient.dat"
    if [ ! -f "$COEFFS_FILE" ]; then
        # Try alternate name
        COEFFS_FILE=$(find postProcessing/forceCoeffs1 -name "*.dat" | head -1)
    fi

    if [ -f "$COEFFS_FILE" ]; then
        # Get last line (converged values)
        LAST=$(tail -1 "$COEFFS_FILE")
        CD=$(echo "$LAST" | awk '{print $2}')
        CS=$(echo "$LAST" | awk '{print $3}')
        CL=$(echo "$LAST" | awk '{print $4}')
        CMR=$(echo "$LAST" | awk '{print $5}')
        CMP=$(echo "$LAST" | awk '{print $6}')
        CMY=$(echo "$LAST" | awk '{print $7}')
        echo "Cd=$CD, Cs=$CS, Cl=$CL, CmYaw=$CMY"
        echo "$ANGLE,$UX,$UY,$CD,$CS,$CL,$CMR,$CMP,$CMY" >> "$RESULTS"
    else
        echo "WARNING: No coefficient file found for yaw=$ANGLE"
        echo "$ANGLE,$UX,$UY,NaN,NaN,NaN,NaN,NaN,NaN" >> "$RESULTS"
    fi

    # Save postProcessing for this angle
    if [ -d postProcessing ]; then
        cp -r postProcessing "postProcessing_yaw${ANGLE}"
    fi
done

echo ""
echo "=== Sweep complete ==="
echo "Results saved to $RESULTS"
cat "$RESULTS"