%Simple Matlab interface for Mor4ansys;
%ypliu26 at gmail.com;
%Developed by Yongpan Liu
%EE Dept., Tsinghua University, P.R.China.
%MOR2MATLAB function is adopted from the simulink interface package
%by Pascal Maglie from ETH.

%read matrix from files;
[A,rows,cols,entries]=MOR2MATLAB('mor.A');
[B,rows,cols,entries]=MOR2MATLAB('mor.B');
[E,rows,cols,entries]=MOR2MATLAB('mor.E');
[C,rows,cols,entries]=MOR2MATLAB('mor.C');
[fix,rows,cols,entries]=MOR2MATLAB('fixstep.txt');

%global matrix transform for differential equation;
Ein = inv(E);
AA = Ein*A;
BB = Ein*B;
global Am;
global Bm;
Am = AA;
Bm = BB;

%transient analysis;
y0=zeros(30,1); %initial temperature; 
atol=ones(30,1)*1e-5;
options = odeset('RelTol',1e-3,'AbsTol',atol);
%[t,Y] = ode45(@thrusys,[0 20],y0,options); %very slow;
[t,Y] = ode15s(@thrusys,[0 20],y0,options); %very fast;

%plot figure;
T=C*Y'; %turn down the output dimension by ode15 call; 
To = fix;
to = To(1,:);
figure(1);
plot(t,T(1,:)','-g');
hold on;
plot(to,To(2,:)','-r');
xlabel('time (s)');
ylabel('temperature (K)')
legend('mor30-FuelTop','ansys-FuelTop',4);
title('mor30 vs ansys transient results');
axis([0 20 150 500]);

figure(2);
plot(t,T(2,:)','-.g');
hold on;
plot(to,To(3,:)','-.r');
xlabel('time (s)');
ylabel('temperature (K)')
legend('mor30-FT-200','ansys-FT-200',4);
title('mor30 vs ansys transient results');
axis([0 20 0 350]);

figure(3);
plot(t,T(3,:)',':g');
hold on;
plot(to,To(4,:)',':r');
xlabel('time (s)');
ylabel('temperature (K)')
legend('mor30-FT-400','ansys-FT-400',4);
title('mor30 vs ansys transient results');
axis([0 20 0 300]);

figure(4);
plot(t,T(4,:)','--g');
hold on;
plot(to,To(5,:)','--r');
xlabel('time (s)');
ylabel('temperature (K)')
legend('mor30-FuelBot ','ansys-FuelBot ',4);
title('mor30 vs ansys transient results');
axis([0 20 0 12]);