equilmoments_jd

0001 function [pshift,pelong,ptrian,pli] = equilmoments_jd(theta,ptx,pty,mfac)
0002 %
0003 if nargin < 4,
0004     mfac = 50;
0005 end
0006 %
0007 nt = length(theta);
0008 %
0009 if any(theta ~= linspace(0,2*pi,nt)),
0010     error('case of nonuniform theta grid not yet implemented')
0011 end
0012 %
0013 pmx = interpft(ptx(:,1:nt-1),(nt-1)*mfac,2); %FFT interpolation
0014 pmy = interpft(pty(:,1:nt-1),(nt-1)*mfac,2); %FFT interpolation
0015 %
0016 pxmin = min(pmx,[],2);
0017 pxmax = max(pmx,[],2);
0018 pymin = min(pmy,[],2);
0019 pymax = max(pmy,[],2);
0020 %
0021 pli = (pxmax - pxmin)/(pxmax(end) - pxmin(end));
0022 pshift = (pxmax + pxmin)/2;
0023 pelong = (pymax - pymin)./(pxmax - pxmin);
0024 %
0025 pmymin = repmat(pymin,[1,(nt-1)*mfac]);
0026 pmymax = repmat(pymax,[1,(nt-1)*mfac]);
0027 %
0028 pxymin = sum(pmx.*(pmy == pmymin),2)./sum(pmy == pmymin,2);
0029 pxymax = sum(pmx.*(pmy == pmymax),2)./sum(pmy == pmymax,2);
0030 %
0031 ptrian = -(pxymin + pxymax - 2*pshift)./(pxmax - pxmin);
0032 
0033 
0034 
0035