function y = smooth_indicator(x,x1,x2,a,mthd)
   // Smooth indicator function of an interval
   // Computes a smooth, analytic version of the 
   // indicator function of an interval [x1, x2]
   // (which is 1 on the interval and 0 otherwise).
   // The transition from 0 to 1 at x1 and from 
   // 1 to 0 at x2 is smooth and indefinetly 
   // derivable and is fastest for smaller a.
   //
   // Usage:
   //        y = smooth_indicator(x,x1,x2,a,mthd)
   // where  
   //        x:    scalar, vector or array
   //        x1:   lower end of interval 
   //        x2:   upper end of interval
   //        a:    smoothness parameter
   //        mthd: method ("atan" or "tanh")
   //
   // The function implements one of the expressions
   //
   // y = atan(%pi*a*(x-x1))/%pi - atan(%pi*a*(x-x2))/%pi;        
   // y = tanh(2*a*(x-x1))/2 - tanh(2*a*(x-x2))/2;
   // 
   // according to argument mthd
   //  
   // NOTES: The default value of a is 1000/(x2-x1),
   //        which yields a rather steep slope at
   //        the intervals ends. Any of x1 or x2 can
   //        be -%inf or %inf. If x1 > x2 the result
   //        is negative. If a = 0, it is equal to 
   //        the characteristic function, but if 
   //        x = x1 or x = x2, a NaN is obtained
   //           
   // -------------------------
   // Author: Federico Miyara
   // Date:   2022-08-28
   //         2022-09-22
   
   rhs = argn(2)
   if rhs<5 
      mthd = "tanh";
   end
   if rhs<4 
      a = 1000/(x2-x1);
   end
   if rhs==2
      x2 = x1 + 1; 
   end
   if rhs==1 
      x1 = 0;
      x2 = 1;
   end
   
   select mthd
   case "atan"
      y = atan(%pi*a*(x-x1))/%pi - atan(%pi*a*(x-x2))/%pi;
   case "tanh"
      y = tanh(2*a*(x-x1))/2 - tanh(2*a*(x-x2))/2;   
   end
   
   
endfunction
