function [nli,mwi]=pseudint(w,f,kmax,kmin,ci,weight)
%PSEUDINT is a utility function for FPSEUDO.
%       [NLI,MWI]=PSEUDINT(W,F,kmax,kmin,ci,Weight) accumulates
%       the frequency responses and builds the
%       required real symmetric matrices.
%       F is the MVFR matrix, W is the associated freq. vector.
%       kmax is the maxinum order of compensator elements.
%       kmin is the minimum order of compensator elements.
%       ci is the column to be minimized.
%       Weight is the frequency dependent weighting.
%       NLI is the real symmetric denominator matrix.
%       MWI is the real symmetric numerator matrix.

%       Dr M.P. Ford 2nd September 1987
% Copyright (c) 1987 by GEC Engineering Research Centre & Cambridge Control Ltd

[mf,nf]=fsize(w,f);
lw=length(w);
j=sqrt(-1);
w=w(:);    % leave the j out until forming the final matrix
wt=[];
for i=2*kmax:-1:2*kmin   % do powers of all frequencies at once
   wt=[wt,w.^i];
end
[temp,lwt]=size(wt);
% form Li and Wi
li=zeros(mf,mf);
li(ci,ci)=1;
wi=diag(ones(1,mf));
wi(ci,ci)=0;

% set up nsum,msum to hold integrals
% set up ni,mi to hold values at current frequency
% set up nlast,mlast to hold last values
nsum=zeros(nf,nf*lwt);
msum=nsum;
ni=nsum;
mi=nsum;
nlast=nsum;
mlast=nsum;

intw=w;
%       Comment the next line out if you want to
%       integrate over a linear frequency range.
intw=log10(w);           % integrate over a log scale

intw=diff(intw)./2;      % form differences/2

k=1:mf;
p=1:nf;
fm=f(k,:);
nli=fm'*(li.*weight(1))*fm;
mwi=fm'*(wi.*weight(1))*fm;
%  First point
for q=1:lwt
   nlast(:,(q-1)*nf+p)=wt(1,q).*nli;
   mlast(:,(q-1)*nf+p)=wt(1,q).*mwi;
end
if lw>1    %   Then integrate
    for i=2:lw           % for each frequency
	fm=f((i-1)*mf+k,:);
	nli=fm'*(li.*weight(i))*fm;
	mwi=fm'*(wi.*weight(i))*fm;
	for q=1:lwt
	    ni(:,(q-1)*nf+p)=wt(i,q).*nli;
	    mi(:,(q-1)*nf+p)=wt(i,q).*mwi;
	end
    %  integrate S*N and S*M using trapezoidal
	nsum=nsum+(nlast+ni).*intw(i-1);
	msum=msum+(mlast+mi).*intw(i-1);
	nlast=ni;
	mlast=mi;
    end
else       %  Else use first point
    nsum=nlast;
    msum=mlast;
end
% form the full real symetric matrix
snli=kmax-kmin+1;      % *nf = size of full square matrix
nli=zeros(nf*snli,nf*snli);
mwi=nli;
% put in the off-diagonals first
% work out sign of (1,1) matrix
j11=1;                 % conj(j^kmax)*(j^kmax) is always +1
% first off diagonal above is -j * j11
for m=1:snli
jmn=j11;            % all diagonals have the same sign
   for n=(m+1):snli
      jmn=-j*jmn;   %  moving along rows multiplies by 1/j = -j
      nli((m-1)*nf+p,(n-1)*nf+p)=real(jmn*nsum(:,((m-1)*2+n-m)*nf+p));
      mwi((m-1)*nf+p,(n-1)*nf+p)=real(jmn*msum(:,((m-1)*2+n-m)*nf+p));
   end  % for n
end  % for m
%  form lower part of nli and mwi
nli=nli+nli';
mwi=mwi+mwi';
% now put in the diagonal matrices
for m=1:snli
    nli((m-1)*nf+p,(m-1)*nf+p)=real(j11*nsum(:,((m-1)*2)*nf+p));
    mwi((m-1)*nf+p,(m-1)*nf+p)=real(j11*msum(:,((m-1)*2)*nf+p));
end  % for m