Jacobian of Massive Equations

4 views (last 30 days)
Ayden Clay on 1 Feb 2021
Edited: Ayden Clay on 1 Feb 2021
As the title states I am trying to find the Jacobian of a few massive equations.
Per the documentation I have set up the problem as follows:
% Define symbols to be used throughout.
syms lat lon u v w h phi theta psi p q r delta1 delta2 delta3 delta4 Hg GM J2 Re e Ixx Ixz Iyy Izz omega m real
% Rotational rate in body axis.
omegabnb = [p;q;r];
% Inertia matrix.
I = [Ixx,0,Ixz;0,Iyy,0;-Ixz,0,Izz];
% Rotational rate of Earth.
omegaeie = [0;0;omega];
RE = Re/sqrt(1-e^2*sin(lat)^2);
RN = (Re*(1-e^2))/((1-e^2*sin(lat)^2)^(3/2));
% Euler orientations.
roll = [1,0,0;0,cos(phi),-sin(phi);0,sin(phi),cos(phi)];
pitch = [cos(theta),0,sin(theta);0,1,0;-sin(theta),0,cos(theta)];
yaw = [cos(psi),-sin(psi),0;sin(psi),cos(psi),0;0,0,1];
% Euler Matrix.
Cbn = simplify(yaw*pitch*roll);
Cnb = Cbn';
% LLA Matrix.
[~,Cne] = nRe(lat,lon);
Cne = simplify(Cne);
Cen = Cne';
% Earth rotation orientation.
[~,Cei] = eRi(Hg);
Cei = simplify(Cei);
Cie = Cei';
% body velocity
vbeb = [u;v;w];
vneb = Cbn'*vbeb;
vneb = simplify(vneb);
% Latitude and longitudinal rates.
latdot = simplify(vneb(1)/(RN+h));
longdot = simplify((vneb(2)/(RE+h))/cos(lat));
omeganen = [longdot*cos(lat);-latdot;-longdot*sin(lat)];
% Forces (outputs six equations determining the various forces and moments, equations in terms of most of the state)
% VERY long equations (126 coefficients).
[Fb,Mb] = MomentFunc;
Rc = [1/(RN+h),0,0;
0,1/((RN+h)*cos(lat)),0 ;
0,0,-1];
% Position calculation.
reeb = [((Re/sqrt(1-e^2*(sin(lat)^2)))+h)*cos(lat)*cos(lon);
((Re/sqrt(1-e^2*(sin(lat)^2)))+h)*cos(lat)*sin(lon);
((Re*(1-e^2)/sqrt(1-e^2*(sin(lat)^2)))+h)*sin(lat)];
riib = Cie*reeb;
Riib = sqrt(riib(1)^2+riib(2)^2+riib(3)^2);
tmpx = (1-5*(riib(3)/Riib)^2)*riib(1);
tmpy = (1-5*(riib(3)/Riib)^2)*riib(2);
tmpz = (3-5*(riib(3)/Riib)^2)*riib(3);
tmp = [tmpx;tmpy;tmpz];
tmp = riib + (3/2)*J2*(Re^2)*tmp/(Riib^2);
% Acceleration due to gravity.
gammaiib = -GM/(Riib^3)*(tmp);
% Latitude, Longitude and Altitude rate.
% Body velocity rate.
vbebdot = Cen*Cnb*(Fb/m) + Cei*gammaiib - Cbn*Cne*SkewMat(Cen*omeganen+Cne*Cnb*omegabnb)*Cen*Cnb*vbeb;
% Angular rotational rate.
omegabnbdot = I\(Mb-SkewMat(Cbn*Cne*omegaeie+Cbn*omeganen+omegabnb)*I*(Cbn*Cne*omegaeie+Cbn*omeganen+omegabnb))-SkewMat(-omegabnb)*Cbn*omeganen;
% Roll rate.
phidot = p*q*sin(phi)*tan(theta)+r*cos(phi)*tan(theta);
% Pitch rate.
% Yaw rate.
psidot = q*sin(phi)*sec(theta)+r*cos(phi)*sec(theta);
% Control deflection rate (DUE TO STATE, NOT DUE TO CONTROL).
% Jacobian for linearization.
This progam runs for several days and then ends with an out of memory error. I wish I'd stored the error, but I will have it once it finishes running again.
For reference I have 64GB ram.
Unclear above is the fact that vbebdot and omegabnbdot are HUGE equations (spanning several pages).