function [Rabbit_Location,Rabbit_Energy] = Harris_Hawks_Optimization1 (fobj,N,Lower,Upper,Maxiter)
## fx input function
N=length(Lower);
## N number of individulas in group
## Lower: Lower limits (xmin) of independent variables
## Upper: Upper limits (xmax) of independent variables
## Maxiter: Maximum number of iterations
## fobj : onjective function 
D=length(Lower);
X=zeros(N,D);
XNEW=zeros(N,D);
X1=zeros(1,D); 
X2=zeros(1,D);
X_rand=zeros(1,D);
fitnessX=zeros(1,N);
Rabbit_Location=zeros(1,D);
mean2Dmat=zeros(1,D);
lvyout1=zeros(1,D);
lvyout2=zeros(1,D);
 
%% Assigning an initial value to the Rabbit Energy
%% This value will be iteratively decreased 
%% until the termination criterion is met 
Rabbit_Energy=1e+200;

%% Constructing the initial population

for i=1:N
  for j=1:D
    X(i,j)=Lower(j)+(Upper(j)-Lower(j))*rand;
  endfor
endfor

iter=1;

while(iter <= Maxiter)
     
     X=boundary(X, Lower,Upper,N,D); %% Applying boundary check mechanism 
     %% Solution update
     %% Updating the rabbit location (best solution) 
     %% and Rabbit Eenrgy (its respective fitness value)
     for i=1:N
       fitnessX(i)=fobj(X(i,:));
       if (fitnessX(i)<Rabbit_Energy)
         Rabbit_Energy=fitnessX(i);
         Rabbit_Location=X(i,:);
       end
     endfor
     
     E0=2*(1-(iter/Maxiter)); 
     
     for i=1:N
        J=2*rand-1; %% Calculating the jumping strength
        Escaping_Energy=E0*J; %% Updating rabbit escape energy
        
        if(abs(Escaping_Energy)>=1.0) %% Exploration phase
        
           q=rand;
           rand_index=ceil(N*rand);
           X_rand=X(rand_index,:);
           if (q<0.5)
              X(i,:)=X_rand-rand.*abs(X_rand-2.0.*rand.*X(i,:));
           else
              mean2Dmat=mean(X);
              X(i,:)=(Rabbit_Location-mean2Dmat)-(rand.*(Lower+((Upper-Lower).*rand)));
           end
        
        elseif(abs(Escaping_Energy)<1.0) %% Exploitation phase
            r=rand;
            if((r>=0.5)&&(abs(Escaping_Energy)<0.5))
               X(i,:)=Rabbit_Location-Escaping_Energy.*abs(Rabbit_Location-X(i,:)); %% Employing Hard besiege attacking strategy by Eq
           
            elseif ((r>=0.5)&&(abs(Escaping_Energy)>=0.5))  %% Employing Hard besiege attacking strategy by Eq
               X(i,:)=(Rabbit_Location-X(i,:))-Escaping_Energy.*abs(2.0.*(1.0-rand).*Rabbit_Location-X(i,:));
            
            elseif ((r<0.5) && (abs(Escaping_Energy)>=0.5))%% Employing soft besiege with rapid dive attacking strategy by Eq
               X1=Rabbit_Location-Escaping_Energy.*abs(2.0*(1.0-rand).*Rabbit_Location-X(i,:));
               if (fobj(X1) < fobj(X(i,:)))
                   X(i,:)=X1;
               else
                   lvyout1=Levy(D);
                   X2=Rabbit_Location-Escaping_Energy.*abs( 2.0*(1.0-rand).*Rabbit_Location-X(i,:)) + rand.*lvyout1;
                   if (fobj(X2) < fobj(X(i,:)))
                       X(i,:)=X2;
                   end
               end
           elseif ((r<0.5)&&(abs(Escaping_Energy)<0.5)) %% Employing hard besiege with rapid dive attacking strategy by Eq
                mean2Dmat=mean(X);
                X1=Rabbit_Location-Escaping_Energy.*abs(2.0*(1.0-rand).*Rabbit_Location-mean2Dmat);
                if (fobj(X1) < fobj(X(i,:)))
                   X(i,:)=X1;
                else
                   lvyout2=Levy(D);
                   mean2Dmat=mean(X);
                   X2=Rabbit_Location-Escaping_Energy.*abs( 2.0*(1.0-rand).*Rabbit_Location-mean2Dmat) + rand.*lvyout2;
                      if (fobj(X2) < fobj(X(i,:)))
                           X(i,:)=X2;
                      end
                end
           end
       end   
        
   end
    iter=iter+1;
end
endfunction
