function [g,dg]=Recon(Var)
% Optimization objective 
% g  :objective to be minimized
% dg :gradient of g at var
% y measurement vector
% P Measurement covarians matrix (inverted)
% U Measurement incident matrix

GlobalVars;

% Make meas data
% y Measurement vector
% Y Measurement covariance matrix
% P Measurement index matrix

% if Unscaled
%    Sc =0.0004;     % Objective scaling value
% else
%    Sc =0.0001;     % Objective scaling value
% end
% Sc=Sc*0.1;

Sc=0.001;
y=Meas;
P=diag(MeasStd.^-2);
U=MeasU;

% Gaussian distribution
% g=Sc*(U*Sv*Var-y)'*P*(U*Sv*Var-y);
% if nargout > 1
%    dg=Sc*2*(U*Sv)'*P*(U*Sv*Var-y);
% end

% Combined Gaussian distribution
g=0;
p=0.3;b=10; %Tuning parameters
[r,ixy]=find(U==1);
dev=Var(ixy)-Sv(ixy,ixy)\y;
Q=Sv(ixy,ixy)\MeasStd(ixy);
for i=1:length(ixy) 
  g=g-log((1-p)*exp(-0.5*dev(i)^2/Q(i)^2)+ p/b*exp(-0.5*dev(i)^2/(Q(i)^2*b^2)));
end
g=Sc*g;
 
if nargout > 1
    dg=zeros(size(Var,1),1);
    for i=1:length(ixy)
      dg(ixy(i),1)=((b^2-1)*p*dev(i)*exp(0.5*dev(i)^2/Q(i)^2)) ...
                   /(b^2*Q(i)^2*(b*(p-1)*exp(0.5*dev(i)^2/(Q(i)^2*b^2))-p*exp(0.5*dev(i)^2/Q(i)^2))) ...
                   +dev(i)/Q(i)^2;
           end
    dg=Sc.*dg;
end
return         

