function [pod,pode]=mypod(R,Q,n);
% compute n first POD modes of covariance matrix 
%	inputs:
% R: pointwise covariance matrix
% Q: inner product matrix
% n: number of desired POD modes
% 	outputs:
% pod: matrix of first n POD modes (as columns)
% pode: vector of variance of POD amplitude 

%%%% eigenmode decomposition of R*Q
[pod,pode]=eig(R*Q);
pode=abs(diag(pode));

%%%% sort in decreasing mean energy and select n first
[t,o]=sort(-real(pode));
pode=pode(o(1:n));
pod=pod(:,o(1:n));
