/* phi_tri_pts.c  generate points uniformly in a triangle       */
/*                input, np, is number of points+1 along a side */
/*                triangle x1,y1 x2,y2 x3,y3                    */
/*                return is the number of points in xs,ys       */

#include <stdio.h>
#include <stdlib.h>

int phi_tri_pts(int np, double x1, double y1, double x2, double y2,
	        double x3, double y3, double xs[], double ys[])
{
  double x1a, y1a, x2a, y2a, x3a, y3a, x12, y12, x32, y32, nf;
  int p, j, nv;

  nv=0;
  nf=(double)np;
  xs[nv]=x1a=((2*np-2)*x1+x2+x3)/(2.0*nf); /* move in three points */
  ys[nv]=y1a=((2*np-2)*y1+y2+y3)/(2.0*nf);
  nv++;
  xs[nv]=x2a=((2*np-2)*x2+x3+x1)/(2.0*nf);
  ys[nv]=y2a=((2*np-2)*y2+y3+y1)/(2.0*nf);
  nv++;
  xs[nv]=x3a=((2*np-2)*x3+x1+x2)/(2.0*nf);
  ys[nv]=y3a=((2*np-2)*y3+y1+y2)/(2.0*nf);
  nv++;

  for(p=1; p<=np-1; p++)     /* 3 edges, n points each */
  { 
    xs[nv]=(p/nf)*x1a+((np-p)/nf)*x2a;
    ys[nv]=(p/nf)*y1a+((np-p)/nf)*y2a;
    nv++;
    xs[nv]=(p/nf)*x2a+((np-p)/nf)*x3a;
    ys[nv]=(p/nf)*y2a+((np-p)/nf)*y3a;
    nv++;
    xs[nv]=(p/nf)*x3a+((np-p)/nf)*x1a;
    ys[nv]=(p/nf)*y3a+((np-p)/nf)*y1a;
    nv++;
  }

  for(j=1; j<=np-2; j++)
  {
    x12=(j/nf)*x2a+((np-j)/nf)*x1a;
    y12=(j/nf)*y2a+((np-j)/nf)*y1a;
    x32=((np-j)/nf)*x3a+(j/nf)*x2a;
    y32=((np-j)/nf)*y3a+(j/nf)*y2a;
    for(p=1; p<=np-1-j; p++)
    {    
      xs[nv]=(p/(nf-j))*x12+((np-p-j)/(nf-j))*x32;
      ys[nv]=(p/(nf-j))*y12+((np-p-j)/(nf-j))*y32;
      nv++;
    }
  }
  return nv;
} /* end phi_tri_pts */