%---------------------------------------------------------------------
%   PURPOSE               - NUMERICAL INTEGRATION OF A TABULATED COMPLEX 
%									 FUNCTION USING A TAYLOR'S EXPANSION OF
%                           THE FIRST  3N + 1  POINTS FOLLOWED BY A 
%                           TRAPEZOIDAL RULE INTEGRATION OF THE 
%                           REMAINING  NP - (3N + 1)  POINTS
%
%   USAGE                 - FUN = CTABINT (X, Y, NP)
%
%   ARGUMENTS    X        - VECTOR OF LENGTH  NP  CONTAINING THE VALUES 
%                           OF THE INDEPENDENT VARIABLE (STORED IN
%                           CONSECUTIVELY INCREASING OR DECREASIN ORDER)
%
%                Y        - VECTOR OF LENGTH  NP  CONTAINING THE VALUES
%                           OF THE DEPENDENT VARIABLE (ORDERED TO
%                           CORRESPOND TO THE INDEPENDENT VARIABLE) 
%
%                NP       - NUMBER OF THE PROVIDED POINTS 
%
%   REFERENCE             - LOS ALAMOS SCIENTIFIC LABORATORIES
%---------------------------------------------------------------------

function y=tabint (X,Y,NP)
      B = 0;
      N = NP - 1;
      K = MOD (N, 3);
      M = N - 2 - K;
      B=0;
      B1=0;
      for I = 1:3:M 
         H1 = X(I + 1) - X(I);
         H2 = X(I + 2) - X(I + 1);
         H3 = X(I + 3) - X(I + 2);
         H4 = H2 + H3;
         A4 = (H1 + H4)*(H4 * (H1 + H4)+(2.0 * H4 - H1)*(H3 - H1 - H2))/ (12.0 * H3 * H4); 
         A3 = (H1 + H4) ^ 3 * (H4 - H1) / (12.0 * H2 * H3 * (H1 + H2));
         A1 = ((H1 * H1 - H4 * H4) / 2.0 + H2 * A3 + H4 * A4) / H1;
         A2 = H1 + H4 - (A1 + A3 + A4); 
         B = B + A1 * Y(I) + A2 * Y(I + 1) + A3 * Y(I + 2) + A4 * Y(I + 3); 
      end    
          if (K==0)
             y=B;
             %B
          else
          M = NP - K;
      for J = M:1:N
          B1=B+(X(J + 1) - X(J)) * (Y(J + 1) + Y(J)) / 2.;
      end
      y=B1;
      %B1
          end