diff --git a/src/weather/ERA5/sample.grib b/src/weather/ERA5/sample.grib
new file mode 100644
index 00000000..d4d43216
Binary files /dev/null and b/src/weather/ERA5/sample.grib differ
diff --git a/src/weather/ERA5/sample.grib.5b7b6.idx b/src/weather/ERA5/sample.grib.5b7b6.idx
new file mode 100755
index 00000000..b7a495e5
Binary files /dev/null and b/src/weather/ERA5/sample.grib.5b7b6.idx differ
diff --git a/src/weather/draw.py b/src/weather/draw.py
new file mode 100644
index 00000000..f7e76563
--- /dev/null
+++ b/src/weather/draw.py
@@ -0,0 +1,49 @@
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+from pyproj import Geod
+import numpy as np
+
+
+def plot_flight_arc(mission):
+    
+    # Set up map
+    
+    fig = plt.figure(figsize=(10, 7))
+    ax = plt.axes(projection=ccrs.PlateCarree())
+    ax.set_extent([-130, -65, 22, 50], crs=ccrs.PlateCarree())
+    
+    # Add map features
+    ax.add_feature(cfeature.STATES, linewidth=0.5)
+    ax.add_feature(cfeature.COASTLINE)
+    ax.add_feature(cfeature.BORDERS, linestyle=':')
+    ax.add_feature(cfeature.LAND, facecolor='lightgray')
+    ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
+    ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5)
+    
+    # Extract OD lat lon position
+    lon_dep, lat_dep, _ = mission["dep_location"]
+    lon_arr, lat_arr, _ = mission["arr_location"]
+    
+    # Use WGS84 ellipse definition
+    geod = Geod(ellps="WGS84")
+    
+    # Get 100 equally spaced points along the arc
+    points = geod.npts(lon_dep, lat_dep, lon_arr, lat_arr, 100)
+    lons = [lon_dep] + [pt[0] for pt in points] + [lon_arr]
+    lats = [lat_dep] + [pt[1] for pt in points] + [lat_arr]
+
+    ax.plot(lons, lats, 'k-', linewidth=0.8, alpha=0.7)
+    
+    # Mark endpoints
+    ax.plot(lon_dep, lat_dep, 'go', markersize=6, label="Departure")
+    ax.plot(lon_arr, lat_arr, 'ro', markersize=6, label="Arrival")
+    
+    # Add legend and title
+    ax.legend(loc='lower left')
+    dep_code = mission['dep_airport']
+    arr_code = mission['arr_airport']
+    ax.set_title(f"Flight Arc: {dep_code} → {arr_code}")
+    
+    plt.tight_layout()
+    plt.savefig("sample.png", dpi=300)
\ No newline at end of file
diff --git a/src/weather/main.py b/src/weather/main.py
new file mode 100644
index 00000000..5f4fbb17
--- /dev/null
+++ b/src/weather/main.py
@@ -0,0 +1,92 @@
+import json
+import draw as dr
+import utils as util
+import trajectory as traj
+
+
+mission_path = "../missions/sample_missions_10.json"
+weather_data = "ERA5/sample.grib"
+
+# Load JSON data
+with open(mission_path, 'r') as file:
+    missions = json.load(file)
+    
+
+# Select a dummy mission for now
+first_mission = missions[0]
+
+
+# Get mission points based on mission def
+trajectory = traj.create_dummy_traj(first_mission)
+
+# Unpack trajectory -->
+lons = trajectory["lons"]
+lats = trajectory["lats"]
+alts = trajectory["H"]
+tas = trajectory["TAS"]
+
+# Initialize lists -->
+gs_list, heading_list, u_list, v_list = [], [], [], []
+
+# Build wind interpolator ->
+u_interp, v_interp, meta = util.build_era5_interpolators(weather_data)
+
+# Loop over all mission points to compute heading and ground speed ->
+for i in range(len(lons)):
+    if i < len(lons) - 1:
+        lon_next, lat_next = lons[i + 1], lats[i + 1]
+    else:
+        # For the last point, repeat previous heading
+        lon_next, lat_next = lons[i - 1], lats[i - 1]
+
+    gs, heading, u, v = util.compute_ground_speed(
+        lon=lons[i],
+        lat=lats[i],
+        lon_next=lon_next,
+        lat_next=lat_next,
+        alt_ft=alts[i],
+        tas_kts=tas[i],
+        u_interp=u_interp,
+        v_interp=v_interp,
+    )
+    
+    # Append to respective lists for sanity check ->
+    gs_list.append(gs)
+    heading_list.append(heading)
+    u_list.append(u)
+    v_list.append(v)
+    
+# Append to trajectory ->
+trajectory["GS"] = gs_list
+trajectory["heading_rad"] = heading_list
+trajectory["u_wind"] = u_list
+trajectory["v_wind"] = v_list
+
+
+# Print the first few points
+print("\n")
+print("#--------------------------------SAMPLE TRAJECTORY (SNAPSHOT) -----------------------------#\n")
+for lon, lat, tas, gs, alt in zip(trajectory["lons"][:5], trajectory["lats"][:5], trajectory["TAS"][:5], trajectory["GS"][:5], trajectory["H"][:5]):
+    print(f"Lon: {lon:.2f}, Lat: {lat:.2f}, TAS: {tas:.2f} kt, GS: {gs:.2f} kt, Alt: {alt:.2f} ft")
+print("#------------------------------------------------------------------------------------------#\n")
+
+
+
+#track, heading, drift, tas, u, v, wind_mag = util.get_tas(trajectory, "era5.grib")
+
+#for i in range(5):
+#    print(f"Pt {i}: Track={track[i]:.1f}°, Heading={heading[i]:.1f}°, Drift={drift[i]:+.1f}°, "
+#          f"TAS={tas[i]:.1f} kt, U={u[i]:.1f}, V={v[i]:.1f}, Wind={wind_mag[i]:.1f} m/s")
+
+
+
+#dr.plot_flight_arc(first_mission)
+
+#util.get_flight_track(first_mission)
+
+
+
+
+
+
+
diff --git a/src/weather/sample.png b/src/weather/sample.png
new file mode 100644
index 00000000..e4aa8acb
Binary files /dev/null and b/src/weather/sample.png differ
diff --git a/src/weather/trajectory.py b/src/weather/trajectory.py
new file mode 100644
index 00000000..d3d4865e
--- /dev/null
+++ b/src/weather/trajectory.py
@@ -0,0 +1,132 @@
+from pyproj import Geod
+import numpy as np
+
+
+
+def get_mission_points(mission):
+    """
+    Generates a discretized set of latitude and longitude points along a great-circle path
+    between departure and arrival locations, assuming constant cruise altitude and ground speed.
+
+    Parameters
+    ----------
+    mission : dict
+        Dictionary containing origin and destination coordinates with the following keys:
+            - 'dep_location': tuple of (longitude, latitude, altitude) for the departure point [degrees, degrees, feet]
+            - 'arr_location': tuple of (longitude, latitude, altitude) for the arrival point [degrees, degrees, feet]
+
+    Returns
+    -------
+    dict
+        A dictionary containing:
+            - 'lons' : list of longitudes along the flight path [degrees]
+            - 'lats' : list of latitudes along the flight path [degrees]
+            - 'GS'   : list of assumed constant ground speed at each point [knots]
+            - 'H'    : list of assumed constant cruise altitude at each point [feet]
+
+    Notes
+    -----
+    - The path is discretized using 100 intermediate points (plus endpoints), resulting in 102 total waypoints.
+    - This function uses `pyproj.Geod` with the WGS84 ellipsoid to compute a geodesic path.
+    - All values for speed and altitude are placeholders; replace them with mission-specific data in actual use.
+    """
+    
+    # Instantiate WGS84 ellipsoid
+    geod = Geod(ellps ="WGS84")
+    
+    # Extract OD lat-lon 
+    lon_dep, lat_dep, _ = mission["dep_location"]
+    lon_arr, lat_arr, _ = mission["arr_location"]
+    
+    # Assume 100 points for discretization (for demo only)
+    # Note: This will change when flying actual missions
+    
+    points = geod.npts(lon_dep, lat_dep, lon_arr, lat_arr, 100)
+    
+    lons = [lon_dep] + [pt[0] for pt in points] + [lon_arr]
+    lats = [lat_dep] + [pt[1] for pt in points] + [lat_arr]
+    
+    # Assign a dummy ground speed
+    ground_speeds = [450] * len(lons)
+    
+    # Assign a dummy cruise altitude
+    altitude_ft = [35000] * len(lons)
+    
+    return {
+        "lons": lons,
+        "lats": lats,
+        "GS": ground_speeds,
+        "H": altitude_ft
+    }
+    
+def create_dummy_traj(mission):
+    """
+    Generates a discretized trapezoidal trajectory profile (altitude and speed)
+    along a geodesic path between departure and arrival locations.
+
+    Parameters
+    ----------
+    mission : dict
+        Dictionary with:
+            - 'dep_location': (longitude, latitude, altitude) of departure [degrees, degrees, feet]
+            - 'arr_location': (longitude, latitude, altitude) of arrival [degrees, degrees, feet]
+
+    Returns
+    -------
+    dict
+        Dictionary containing:
+            - 'lons': list of longitudes [degrees]
+            - 'lats': list of latitudes [degrees]
+            - 'TAS'  : list of True Air Speed [knots]
+            - 'H'   : list of altitudes [feet]
+    """
+    # Instantiate GGS84 ellipsoid for geodesic calculation
+    geod = Geod(ellps="WGS84")
+    
+    # Extract dept + Arrival lat-lon
+    lon_dep, lat_dep, _ = mission["dep_location"]
+    lon_arr, lat_arr, _ = mission["arr_location"]
+    
+    # Assume 100 points for trajectory discretization (demo only)
+    n_total = 100
+    points = geod.npts(lon_dep, lat_dep, lon_arr, lat_arr, n_total)
+    lons = [lon_dep] + [pt[0] for pt in points] + [lon_arr]
+    lats = [lat_dep] + [pt[1] for pt in points] + [lat_arr]
+    
+    # Define climb and descent length as 25 % of total trajectory
+    n_climb = n_descent = n_total // 4
+    n_cruise = len(lons) - n_climb - n_descent
+    
+    # Define end points for speed and altitude 
+    alt_start, alt_cruise = 1000, 35000  # feet
+    spd_start, spd_cruise = 140, 450 # knots
+    
+    # Define climb profile
+    climb_alt = np.linspace(alt_start, alt_cruise, n_climb)
+    climb_spd = np.linspace(spd_start, spd_cruise, n_climb)
+    
+    # Define cruise profile
+    cruise_alt = np.full(n_cruise, alt_cruise)
+    cruise_spd = np.full(n_cruise, spd_cruise)
+    
+    # Define descent profile
+    descent_alt = np.linspace(alt_cruise, alt_start, n_descent)
+    descent_spd = np.linspace(spd_cruise, spd_start, n_descent)
+    
+    
+
+    # Collect height and TAS
+    H = np.concatenate([climb_alt, cruise_alt, descent_alt])
+    TAS = np.concatenate([climb_spd, cruise_spd, descent_spd])
+    
+    H = H[:len(lons)]
+    TAS = TAS[:len(lons)]
+    
+    return {
+        "lons": lons,
+        "lats": lats,
+        "TAS": TAS,
+        "H": H
+    }
+    
+    
\ No newline at end of file
diff --git a/src/weather/utils.py b/src/weather/utils.py
new file mode 100644
index 00000000..788324d1
--- /dev/null
+++ b/src/weather/utils.py
@@ -0,0 +1,331 @@
+from pyproj import Geod
+import numpy as np
+from scipy.interpolate import RegularGridInterpolator
+import xarray as xr
+
+def altitude_to_pressure_hpa(alt_ft):
+    """Convert altitude in feet to pressure in hPa using ISA approximation."""
+    alt_m = np.array(alt_ft) * 0.3048
+    pressure = 1013.25 * (1 - (0.0065 * alt_m) / 288.15) ** 5.2561
+    return pressure
+    
+
+def get_wind_at_points(mission_data, era5_path):
+    """
+    Interpolates U and V wind components from ERA5 at each lat/lon/alt point.
+
+    Parameters
+    ----------
+    mission_data : dict
+        Output from get_mission_points(); must contain 'lons', 'lats', and 'H'.
+    era5_path : str
+        Path to ERA5 GRIB file.
+
+    Returns
+    -------
+    u_vals : np.ndarray
+        Zonal wind (eastward) component at each point [m/s]
+    v_vals : np.ndarray
+        Meridional wind (northward) component at each point [m/s]
+    wind_speed : np.ndarray
+        Wind speed magnitude [m/s]
+    """
+    
+    # Load ERA4 GRIB file (single time slice)
+    ds = xr.open_dataset(era5_path, engine="cfgrib")
+    
+    # Ensure proper dimension order
+    lats_era = ds.latitude.values
+    lons_era = ds.longitude.values
+    levels_era = ds.isobaricInhPa.values
+    
+    # Sort lattitude in ascending order
+    if lats_era[0] > lats_era[-1]:
+        lats_era = lats_era[::-1]
+        ds['u'] = ds['u'][::-1, :, :]
+        ds['v'] = ds['v'][::-1, :, :]
+
+    # Prepare 3D interpolators
+    u_interp = RegularGridInterpolator(
+        (levels_era, lats_era, lons_era),
+        ds['u'].values,
+        bounds_error=False,
+        fill_value=np.nan
+    )
+    
+    v_interp = RegularGridInterpolator(
+        (levels_era, lats_era, lons_era),
+        ds['v'].values,
+        bounds_error=False,
+        fill_value=np.nan
+    )
+    
+    # Convert altitude (ft) to pressure (hPa)
+    pressures = altitude_to_pressure_hpa(mission_data['H'])
+    
+    # Convert lon to [0, 360] for ERA5 grid compatibility
+    lons_360 = [(lon + 360) if lon < 0 else lon for lon in mission_data['lons']]
+    
+    # Form (pressure, lat, lon) points for interpolation
+    points = np.array([
+        (p, lat, lon) for p, lat, lon in zip(pressures, mission_data['lats'], lons_360)
+    ])
+    
+    # Interpolate
+    u_vals = u_interp(points)
+    v_vals = v_interp(points)
+    
+    wind_speed = np.sqrt(u_vals**2 + v_vals**2)
+    
+    return u_vals, v_vals, wind_speed
+
+def get_tas(mission_data, era5_path):
+    """
+    Computes the true airspeed (TAS), heading, and drift angle along a flight path using geodesic 
+    track angles and interpolated wind components from ERA5 reanalysis data.
+
+    Parameters
+    ----------
+    mission_data : dict
+        Dictionary containing flight path data with the following keys:
+            - 'lats': array-like, latitudes of flight path points [degrees]
+            - 'lons': array-like, longitudes of flight path points [degrees]
+            - 'GS': array-like, ground speed at each point [knots]
+
+    era5_path : str
+        File path to the ERA5 dataset used for wind interpolation.
+
+    Returns
+    -------
+    track_angles : np.ndarray
+        Array of geodesic track angles (course) between consecutive lat-lon points [degrees].
+
+    headings : np.ndarray
+        Array of aircraft headings (direction of motion through the air) [degrees].
+
+    drifts : np.ndarray
+        Array of drift angles (heading - track), indicating crosswind correction required [degrees].
+
+    tas_vals : np.ndarray
+        Array of computed true airspeeds (TAS) along the path [knots].
+
+    u_vals : np.ndarray
+        Interpolated eastward wind components (u-wind) at each flight point [m/s].
+
+    v_vals : np.ndarray
+        Interpolated northward wind components (v-wind) at each flight point [m/s].
+
+    wind_speed : np.ndarray
+        Magnitude of interpolated wind speed at each point [m/s].
+
+    Notes
+    -----
+    - TAS is computed using the wind triangle from the vector difference between ground speed vector and wind vector.
+    - Drift angle is signed and bounded to the range [-180, 180] degrees.
+    - Requires `pyproj.Geod` and NumPy.
+    """
+
+    # Get wind components based on path
+    u_vals, v_vals, wind_speed = get_wind_at_points(mission_data, era5_path)
+    
+    # Compute geodesic track angles between consecutive flight segments
+    geod = Geod(ellps="WGS84")
+    lons = mission_data["lons"]
+    lats = mission_data["lats"]
+    gs = mission_data["GS"]
+    
+    track_angles = []
+    headings = []
+    drifts = []
+    tas_vals = []
+    
+    # Loop over all points along the arc
+    for i in range(len(lons) - 1):
+        # Compute track angle for the two lat-lon pairs
+        lon1, lat1 = lons[i], lats[i]
+        lon2, lat2 = lons[i + 1], lats[i + 1]
+        az_fwd, _, _ = geod.inv(lon1, lat1, lon2, lat2)
+        track = az_fwd % 360
+        track_angles.append(track)
+        
+    # Use wind triangle to compute heading and drift
+        track_rad = np.deg2rad(track)
+        vgx = gs[i] * np.cos(track_rad)  # ground speed-x along track
+        vgy = gs[i] * np.sin(track_rad)  # ground speed-y along track
+        
+        
+        # Substract wind vector (souce ERA-5 dummy weather)
+        vax = vgx - u_vals[i]   # TAS x-component
+        vay = vgy - v_vals[i]   # TAS y-component
+        
+        tas_mps = np.sqrt(vax**2 + vay**2)
+        tas_knots = tas_mps / 0.514444  # m/s -> knots
+        
+        tas_vals.append(tas_knots)
+        
+        heading_rad = np.arctan2(vay, vax)
+        heading_deg = np.rad2deg(heading_rad) % 360
+        
+        # Compute path drift
+        drift = (heading_deg - track + 360) % 360
+        if drift > 180:
+            drift -= 360  # Convert to signed drift
+        headings.append(heading_deg)
+        drifts.append(drift)
+        
+    return (
+        np.array(track_angles),
+        np.array(headings),
+        np.array(drifts),
+        np.array(tas_vals),
+        u_vals,
+        v_vals,
+        wind_speed
+    )
+    
+
+def apply_drift_to_mission_points(mission_data, drift_angles):
+    """
+    Computes drifted lat/lon positions by adjusting the heading from the track angle.
+
+    Parameters
+    ----------
+    mission_data : dict
+        Output from get_mission_points(), containing 'lons', 'lats'.
+    drift_angles : list or np.ndarray
+        Drift angles in degrees for each segment.
+
+    Returns
+    -------
+    drifted_lons : list
+        Adjusted longitude values after drift.
+    drifted_lats : list
+        Adjusted latitude values after drift.
+    """
+    geod = Geod(ellps="WGS84")
+    lons = mission_data["lons"]
+    lats = mission_data["lats"]
+
+    drifted_lons = [lons[0]]
+    drifted_lats = [lats[0]]
+
+    for i in range(len(drift_angles)):
+        lon1 = drifted_lons[-1]
+        lat1 = drifted_lats[-1]
+        lon2 = lons[i + 1]
+        lat2 = lats[i + 1]
+
+        # Compute track azimuth and distance from original points
+        fwd_azimuth, _, distance_m = geod.inv(lon1, lat1, lon2, lat2)
+
+        # Apply drift to adjust heading
+        drifted_heading = (fwd_azimuth + drift_angles[i]) % 360
+
+        # Project new position using drifted heading and same distance
+        lon_drifted, lat_drifted, _ = geod.fwd(lon1, lat1, drifted_heading, distance_m)
+
+        drifted_lons.append(lon_drifted)
+        drifted_lats.append(lat_drifted)
+
+    return drifted_lons, drifted_lats
+
+
+def build_era5_interpolators(era5_path):
+    """
+    Load GRIB file and build interpolators for u/v wind components.
+
+    Parameters
+    ----------
+    era5_path : str
+        Path to ERA5 GRIB file.
+
+    Returns
+    -------
+    tuple: (u_interp, v_interp, metadata)
+        Interpolators and grid metadata (lats, lons, levels).
+    """
+    ds = xr.open_dataset(era5_path, engine="cfgrib")
+
+    # ERA5 longitudes are in [0, 360]; convert to [-180, 180]
+    if (ds.longitude > 180).any():
+        ds = ds.assign_coords(longitude=((ds.longitude + 180) % 360) - 180)
+        ds = ds.sortby("longitude")
+
+    # Ensure latitudes are increasing
+    if np.any(np.diff(ds.latitude.values) < 0):
+        ds = ds.sortby("latitude")
+
+    # Extract dimensions
+    lats = ds.latitude.values
+    lons = ds.longitude.values
+    levels = ds.isobaricInhPa.values  # Pressure levels
+
+    # Extract u and v wind components (assumes 3D: level x lat x lon)
+    u_data = ds["u"].values
+    v_data = ds["v"].values
+
+    # Interpolator over (level, lat, lon) (Filling with NaNs for safety)
+    u_interp = RegularGridInterpolator((levels, lats, lons), u_data, bounds_error=False, fill_value=4.0)
+    v_interp = RegularGridInterpolator((levels, lats, lons), v_data, bounds_error=False, fill_value=4.0)
+
+    return u_interp, v_interp, {"levels": levels, "lats": lats, "lons": lons}
+
+def compute_ground_speed(lon, lat, lon_next, lat_next, alt_ft, tas_kts, u_interp, v_interp):
+    
+    """
+    Computes ground speed for a single point using TAS, heading, and interpolated winds.
+
+    Parameters
+    ----------
+    lon, lat : float
+        Longitude and latitude of the point [degrees]
+    lon_next, lat_next : float
+        Next point [deg] to compute heading
+    alt_ft : float
+        Altitude in feet
+    tas_kts : float
+        True airspeed in knots
+    heading_rad : float
+        Aircraft heading (radians from true north)
+    u_interp, v_interp : interpolator functions
+        ERA-5 wind interpolators (from build_era5_interpolators)
+
+    Returns
+    -------
+    gs_kts : float
+        Ground speed at this point [knots]
+    u, v : float
+        Interpolated wind components [m/s]
+    """
+    
+    # Convert altitude and TAS to S.I units
+    pressure_level = float(altitude_to_pressure_hpa(alt_ft))  # hPa
+    tas_ms = tas_kts * 0.514444                       # m/s
+    
+    
+    # Compute heading in radians from current to next point
+    geod = Geod(ellps="WGS84")
+    azimuth_deg, _, _ = geod.inv(lon, lat, lon_next, lat_next)
+    heading_rad = np.deg2rad(azimuth_deg)
+    
+    # Interpolate wind at this location
+    u = u_interp([[pressure_level, lat, lon]])[0]
+    v = v_interp([[pressure_level, lat, lon]])[0]
+    
+    
+    
+    # Airspeed vector in Earth frame
+    u_air = tas_ms * np.cos(heading_rad)
+    v_air = tas_ms * np.sin(heading_rad)
+    
+
+    
+    # Ground speed = air vector + wind vector
+    u_ground = u_air + u
+    v_ground = v_air + v
+    
+    gs_ms = np.sqrt(u_ground**2 + v_ground**2)
+    gs_kts = gs_ms / 0.514444
+    
+    return gs_kts, heading_rad, u, v
+                 
\ No newline at end of file