function [wq,xq,yq] = tri_quadcofs (coord,flag)
%
d = size(coord,2);
if d ~= 2 return; end;
if ~ismember(flag,[1 3 7 13]) return;end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set coordinates and area
x1 = coord(1,1);x2 = coord(2,1);x3 = coord(3,1);
y1 = coord(1,2);y2 = coord(2,2);y3 = coord(3,2);

vol = x2*y3+x3*y1+x1*y2-x2*y1-x3*y2-x1*y3;  
vol = abs(vol)/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% one point 
if flag == 1 %% one point 

wq(1) = 1.0E+00;

  xq(1) = ( x1 + x2 + x3) /3;
  yq(1) = ( y1 + y2 + y3) /3;

elseif flag == 3 %% three points 

  wq(1) = 1.0E+00 / 3.0E+00;
  wq(2) = wq(1);
  wq(3) = wq(1);

  xq(1) = 0.5E+00 * ( x1 + x2 );
  xq(2) = 0.5E+00 * ( x2 + x3 );
  xq(3) = 0.5E+00 * ( x1 + x3 );

  yq(1) = 0.5E+00 * ( y1 + y2 );
  yq(2) = 0.5E+00 * ( y2 + y3 );
  yq(3) = 0.5E+00 * ( y1 + y3 );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% seven point (for polyn. of order 5)
elseif flag == 7 %%  7 point 

  AL1 = 0.0597158717;
BE1 = 0.4701420641;
AL2 = 0.7974269853;
BE2 = 0.1012865073;
WC1 = 0.2250000000;
WC2 = 0.1323941527;
WC3 = 0.1259391805;
ONETRI = 0.33333333;

  sp1 = [ONETRI,AL1,BE1,BE1,AL2,BE2,BE2];
  sp2 = [ONETRI,BE1,AL1,BE1,BE2,AL2,BE2];
  sp3 = [ONETRI,BE1,BE1,AL1,BE2,BE2,AL2];
  wq =[WC1,WC2,WC2,WC2,WC3,WC3,WC3];

  xq = x1*sp1 + x2*sp2 + x3*sp3;
  yq = y1*sp1 + y2*sp2 + y3*sp3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13 point
elseif flag == 13 

  for ( i = 1 : 3 )
    wq(i) = 0.175615257433204E+00;
    ii = i + 3;
    wq(ii) = 0.053347235608839E+00;
    ii = i + 6;
    iii = ii + 3;
    wq(ii) = 0.077113760890257E+00;
    wq(iii) = wq(ii);
  end

  wq(13) = -0.14957004446767E+00;

  z1 = 0.479308067841923E+00;
  z2 = 0.260345966079038E+00;
  z3 = 0.869739794195568E+00;
  z4 = 0.065130102902216E+00;
  z5 = 0.638444188569809E+00;
  z6 = 0.312865496004875E+00;
  z7 = 0.048690315425316E+00;

  xq( 1) = z1 * x1 + z2 * x2 + z2 * x3;
  yq( 1) = z1 * y1 + z2 * y2 + z2 * y3;
  xq( 2) = z2 * x1 + z1 * x2 + z2 * x3;
  yq( 2) = z2 * y1 + z1 * y2 + z2 * y3;
  xq( 3) = z2 * x1 + z2 * x2 + z1 * x3;
  yq( 3) = z2 * y1 + z2 * y2 + z1 * y3;
  xq( 4) = z3 * x1 + z4 * x2 + z4 * x3;
  yq( 4) = z3 * y1 + z4 * y2 + z4 * y3;
  xq( 5) = z4 * x1 + z3 * x2 + z4 * x3;
  yq( 5) = z4 * y1 + z3 * y2 + z4 * y3;
  xq( 6) = z4 * x1 + z4 * x2 + z3 * x3;
  yq( 6) = z4 * y1 + z4 * y2 + z3 * y3;
  xq( 7) = z5 * x1 + z6 * x2 + z7 * x3;
  yq( 7) = z5 * y1 + z6 * y2 + z7 * y3;
  xq( 8) = z5 * x1 + z7 * x2 + z6 * x3;
  yq( 8) = z5 * y1 + z7 * y2 + z6 * y3;
  xq( 9) = z6 * x1 + z5 * x2 + z7 * x3;
  yq( 9) = z6 * y1 + z5 * y2 + z7 * y3;
  xq(10) = z6 * x1 + z7 * x2 + z5 * x3;
  yq(10) = z6 * y1 + z7 * y2 + z5 * y3;
  xq(11) = z7 * x1 + z5 * x2 + z6 * x3;
  yq(11) = z7 * y1 + z5 * y2 + z6 * y3;
  xq(12) = z7 * x1 + z6 * x2 + z5 * x3;
  yq(12) = z7 * y1 + z6 * y2 + z5 * y3;
  xq(13) = ( x1 + x2 + x3 ) / 3.0E+00;
  yq(13) = ( y1 + y2 + y3 ) / 3.0E+00;


end;

wq = vol*wq;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


