function y = thr(p,a,b,n,ind)

% Variables:
% p - period
% a - parameter (a) of uniform(a,b)
% b - parameter (b) of uniform(a,b)
% n - number of flows
% ind - lower or upper bound

if n<1
%impossible
 y=-1;
 y_lb=y;
 y_ub=y;

elseif b<a 
% impossible
 y=-1;
 y_lb=y;
 y_ub=y;

elseif a==b
% deterministic scenario
 if p >= a
  y=(p-a)/p;
 else
  y=(ceil(a/p)*p-a)/p;  
 end
 y_lb=y;
 y_ub=y;

else
% b>a
 if p>=b
  y=(p-(a+(b-a)/(n+1)))/p;
  y_lb=y;
  y_ub=y;

 elseif p>=a
% p>=a and p<b
  k=floor(b/p);
  p1=(((p-a)/(b-a))^n)*(p-(a+(p-a)/(n+1)))/p;
  p2_lb=0;
  p2_ub=0;
  if (k>=2)
   for i=1:(k-1)
     p2_lb=p2_lb+(((((i+1)*p-a)/(b-a))^n)-(((i*p-a)/(b-a))^n))*(((i+1)*p-(i*p+((i+1)*p-i*p)/(n+1)))/((i+1)*p));
     p2_ub=p2_ub+(((((i+1)*p-a)/(b-a))^n)-(((i*p-a)/(b-a))^n))*(((i+1)*p-(a+((i+1)*p-a)/(n+1)))/((i+1)*p));
   end
  end
  %p3_lb=p3_lb*(((k+1)*p-k*p-((k+1)*p-kp)/(n+1))/((k+1)*p))
  p3_lb=(1-((k*p-a)/(b-a))^n)*((p-p/(n+1))/((k+1)*p));
  p3_ub=(1-((k*p-a)/(b-a))^n)*(((k+1)*p-a-((k+1)*p-a)/(n+1))/((k+1)*p));
  y_lb=p1+p2_lb+p3_lb;
  y_ub=p1+p2_ub+p3_ub;

 else
 %p<a
   k=floor(b/p);
   if k==1
    y=(ceil(a/p)*p-(a+(b-a)/(n+1)))/(ceil(a/p)*p);
    y_lb=y;
    y_ub=y;

   else
   %k>=2
    p1=(((ceil(a/p)*p-a)/(b-a))^n)*((ceil(a/p)*p-a-(ceil(a/p)*p-a)/(n+1))/(ceil(a/p)*p));
    p2_lb=0;
    p2_ub=0;
    if (k>=2)
     %ovo ovde dole ne valja -> ne sme da ide od jedan
     for i=ceil(a/p):(k-1)
      p2_lb=p2_lb+(((((i+1)*p-a)/(b-a))^n)-(((i*p-a)/(b-a))^n))*(((i+1)*p-(i*p+((i+1)*p-i*p)/(n+1)))/((i+1)*p));
      p2_ub=p2_ub+(((((i+1)*p-a)/(b-a))^n)-(((i*p-a)/(b-a))^n))*(((i+1)*p-(a+((i+1)*p-a)/(n+1)))/((i+1)*p));
     end
    end
    %p3_lb=(1-((k*p-a)/(b-a))^n)*(((k+1)*p-k*p-((k+1)*p-kp)/(n+1))/((k+1)*p));
    p3_lb=(1-((k*p-a)/(b-a))^n)*((p-p/(n+1))/((k+1)*p));
    p3_ub=(1-((k*p-a)/(b-a))^n)*(((k+1)*p-a-((k+1)*p-a)/(n+1))/((k+1)*p));
    y_lb=p1+p2_lb+p3_lb;
    y_ub=p1+p2_ub+p3_ub;
   end
 end
end


if ind==0
 y=y_lb;
else
 y=y_ub;
end