/* Autor: Amaro R. Silva, Prof. IST- Abr 1994
* ponha libML.a e mathlink.h em locais acessiveis ao compilador e linker
* uso geral em unix: 
*
* mprep templatefile.tm -o templatefile.tm.c  
* gcc -o exefile  codefile.c  templatefile.tm.c  -lML  -lm
*
*  (veja as notas de ajuda para compilar noutros sistemas)
*
* No Mathematica designe lp=Install[.../arrow3D] porque vai precisar de
* escrever nesse link, para que o default Viewpoint seja usado nos cálculos.  
*/

#include "mathlink.h"
#include "mdefs.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define   pi   		4*atan(1.) 
#define	  YES		1
#define   NO		0

/*==========================Global Vars============================*/

double  *arrowlist,*VP,INIT=YES;

	
arrow3D(double *vect, int n, double *point, int m, double s, double h)
{
	int	i,j,k=0;
	int count;
	double	nv[3];
	double	uv[3];
	double  x,y,z,L,l;
	double  tip[] = {-.2, 0, -.35, 2., 0., 0., -.35, -1.};
	
	if(INIT==YES){
		MLPutFunction(stdlink,"Evaluate",1);
		MLPutFunction(stdlink,"ToExpression",1);
		MLPutString(stdlink,
		"LinkWrite[lp,(Cases[Options[Graphics3D],z_/;!FreeQ[z,ViewPoint]])[[1,2]]];LinkRead[lp]");
		MLEndPacket(stdlink);
		MLFlush(stdlink);
		MLGetRealList(stdlink,&VP,&count);
		MLNewPacket(stdlink);
		MLClearError(stdlink);
		MLFlush(stdlink);
	INIT=NO;}
	
	{x=vect[0];  y=vect[1];  z=vect[2];}
	
	L=sqrt(x*x + y*y + z*z+.01);	
	
	{nv[0]=x/L;  nv[1]=y/L;  nv[2]=z/L;}
	
	if(s==0 ) {s=1/L;}   
	
	{vect[0]=s*x;   vect[1]=s*y;  vect[2]=s*z;  }
	
	l=.05*h/sqrt(VP[1]*VP[1] + VP[2]*VP[2] + VP[3]*VP[3]);
	
	uv[0] = l*(nv[2]*VP[1] - nv[1]*VP[2]); 
	uv[1] = l*(nv[0]*VP[2] - nv[2]*VP[0]);
	uv[2] = l*(nv[1]*VP[0] - nv[0]*VP[1]);
	
	MLFlush(stdlink);	
	MLClearError(stdlink);
	
	MLPutFunction(stdlink,"List",2);
	MLPutFunction(stdlink,"Line",1);
	MLPutFunction(stdlink,"List",2);	
	MLPutRealList(stdlink,point,m);
		k = 0;
		while (k <= n) {
			arrowlist[k] =point[k] + vect[k] + tip[0] * nv[k];
			k++;
		}
	MLPutRealList(stdlink,arrowlist,n);
	
	MLPutFunction(stdlink,"Polygon",1);
	MLPutFunction(stdlink,"List",4);
	
	for(i=0;i<8;i+=2){
		k = 0;
		while (k <= n) {
			arrowlist[k] = point[k] + vect[k] + tip[i] * nv[k] + tip[i + 1] * uv[k];
			k++;
		}
		MLPutRealList(stdlink,arrowlist,3);
	}

	MLEndPacket(stdlink);
	MLFlush(stdlink);
	MLClearError(stdlink);
	
	return;
	}


/*==========================main============================*/

int   main(int argc, char *argv[])
{

arrowlist=(double *)malloc(sizeof(double)*255);

MLMain(argc, argv);

free(arrowlist);

return ; 
}