function ZZ= NLLS_LUDENSITY_K_SIM(k,landr,landdable,landd,landb,landc,sumHHzone,HHactual,dimZone,sumEMPzone,BASzone,COMMzone);


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

%k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);  k7=k(8,1); k8=k(9,1); k10=k(10,1); k11=k(11,1); k12=k(12,1); k13=k(13,1); k14=k(14,1); k15=k(15,1); 

%Calculate Lr
for i=1:dimZone
if(landdable(i) < 0.00001); landdable(i) = 0.00001;end;
if(sumHHzone(i) < 0.00001); sumHHzone(i) = 0.00001;end;
if(landb(i) < 0.00001); landb(i) = 0.00001; end;
if(landd(i) < 0.00001); landd(i) = 0.00001; end;
end;

if hhclass==1
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1);
Lr =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzone./sumEMPzone).^k3).*((COMMzone./sumEMPzone).^k4).*((HH_1./sumHHzone).^k5));
elseif hhclass==2
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);
Lr =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzone./sumEMPzone).^k3).*((COMMzone./sumEMPzone).^k4).*((HH_1./sumHHzone).^k5).*((HH_2./sumHHzone).^k6).*sumHHzone) ;
elseif hhclass==3
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);  k7=k(8,1);
Lr =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzone./sumEMPzone).^k3).*((COMMzone./sumEMPzone).^k4).*((HH_1./sumHHzone).^k5).*((HH_2./sumHHzone).^k6).*((HH_3./sumHHzone).^k7).*sumHHzone) ;
elseif hhclass==4
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);  k7=k(8,1); k8=k(9,1);
Lr =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzone./sumEMPzone).^k3).*((COMMzone./sumEMPzone).^k4).*((HH_1./sumHHzone).^k5).*((HH_2./sumHHzone).^k6).*((HH_3./sumHHzone).^k7).*((HH_4./sumHHzone).^k8).*sumHHzone) ;
elseif hhclass==5 
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);  k7=k(8,1); k8=k(9,1); k9=k(10,1);
Lr =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzone./sumEMPzone).^k3).*((COMMzone./sumEMPzone).^k4).*((HH_1./sumHHzone).^k5).*((HH_2./sumHHzone).^k6).*((HH_3./sumHHzone).^k7).*((HH_4./sumHHzone).^k8).*((HH_5./sumHHzone).^k9).*sumHHzone) ;
elseif hhclass==6 
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);  k7=k(8,1); k8=k(9,1); k9=k(10,1); k10=k(11,1);
Lr =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzone./sumEMPzone).^k3).*((COMMzone./sumEMPzone).^k4).*((HH_1./sumHHzone).^k5).*((HH_2./sumHHzone).^k6).*((HH_3./sumHHzone).^k7).*((HH_4./sumHHzone).^k8).*((HH_5./sumHHzone).^k9).*((HH_6./sumHHzone).^k10).*sumHHzone) ;
end    
% Scaling ...
%sum_landr = sum(landr);
%sum_Lr = sum(Lr);
%if(sum_Lr > sum_landr)
%k0_temp = k0.*(sum_landr/sum_Lr)    
%Lr = Lr.*(sum_landr/sum_Lr);
%end;

for i=1:dimZone
    if(Lr(i) < 0.0001) 
        Lr(i) = 0.0001;
    end
end

% Residual calculation
epsilon_k = landr -Lr;
csvwrite('C:\G-LUM Finale\Calibration\Output\epsilon_k.dat',epsilon_k); % Output path must be specified


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


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