diff --git a/CHANGELOG.md b/CHANGELOG.md index da13043..24b3056 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,3 +15,5 @@ This file records dated repository-level documentation and implementation change - `2026-03-28`: Documented professor guidance that the atmosphere model must be integrated into the simulation, must include the stratosphere branch for altitude `> 11 km`, and will not run until the `XX` placeholders are replaced. - `2026-03-28`: Clarified file roles so `MEMORY.md` is the durable handoff and `TODO.md` is the active execution tracker. - `2026-03-28`: Confirmed the stratosphere equations from `docs/Addendum to Atomosphere Model Feb 26.pdf` and recorded the required `z > 11000 m` branch for the future `getRho.m` implementation. +- `2026-03-28`: Implemented [`code/getRho.m`](./code/getRho.m), filled the known scalar assignment constants in [`code/Interceptor_3DOF.m`](./code/Interceptor_3DOF.m) and [`code/getGuidance.m`](./code/getGuidance.m), assumed `Maz = 0 deg` for the current coplanar setup, validated `getRho.m` in Octave, and confirmed that the next blocker is the unresolved target/missile data arrays in [`code/T3dxdt.m`](./code/T3dxdt.m) and [`code/M3dxdt.m`](./code/M3dxdt.m). +- `2026-03-28`: Transcribed the target and missile thrust, mass, area, and drag data into [`code/T3dxdt.m`](./code/T3dxdt.m) and [`code/M3dxdt.m`](./code/M3dxdt.m), fixed row/column handling in [`code/getGuidance.m`](./code/getGuidance.m), verified the dynamics helpers in Octave, and found that the remaining issue is slow end-to-end runtime in the nested `ode45` loop. diff --git a/MEMORY.md b/MEMORY.md index 805720d..40cdfdd 100644 --- a/MEMORY.md +++ b/MEMORY.md @@ -41,6 +41,9 @@ This file is the durable handoff for AI agents working in this repository. Keep - Compute `P11k = P0*(T/T0)^(-g/(lapserate*R))`. - Compute `P = P11k*exp(-g*(z - 11000)/(R*T))`. - Then compute `rho = P/(R*T)` and `acousticSpeed = sqrt(gamma*R*T)`. +- Active implementation assumptions: + - `Maz = 0 deg` is currently assumed in `code/Interceptor_3DOF.m` to preserve the coplanar `y = 0` engagement geometry until a source document says otherwise. + - The target and missile thrust, mass, and drag arrays were manually transcribed from the image-based plots in `docs/Missile and Target Data for 3DOF Spring 26.pdf` on `2026-03-28`. - Assignment baseline confirmed from the handout: - Velocity pursuit guidance gain = `0.5` - Blind range = `2 m` @@ -60,18 +63,20 @@ This file is the durable handoff for AI agents working in this repository. Keep ## Code Status - [`code/Interceptor_3DOF.m`](./code/Interceptor_3DOF.m) drives the integration loop, plotting, and CSV export. -- [`code/T3dxdt.m`](./code/T3dxdt.m), [`code/M3dxdt.m`](./code/M3dxdt.m), and [`code/getGuidance.m`](./code/getGuidance.m) still contain `XX` placeholders. -- `getRho.m` is referenced by the driver and both dynamics functions, but the file is missing. -- Because `getRho.m` does not exist yet, neither the lower-atmosphere branch nor the addendum-defined stratosphere branch for `z > 11 km` is implemented in the repo. +- [`code/getRho.m`](./code/getRho.m) now exists and implements both the lower-atmosphere model and the addendum-defined stratosphere branch for `z > 11 km`. +- [`code/Interceptor_3DOF.m`](./code/Interceptor_3DOF.m) and [`code/getGuidance.m`](./code/getGuidance.m) now contain the known scalar assignment constants. +- [`code/T3dxdt.m`](./code/T3dxdt.m) and [`code/M3dxdt.m`](./code/M3dxdt.m) now contain manually transcribed target/missile reference areas, thrust curves, mass curves, and drag tables. +- [`code/getGuidance.m`](./code/getGuidance.m) now normalizes input shapes to column vectors so `M3dxdt` can call it without row/column mismatch failures. +- Octave now evaluates `getRho.m`, `T3dxdt.m`, and `M3dxdt.m` successfully, and the full script begins integrating without placeholder errors. +- End-to-end validation is still incomplete because the current nested `ode45` loop runs very slowly in Octave and was not observed through final plot/CSV generation within the available runtime. - [`code/getMD.m`](./code/getMD.m) exists and computes miss distance at the point of closest approach. -- The missile/target data PDF appears image-based; a future agent should log the extraction method used when those tables are transcribed. ## Next Actions -1. Implement `getRho.m` and integrate it into the simulation, including the required stratosphere branch for altitude `> 11 km`. -2. Transcribe missile and target geometry, thrust, mass, and drag data from the missile/target packet. -3. Replace all remaining `XX` placeholders with sourced values. -4. Run the MATLAB model and verify the required plots and CSV outputs. +1. Complete an end-to-end Octave validation run and confirm that the required plots and CSV files are produced. +2. Investigate the runtime cost of the nested `ode45` loop in [`code/Interceptor_3DOF.m`](./code/Interceptor_3DOF.m) if full validation remains too slow. +3. Confirm whether the temporary `Maz = 0 deg` assumption should remain. +4. Decide whether the manually transcribed target/missile curves need a second pass for tighter fidelity. ## Update Protocol For Agents diff --git a/README.md b/README.md index 1c5e36a..251aa14 100644 --- a/README.md +++ b/README.md @@ -42,16 +42,16 @@ The required outputs are: ## Current Code State - [`code/Interceptor_3DOF.m`](./code/Interceptor_3DOF.m) is the top-level driver script and writes `Ttrajectory.csv` and `Mtrajectory.csv`. -- [`code/T3dxdt.m`](./code/T3dxdt.m) and [`code/M3dxdt.m`](./code/M3dxdt.m) contain the target and interceptor dynamics models. +- [`code/T3dxdt.m`](./code/T3dxdt.m) and [`code/M3dxdt.m`](./code/M3dxdt.m) now include transcribed target and interceptor area, thrust, mass, and drag data from the course plots. - [`code/getGuidance.m`](./code/getGuidance.m) holds the guidance-law logic. +- [`code/getRho.m`](./code/getRho.m) now implements the lower-atmosphere model and the addendum-defined stratosphere branch. - [`code/getMD.m`](./code/getMD.m) computes miss distance near the point of closest approach. At the moment, the repository is a partially completed assignment skeleton rather than a finished simulation: -- Multiple required values are still marked as `XX`. -- `getRho.m` is referenced by the code but is not present in the repository, so the atmosphere model has not been integrated yet. -- When `getRho.m` is added, it must support both the lower-atmosphere branch and the addendum-defined stratosphere branch for `z > 11 km`. -- The missile/target data packet still needs to be fully transcribed into MATLAB arrays and constants. +- The executable `XX` placeholders have been removed from the current MATLAB code. +- The target and missile curves were manually transcribed from the image-based course plots, so those values should be treated as sourced approximations rather than machine-extracted ground truth. +- The model starts running in Octave, but the nested fixed-step `ode45` loop is slow enough that end-to-end validation is still incomplete. ## Working Expectations diff --git a/TODO.md b/TODO.md index 531cc6d..b35876f 100644 --- a/TODO.md +++ b/TODO.md @@ -6,19 +6,17 @@ Last updated: `2026-03-28` ## Now -- [ ] Implement `getRho.m` and include the required stratosphere branch for altitude `> 11 km`. -- [ ] Fill the known assignment constants in [`code/Interceptor_3DOF.m`](./code/Interceptor_3DOF.m) and [`code/getGuidance.m`](./code/getGuidance.m). -- [ ] Replace the current syntax-breaking `XX` placeholders that prevent the model from starting in Octave. -- [ ] Use [`docs/Addendum to Atomosphere Model Feb 26.pdf`](./docs/Addendum%20to%20Atomosphere%20Model%20Feb%2026.pdf) when implementing the `z > 11000 m` atmosphere branch. +- [ ] Complete an end-to-end Octave validation run and confirm that the required plots and CSV files are produced. +- [ ] Investigate why the nested `ode45` loop is slow in Octave if the full run remains impractical. +- [ ] Confirm whether the temporary `Maz = 0 deg` assumption should remain. ## Blocked -- [ ] Transcribe missile and target reference areas, thrust curves, mass curves, and drag tables from [`docs/Missile and Target Data for 3DOF Spring 26.pdf`](./docs/Missile%20and%20Target%20Data%20for%203DOF%20Spring%2026.pdf). +- [ ] Full validation is currently limited by Octave runtime, not by missing source data. ## Next -- [ ] Replace the remaining `XX` arrays in [`code/T3dxdt.m`](./code/T3dxdt.m) and [`code/M3dxdt.m`](./code/M3dxdt.m) with sourced data. -- [ ] Run [`code/Interceptor_3DOF.m`](./code/Interceptor_3DOF.m) in Octave and fix the next runtime errors. +- [ ] Decide whether the manually transcribed thrust, mass, and drag curves need a second pass for tighter fidelity. - [ ] Verify the required plots and inspect `Ttrajectory.csv` and `Mtrajectory.csv`. ## Done Criteria diff --git a/code/Interceptor_3DOF.m b/code/Interceptor_3DOF.m index 0875251..2b943c8 100644 --- a/code/Interceptor_3DOF.m +++ b/code/Interceptor_3DOF.m @@ -1,8 +1,8 @@ %% * Interceptor 3DOF* %% Executive clear all; clc; -t0 = 0; % Initial time. Target launch time. -tCommit = XX.; % Commit time (sec) +t0 = 0; % Initial time. Target launch time. +tCommit = 50.; % Commit time (sec) sigEl = 0; % Missile elevation angle standard deviation maxTime = 150; % estimated maximum engagement time in sec dt = 0.01; % integration interval (sec) @@ -21,9 +21,10 @@ Toutput_vector(1,:)= [t0,targetState',TxDot(4:6)']; %% Missile initial conditions Mp0 = [0; 0; 0]; % Missile initial position in ENU frame (m) -Mspeed0 = 1.e-06;% Missile initial speed in ENU frame (m/sec); must be non-zero for unit vector calculation. -Maz = XX; Mel = XX; % Missile launch azimuth and elevation (deg). -Mv0 = Mspeed0*[cosd(Mel)*cosd(Maz); cosd(Mel)*sind(Maz); sind(Mel)];%initial velocity vector in ENU frame +Mspeed0 = 1.e-06;% Missile initial speed in ENU frame (m/sec); must be non-zero for unit vector calculation. +Maz = 0; % Assumed coplanar engagement geometry until a source doc says otherwise. +Mel = 40; % Missile launch elevation (deg). +Mv0 = Mspeed0*[cosd(Mel)*cosd(Maz); cosd(Mel)*sind(Maz); sind(Mel)];%initial velocity vector in ENU frame missileState = [Mp0; Mv0]; % Missile state vector in ENU frame (m. and m/sec.) Moutput_vector(1,:)= [t0,missileState', 0, 0, 0]; missileSpeed(1,:) = [t0,norm(missileState(4:6))]; diff --git a/code/M3dxdt.m b/code/M3dxdt.m index b1fed04..5f626be 100644 --- a/code/M3dxdt.m +++ b/code/M3dxdt.m @@ -1,22 +1,23 @@ function MxDot = M3dxdt(Mt,missileState,targetState) %% Missile is a 120 mm Rocket with an IR Seeker % Compute Missile xDot = [velocity; acceleration]' -g = 9.8066; % acceleration due to gravity -gravity = [0; 0; -g]; % gravity vector -diameter = 0.120; % missile diameter in m. -S = XX; % reference area in m^2 -[acousticSpeed,rho] = getRho(missileState); -Mspeed = norm(missileState(4:6)); -mach = Mspeed/acousticSpeed; -timeCurve = [XX]'; %seconds -thrustCurve = [XX]'; % Newtons -massCurve = [XX]'; % kg - -machCurve = [XX]; -dragCurve = [XX]; - -machCurveB = [XX]; -dragCurveB = [XX]; +g = 9.8066; % acceleration due to gravity +gravity = [0; 0; -g]; % gravity vector +diameter = 0.120; % missile diameter in m. +S = pi*(diameter/2)^2; % reference area in m^2 +[acousticSpeed,rho] = getRho(missileState); +Mspeed = norm(missileState(4:6)); +mach = Mspeed/acousticSpeed; +% Transcribed from the provided missile-data plots. +timeCurve = [0 2.8 2.8001 11.6 11.6001 25]'; %seconds +thrustCurve = [6200 6200 1200 1200 0 0]'; % Newtons +massCurve = [31.8 25.1 25.1 20.9 20.9 20.9]'; % kg + +machCurve = [0 0.5 0.7 0.9 1.0 1.05 1.1 1.2 1.5 2.0 2.5 3.0]; +dragCurve = [0.185 0.180 0.220 0.280 0.330 0.360 0.355 0.330 0.250 0.165 0.110 0.090]; + +machCurveB = [0 0.5 0.7 0.8 0.9 1.0 1.05 1.1 1.2 1.5 2.0 2.5 3.0]; +dragCurveB = [0.150 0.145 0.145 0.160 0.200 0.300 0.315 0.305 0.275 0.205 0.130 0.080 0.045]; if(Mt<=timeCurve(end-1)) Cd = interp1(machCurveB,dragCurveB,mach,'linear','extrap'); diff --git a/code/T3dxdt.m b/code/T3dxdt.m index 1363e6d..40cb373 100644 --- a/code/T3dxdt.m +++ b/code/T3dxdt.m @@ -1,23 +1,24 @@ function TxDot = T3dxdt(t,x) %% Threat is a 240 mm Rocket % Compute Rocket xDot = [velocity; acceleration]' -g = 9.8066; % acceleration due to gravity -gravity = [0; 0; -g]; % gravity vector -diameter = 0.24; % missile diameter in m. -S = XX; % reference area in m^2 -[acousticSpeed,rho] = getRho(x); -mach = norm(x(4:6))/acousticSpeed; - -timeCurve = [XX]'; %seconds -thrustCurve = [XX]'; % Newtons - -massCurve = [XX]'; % kg - -machCurve = [XX]; -dragCurve = [XX]; - -machCurveB = [XX]; -dragCurveB = [XX]; +g = 9.8066; % acceleration due to gravity +gravity = [0; 0; -g]; % gravity vector +diameter = 0.24; % missile diameter in m. +S = pi*(diameter/2)^2; % reference area in m^2 +[acousticSpeed,rho] = getRho(x); +mach = norm(x(4:6))/acousticSpeed; + +% Transcribed from the provided target-data plots. +timeCurve = [0 7.5 7.5001 45]'; %seconds +thrustCurve = [77000 77000 0 0]'; % Newtons + +massCurve = [410 210 210 210]'; % kg + +machCurve = [0 0.5 0.7 0.9 1.0 1.1 1.2 1.3 1.5 2.0 2.5 3.0 3.5 4.0]; +dragCurve = [0.40 0.40 0.43 0.60 0.75 0.87 0.86 0.80 0.74 0.60 0.50 0.44 0.40 0.39]; + +machCurveB = [0 0.5 0.7 0.9 1.0 1.1 1.2 1.3 1.5 2.0 2.5 3.0 3.5 4.0]; +dragCurveB = [0.31 0.31 0.34 0.47 0.62 0.78 0.77 0.72 0.64 0.50 0.40 0.34 0.30 0.27]; if(t<=timeCurve(end-1)) Cd = interp1(machCurveB,dragCurveB,mach,'linear','extrap'); diff --git a/code/getGuidance.m b/code/getGuidance.m index 9b4b9b9..6daabdd 100644 --- a/code/getGuidance.m +++ b/code/getGuidance.m @@ -1,41 +1,44 @@ -function [guideAccel] = getGuidance(Mt,missileState,targetState) -% This function computes the acceleration command to the missile required to -% intercept a target -% Compute proportional navigation (PN) guidance command -PNGain = 0; % PN guidance gain -KVP = XX; -blindRng = XX; %Seeker blind range (m) -acqRng = XX; % Acquisition Range (m) -g = 9.8066; % Acceleration due to gravity (m/sec^2) -gLimit = XX.; % Lateral acceleration limit in gees -relState = targetState' - missileState; -relRng = relState(1:3); -relVel = relState(4:6); -mslVel = missileState(4:6); -% Velocity pursuit guidance -if Mt>XX - nHat = cross(mslVel,relRng); - accelCmdVP = KVP*cross(nHat,mslVel/norm(mslVel))/norm(relRng); -else +function [guideAccel] = getGuidance(Mt,missileState,targetState) +% This function computes the acceleration command to the missile required to +% intercept a target +missileState = missileState(:); +targetState = targetState(:); + +% Compute proportional navigation (PN) guidance command +PNGain = 0; % PN guidance gain +terminalPNGain = 4; +KVP = 0.5; +blindRng = 2; %Seeker blind range (m) +acqRng = 5000; % Acquisition Range (m) +g = 9.8066; % Acceleration due to gravity (m/sec^2) +gLimit = 40.; % Lateral acceleration limit in gees +relState = targetState - missileState; +relRng = relState(1:3); +relVel = relState(4:6); +mslVel = missileState(4:6); +% Velocity pursuit guidance +if Mt>0.5 + nHat = cross(mslVel,relRng); + accelCmdVP = KVP*cross(nHat,mslVel/norm(mslVel))/norm(relRng); +else accelCmdVP = [0; 0; 0]; end if norm(accelCmdVP)>gLimit*g accelCmdVP = gLimit*g*accelCmdVP/norm(accelCmdVP); end -% LOSR is line-of-sight rate -LOSR = cross(relRng,relVel)/dot(relRng,relRng); -if norm(relRng)gLimit*g accelCmd = gLimit*g*accelCmd/norm(accelCmd); -end -guideAccel = accelCmd + accelCmdVP; -disp(norm(accelCmdVP)); -if norm(relRng)