% Function Evaluation
function ZZ=NLLS_LUDENSITY_P(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);


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

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