function [mstnew, vstnew] = upd_marg_norm_moments(mstold, ...
vstold, mstsp, vstsp, inda, indf, xt, ti, gm, nodes, ndps, w, mthd)

if nargin < 11, ndps =[]; end
% if nargin < 12, w =[]; end
if nargin < 13, mthd = 1; end

global pi2 f0 nf e ubnd;

% step B: compute the conditional distr of amplitude a given freq. f
    % assume amp a is indep of freq. f
    %    matf = mstold(inda) ; 
    %    vatf = vstold(inda, inda) ;
    
    % assume amp dep on freq. f
    mftsp =  mstsp(indf);
    vftsp = vstsp(indf, indf);
    mat = mstold(inda) ; 
    mft =  mstold(indf);
    vft = vstold(indf, indf);
    vaft = vstold(inda, indf); % in my notation, vaft is cross-covariance; 
                               % while vatf is conditional cov.
    vafivf = vaft * inv(vft);    
    % compute the conditional variance of a
    vatf = vstold(inda, inda) - vafivf * vaft';
    
 % Step C: the nodes of quadrature rules or sampling points
 if mthd == 0
     chovf = chol(vftsp)';
     chovfnd = chovf*ndps;
     ndf = repmat(mftsp,1,nodes) + chovfnd; % is a matrix. each column  corresponds to a node
     nj = 0;
       for i = 1:nodes
           if chovfnd(:,i) > pi | chovfnd(:,i) < -pi
               ndf(:,i) = f0;
               nj = nj+1;
           end
       end
       % ndf = sort(ndf);
       if nj == nodes,
           error('All nodes are pruned by QR');
       end
   % Step B (Continue)
   % compute the conditional mean of a
   matf = repmat(mat,1,nodes) +  ...
     vafivf*( ndf - repmat(mft,1,nodes) );
   wpr = colize(gauss_pdf(ndf, mft, vft) ...
          ./ gauss_pdf(ndf, mftsp, vftsp) );
   w = w .* wpr;
elseif mthd == 1
      % vftsp = vft; mftsp = mft;
       ndf0 = normrnd_mt(nodes, f0, vftsp);
       ndf = ndf0 + repmat(mftsp,1,nodes);
       mstpti = mstsp(3)*ti*2*pi;
       intvl = pi;
       if 0
           if mstpti > intvl
               nlow = -intvl+mstsp(2);
               nhigh = intvl+mstsp(2);
           else
               nlow = mstsp(2)-mstpti;
               nhigh = 2*intvl+mstsp(2)-mstpti;
           end  
           % nhigh = intvl+mstsp(2);
       else
           nlow = -intvl+mstsp(indf);
           nhigh = intvl+mstsp(indf);
         %  nlow = -intvl+mft; nhigh = intvl+mft;
       end
       
       g =0; null_ind=[];        
       for i = 1:nodes 
           if max(ndf(:,i) < nlow) | max( ndf(:,i) > nhigh)
               g = g+1;
               null_ind(g) = i;
           end
       end
       ndf(:,null_ind) = [];
       nodes = nodes - g;
       if nodes == 0, error('alll samples rejected');end
      % ndf = sort(ndf);
       w(null_ind) = [];
       w = colize(exp( gauss_logProb(ndf, mft, vft) ...
         - gauss_logProb(ndf, mftsp, vftsp)) );
       % ndf0 = tttrnd(nodes, 3, mftsp, vftsp);
       % ndf0prob = tmpdf(ndf0, 3,mftsp, vftsp);
       % w = colize(gauss_pdf(ndf0, mftsp, vftsp) ./ ndf0prob);
   
       % Step B (Continue)
       % compute the conditional mean of a
       matf = repmat(mat,1,nodes) + ...
          vafivf*(ndf-repmat(mft,1,nodes));     
elseif mthd == 3
     mstpti = mstsp(3)*ti;
     intvl = pi;
     if mstpti > intvl
         ndf =linspace(-intvl+mstsp(2),...
          intvl+mstsp(2), nodes);    
     else
        ndf =linspace(mstsp(2)-mstpti,...
           mstsp(2)-mstpti+2*intvl, nodes);    
     end         
     % Step B (Continue)
     % compute the conditional mean of a
     matf = repmat(mat,1,nodes) + ...
      vafivf*(ndf-repmat(mft,1,nodes));     
end
   
%   h = obs_h(ndf, ti);
h = sin(ndf);
  
mxt = sum(h .* matf,1)'; % same as  mxt(j) = h(:,j)' * matf(:,j); % is a column vector, each element corresponds to a node
vxt = sum( (vatf * h) .* h, 1)' + gm;
   
% compute the normalization factor and the mean by quard rule.
nmpdfxi = normpdf(xt, mxt, sqrt(vxt)); % a col vector

mfig = 0;
if mthd ~=3
   c = w' * nmpdfxi;  
end
if mthd ==3 | mfig == 1
    a=[]; d=[];
    for j = 1:nodes
        a(j) =  gauss_pdf(ndf(:,j), mft, vft );
        d(j) =  nmpdfxi(j) * a(j);
    end
   if mfig == 1
        figure(2),subplot(411),plot(w.*nmpdfxi);
              subplot(412),plot(w);
              subplot(413),plot(d);              
              subplot(414),plot(a);              
              pause;
   end
   if mthd ==3
       c = sum(d);
   end
end

   vatfh = vatf * h; % is matrix. each column  corresponds to a node
   
   K = scale_cols(vatfh, 1./vxt); % each column corresponds to a node
   Kdx = scale_cols(K, xt - mxt );
   
  % the conditional mean of a given x and f 
   umat = matf + Kdx; % each column corresponds to a node
   
   % compute the mean of a and b
   wnmc = w .* nmpdfxi / c;
%   mst(inda,i) =  umat * wnmc;  mst(indf,i) =  ndf * wnmc;
   mstnew =  [umat;ndf] * wnmc; 
   
  % compute the variance of st, i.e., a and b
   uvaasum = zeros(nf);
   uvffsum = zeros(nf,nf);
   if 1        
       uvafsum = zeros(nf,nf);
       for j = 1:nodes
         uvat  = (e - K(:,j)*h(:,j)') * vatf + umat(:,j)*umat(:,j)';
         uvaasum = wnmc(j) * uvat   + uvaasum;
         uvafsum = wnmc(j)* umat(:,j)*ndf(:,j)' + uvafsum;     
         uvffsum = wnmc(j)* ndf(:,j)* ndf(:,j)' + uvffsum;
      end
      vstnew = [uvaasum uvafsum; uvafsum' uvffsum]- mstnew * mstnew';
   else        
       uvafsum = zeros(nf,nf);
       for j = 1:nodes
         uvat  = (e - K(:,j)*h(:,j)') * vatf + umat(:,j)*umat(:,j)';
         uvaasum = wnmc(j) * uvat   + uvaasum;
         uvffsum = wnmc(j)* ndf(:,j)* ndf(:,j)' + uvffsum;
       end
       vstnew = [uvaasum-mstnew(inda)*mstnew(inda)' uvafsum; uvafsum' uvffsum-mstnew(indf)*mstnew(indf)'];
    end
    
    if min(eig(vstnew)) <= 0
       warning('vstnew is not positive definite');
       mstnew = mstsp([inda;indf]);
       vstnew = vstsp([inda;indf],[inda;indf]);
   end
