function [e,q0,qm,tvec]=tg(A,Q,nt,tmax);
% compute envelope of maximum energy growth using SVD
% 	inputs:
% A: dynamic matrix
% Q: energy measure matrix
% nt: number of time for saving
% tmax: last time
% 	outputs:
% e: envelope of max growth
% q0: worst case initial condition for each time
% qm: optimal response
% tvec: time vector 

%%%% compute exponential matrix for dt
tvec=linspace(0,tmax,nt);
dt=tvec(2)-tvec(1);
expmat=expm(A*dt); % state transition matrix for dt

%%%% factorize energy measure matrix F'*F=Q
[U,S,V]=svd(Q);
s=sqrt(diag(S));
F=diag(s)*U'; % the square root factor F
Finv=U*diag(ones(size(s))./s); % inverse of F

%%%% compute optimal growth
n=size(A,1);
env=zeros(nt,1);env(1)=1;% for saving of max energy
q0=zeros(n,nt); % for saving of worst case initial condition
qm=zeros(n,nt); % for saving of optimal response
H=eye(n);
for ind=2:nt;
  H=H*expmat; % march state transition matrix
  [u,s,v]=svd(F*H*Finv); % SVD decomposition
  q0(:,ind)=v(:,1); % worst case initial condition
  qm(:,ind)=u(:,1); % optimal response
  e(ind)=abs(s(1,1))^2;% max energy
  fprintf('%% it %3.0f / %3.0f %8.3f || %8.3e  \n',ind,nt,tvec(ind),e(ind));
end
q0=Finv*q0;% project back
qm=Finv*qm;
