% Function Evaluation
function ZZ=NLLS_LUDENSITY_P_SIM(p,landr,landdable,landd,landc,sumEMPzone,BASzone,COMMzone,dimZone,HHactual,sumHHzone)


%Initialization
%p0=p(1,1); p1=p(2,1); p2=p(3,1); p3=p(4,1); p4=p(5,1); p5=p(6,1);p6=p(7,1); p7=p(8,1); p8=p(9,1); p9=p(10,1);
%Calculate Lc

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
%Calculate Lc

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
    p0=p(1,1); p1=p(2,1); p2=p(3,1); p3=p(4,1); p4=p(5,1);
Lc = p0.*(landd./landdable).^p1.*(BASzone./sumEMPzone).^p2.*((COMMzone./sumEMPzone).^p3).*((HH_1./sumHHzone).^p4).*COMMzone;
elseif hhclass==2
Lc = p0.*(landd./landdable).^p1.*(BASzone./sumEMPzone).^p2.*((COMMzone./sumEMPzone).^p3).*((HH_1./sumHHzone).^p4).*((HH_2./sumHHzone).^p5).*COMMzone;
elseif hhclass==3
Lc = p0.*(landd./landdable).^p1.*(BASzone./sumEMPzone).^p2.*((COMMzone./sumEMPzone).^p3).*((HH_1./sumHHzone).^p4).*((HH_2./sumHHzone).^p5).*((HH_3./sumHHzone).^p6).*COMMzone;
elseif hhclass==4
Lc = p0.*(landd./landdable).^p1.*(BASzone./sumEMPzone).^p2.*((COMMzone./sumEMPzone).^p3).*((HH_1./sumHHzone).^p4).*((HH_2./sumHHzone).^p5).*((HH_3./sumHHzone).^p6).*((HH_4./sumHHzone).^p7).*COMMzone;
elseif hhclass==5
Lc = p0.*(landd./landdable).^p1.*(BASzone./sumEMPzone).^p2.*((COMMzone./sumEMPzone).^p3).*((HH_1./sumHHzone).^p4).*((HH_2./sumHHzone).^p5).*((HH_3./sumHHzone).^p6).*((HH_4./sumHHzone).^p7).*((HH_5./sumHHzone).^p8).*COMMzone;
elseif hhclass==6
Lc = p0.*(landd./landdable).^p1.*(BASzone./sumEMPzone).^p2.*((COMMzone./sumEMPzone).^p3).*((HH_1./sumHHzone).^p4).*((HH_2./sumHHzone).^p5).*((HH_3./sumHHzone).^p6).*((HH_4./sumHHzone).^p7).*((HH_5./sumHHzone).^p8).*((HH_6./sumHHzone).^p9).*COMMzone;
end

% Only for likelihood functions ..
if(sum(Lc) > sum(landc) )
    p_temp = p0.*(sum(landc)/sum(Lc))
     Lc = Lc.*(sum(landc)/sum(Lc));
end;

for i=1:dimZone
    if(Lc(i) < 0.001) 
        Lc(i) = 0.001;
    end
end

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


%Likelihood calculation 
%L = 0;
%for i=1:dimZone;
%    L = L + (landc(i)*log(Lc(i)));
%end;

% TELUM Objective

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