function ZZ=NLLS_LUDENSITY_G_SIM(g,landr,landdable,landd,landb,sumEMPzone,BASzone,COMMzone,dimZone,HHactual,sumHHzone)

hhclass=size(HHactual,2);
%Calculate HH by type
if hhclass>0 HH_1=HHactual(:,1); end
if hhclass>1 HH_2=HHactual(:,2); end
if hhclass>2 HH_3=HHactual(:,3); end
if hhclass>3 HH_4=HHactual(:,4); end
if hhclass>4 HH_5=HHactual(:,5); end
if hhclass>5 HH_6=HHactual(:,6); end
%Initialization
%g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);g6=g(7,1);g7=g(8,1);g8=g(9,1);g9=g(10,1);

%Calculate Lb
for i=1:dimZone
if(landdable(i) < 0.001) landdable(i) = 0.001;end;
if(sumEMPzone(i) < 0.001) sumEMPzone(i) = 0.001;end;
if(landd(i) < 0.001) landd(i) = 0.001;end;
end;

if hhclass==1
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1);
Lb = g0.*(landd./landdable).^g1.*(BASzone./sumEMPzone).^g2.*((COMMzone./sumEMPzone).^g3).*((HH_1./sumHHzone).^g4).*BASzone;
elseif hhclass==2
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);
Lb = g0.*(landd./landdable).^g1.*(BASzone./sumEMPzone).^g2.*((COMMzone./sumEMPzone).^g3).*((HH_1./sumHHzone).^g4).*((HH_2./sumHHzone).^g5).*BASzone;
elseif hhclass==3
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);g6=g(7,1);
Lb = g0.*(landd./landdable).^g1.*(BASzone./sumEMPzone).^g2.*((COMMzone./sumEMPzone).^g3).*((HH_1./sumHHzone).^g4).*((HH_2./sumHHzone).^g5).*((HH_3./sumHHzone).^g6).*BASzone;
elseif hhclass==4
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);g6=g(7,1);g7=g(8,1);
Lb = g0.*(landd./landdable).^g1.*(BASzone./sumEMPzone).^g2.*((COMMzone./sumEMPzone).^g3).*((HH_1./sumHHzone).^g4).*((HH_2./sumHHzone).^g5).*((HH_3./sumHHzone).^g6).*((HH_4./sumHHzone).^g7).*BASzone;
elseif hhclass==5
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);g6=g(7,1);g7=g(8,1);g8=g(9,1);
Lb = g0.*(landd./landdable).^g1.*(BASzone./sumEMPzone).^g2.*((COMMzone./sumEMPzone).^g3).*((HH_1./sumHHzone).^g4).*((HH_2./sumHHzone).^g5).*((HH_3./sumHHzone).^g6).*((HH_4./sumHHzone).^g7).*((HH_5./sumHHzone).^g8).*BASzone;
elseif hhclass==6
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);g6=g(7,1);g7=g(8,1);g8=g(9,1);g9=g(10,1);
Lb = g0.*(landd./landdable).^g1.*(BASzone./sumEMPzone).^g2.*((COMMzone./sumEMPzone).^g3).*((HH_1./sumHHzone).^g4).*((HH_2./sumHHzone).^g5).*((HH_3./sumHHzone).^g6).*((HH_4./sumHHzone).^g7).*((HH_5./sumHHzone).^g8).*((HH_6./sumHHzone).^g9).*BASzone;
end 
% Scaling the values.
%if(sum(Lb) ~= sum(landb)) 
%    Lb = Lb.*(sum(landb)/sum(Lb));
%end;

for i=1:dimZone
    if(Lb(i) < 0.001) 
        Lb(i) = 0.001;
    end
end
% Residual calculation ..
epsilon_g = landb - Lb;
csvwrite('C:\G-LUM Finale\Calibration\Output\epsilon_g.dat',epsilon_g); % Output path must be specified

%Likelihood calculation ...
%L = 0;
%for i=1:dimZone;
%    L = L + (landb(i)*log(Lb(i)));
%end;

%ZZ = -L
Z=(Lb-landb).^2;
ZZ=sum(Z,1)./(2*dimZone)