[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
help with adams.m -> lsode conversion
From: |
Haisam K. Ido |
Subject: |
help with adams.m -> lsode conversion |
Date: |
Fri, 25 Jul 1997 14:01:54 +0000 (GMT) |
I'm trying to convert the adams ODE integrator call in matlab
[t_int,x_int] = adams('xprime',[t_k;t_k1],x_in,1e-10);
In a file called od_main.m (below), to the lsode one below with no success
[t_int,x_int] = lsode( xprime, x_0, [t_k;t_k1] );
xprime has this:
function [sys,x0] = xprime(t,x,u,flag)
However, I get the following:
Initializing variables.
Loading observation data.
Iteration number is ...
mloop = 1
error: `flag' undefined near line 11 column 8
error: evaluating expression near line 11, column 8
error: evaluating argument list element number 0
error: evaluating index expression near line 11, column 4
error: evaluating binary operator `==' near line 11, column 13
error: if: error evaluating conditional expression
error: evaluating if command near line 11, column 1
error: called from `xprime' in file
`/disk3/users/idoh/projects/Kalman/kf_labs/od/xprime.m'
error: evaluating argument list element number 0
error: near line 122 of file
`/disk3/users/idoh/projects/Kalman/kf_labs/od/od_main.m'
od_main.m is below:
%
% OD_MAIN.M
% Orbit determination using batch least squares.
%
% Initialize variables.
clear
%close all
format short;
disp(' Initializing variables.')
const;
init;
summary=0;
% Read in observation data.
x_k1 = x_k;
phi_k1 = phi_k;
else %if nloop~=1,
x_in = [x_k; reshape(phi_k,225,1)];
% [t_int,x_int] = adams('xprime',[t_k;t_k1],x_in,1e-10);
[t_int,x_int] = lsode( xprime, x_0, [t_k;t_k1] );
[nsize,msize] = size(t_int);
x_k1 = x_int(nsize,1:15)';
phi_k1 = reshape(x_int(nsize,16:240)',15,15);
end
% Solve for H-tilda and calculated range.
[h_tilda,rho_calc] = h_calc(x_k1,t_k1,sta_id);
% Calculate measurement matrix, H, for this iteration.
h_k1 = h_tilda * phi_k1;
% Calculate observation error.
delta_y = rho_k1 - rho_calc;
summary = [summary
delta_y];
% Accumulate least-square sums.
sum1 = sum1 + h_k1' * inv(r_obs) * h_k1;
sum2 = sum2 + h_k1' * inv(r_obs) * delta_y;
sum_rms = sum_rms + delta_y^2;
% Update variables for next iteration. End inner loop.
phi_k = phi_k1;
x_k = x_k1;
t_k = t_k1;
end
% Compute state error.
[n_sum,m_sum] = size(sum1);
[u_svd,s_svd,v_svd] = svd( inv(p_0) + sum1);
delta_x = zeros(n_sum,1);
for i_sum=1:n_sum,
delta_x = delta_x + ( u_svd(:,i_sum)' * sum2 / s_svd(i_sum,i_sum) )
*...
v_svd(:,i_sum);
end
% Check for convergence. (If desired, can keep iterating until
% rcon is less than some small value.)
rcon = sqrt(sum_rms/i_size)
sum_plot = [sum_plot summary];
summary=[0];
% Update epoch state.
x_k = x_0 + delta_x;
x_0 = x_k;
% Reset variables for next iteration. End outer loop.
phi_k = phi_0;
t_k = t_0;
sum1 = sum1_init;
sum2 = sum2_init;
sum_rms = 0;
end
%end