function [ imc, imr ] = imc( lambdaCt, kpn_x, kpn_y, betaC, TRagR )
% Doloï¿½itev ciljne intenzitete in pripadajoï¿½ega potresnega scenarija

coder.varsize('lambdaCa1');
coder.varsize('imCtM');
coder.varsize('Htarget');
coder.varsize('SaTarget');

x = kpn_x; %intenzitete
y = kpn_y; %frekvence

% Izraï¿½un mediane ciljne intenzitete z numeriï¿½no integracijo enaï¿½be tveganja
dim  = 0.001;
im1  = dim:dim:max(x);
im2  = dim/2:dim:max(x)+dim;
%Hint = exp(pchip_(log(x),log(y),log(im2)));
%Hint = exp(pchip_(log(x),log(y),log(im2')));
Hint=exp(interp1(log(x),log(y),log(im2),'linear','extrap'));
m    = abs(diff(Hint)./diff(im2));

extrap   = 0;
j        = 1;
imCtM    = [1];
imCheck  = 1;
lambdaCa = [];
err      = [];
imCt     = 0;
%%
while imCheck>0.001
    imMax = max(x);
    PC1   = logncdf_(imMax,log(imCtM(j)),betaC);
    while PC1<0.99
        imMax  = imMax+1;
        PC1    = logncdf_(imMax,log(imCtM(j)),betaC);
        extrap = 1;
        im1    = dim:dim:imMax;
        im2    = dim/2:dim:imMax+dim;
        %Hint   = exp(pchip_(log(x),log(y),log(im2)));
        %m      = abs(diff(Hint)./diff(im2));
%         im2a=max(x)+dim/2:dim:imMax+dim;
         Hint=exp(interp1(log(x),log(y),log(im2),'linear','extrap'));
         m=abs(diff(Hint)./diff(im2));
%         m2=[m ma];
    end
%     if extrap==
%         m2=m;
%     end
    lambdaCa1 = [];
    PC        = zeros(length(im1),1);
    
    for i=1:length(im1)
        PC(i)     = logncdf_(im1(i),log(imCtM(j)),betaC);
        lambdaCa1 = [lambdaCa1,PC(i)*m(i)*dim];
        %lambdaCa1(i)=PC(i)*m2(i)*dim;
    end
    
    lambdaCa = [lambdaCa,sum(lambdaCa1)];
    %lambdaCa(j)=sum(lambdaCa1);
    % lambdaCa=PC*m'*dim;

    err = [err, lambdaCa(j)-lambdaCt];
    if j==1 || abs(sum(sign(err)))==j
        if err(j)>0
            imCtM = [imCtM, imCtM(j)*2];
        elseif err(j)<0
            imCtM = [imCtM, imCtM(j)/2];
        else
            imCtM = [imCtM, imCtM(j)];
        end
    else
        if err(j)>0
            imCtM = [imCtM,imCtM(j)+abs((imCtM(j)-imCtM(j-1)))/2];
        elseif err(j)<0
            imCtM = [imCtM,imCtM(j)-abs((imCtM(j)-imCtM(j-1)))/2];
        else
            imCtM = [imCtM,imCtM(j)];
        end
    end
    imCheck = abs(imCtM(j+1)-imCtM(j));
    j       = j+1;
end

imc=imCtM(end);

indagR_ = find(y<1/TRagR,1,'first');
indagR = indagR_(1,1);
if indagR > 1
    imr_ = interp1(y([indagR-1 indagR]), x([indagR-1 indagR]), 1/TRagR);    
    imr = imr_(1,1);
else 
    imr = x(1);
end


%
% %
% kvanAll = zeros(nAcc,1);
% for i=1:nAcc
%    kvanAll(i) = (i-0.5)/nAcc;
% end
% [~,index] = min(abs(kvanAll-0.1587));
% 
% % Izracun ciljne intenzitete za izbran kvantil (0.1625 za skupino s 40 akcelerogrami)
% ss = size(imCtM,2);
% if ss ~= 1
%     imCt = imCtM(ss)*exp(norminv_(kvanAll(index),0,1)*betaC);
% end
% 
% % Doloï¿½itev epsilona za pogojni spekter
% %[median_sa,sigma]=SP_1996(T1,M,R,1,0,4);
% %epsilon=(log(imCt16)-log(median_sa))/sigma;
% 
% % Doloï¿½itev povratne dobe, ki pripada ciljni intenziteti
% %RP=1/exp(pchip_(log(x),log(y),log(imCt16)));
% 
% %Himct
% Himct  = exp(pchip_(log(x),log(y),log(imCt)));
% PCimct = logncdf_(imCt,log(imCtM(ss)),betaC);
% 
% saUser         = kpn_x';
% hUser          = kpn_y';
% 
% isExtrap = extrap;
% if(isExtrap == 1)
% 	% Shorten interpolated HF to max 200 points (or greater than user defined function)
% 	while(size(Hint,2) > 200 && size(Hint,2) > size(hUser,2))
% 		Hint_tmp = Hint(2:2:end-1);
% 		im2_tmp  = im2(2:2:end-1);
% 		Hint     = [Hint(1,1) Hint_tmp Hint(1,end)];
% 		im2      = [im2(1,1) im2_tmp im2(1,end)];
%     end
% 
% 	% Modify interpolated Hazard function to include(Sact)
% 	if ~any(Hint == Himct)
% 	    hInterpolated  = [Hint(Hint > Himct) Himct Hint(Hint < Himct)];
% 	    saInterpolated = [im2(im2 < imCt) imCt im2(im2 > imCt)];
% 	else
% 	    hInterpolated  = Hint;
% 	    saInterpolated = im2;
%     end
% 	% Calc PC
% 	PC_ = zeros(1,size(saInterpolated,2));
% 	for i = 1:size(saInterpolated,2)
% 		PC_(1,i) = logncdf_(saInterpolated(1,i),log(imCtM(ss)),betaC);
% 	end
% else
% 	saInterpolated = [];
% 	hInterpolated  = [];
% 	if ~any(hUser == Himct)
% 	    hUser  = [hUser(hUser > Himct) Himct hUser(hUser < Himct)];
% 	    saUser = [saUser(saUser < imCt) imCt saUser(saUser > imCt)];
% 	end
% 
% 	% Calc PC
% 	PC_ = zeros(1,size(saUser,2));
% 	for i = 1:size(saUser,2)
% 		PC_(1,i) = logncdf_(saUser(1,i),log(imCtM(ss)),betaC);
% 	end
% end



end

