#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#define NMAX 10000
#define NNMAX 500

double * x = NULL;
double * x2 = NULL;
 double * y = NULL;
double * y2 = NULL;
 double * vx = NULL;
double * vx2 = NULL;
 double * vy = NULL;
double * vy2 = NULL;
 double * r = NULL;
double * r2 = NULL;

void initmem()
{
size_t sizeInBytes = NMAX * sizeof(double);
  x=(double *) malloc(sizeInBytes);
  vx=(double *) malloc(sizeInBytes);
  vy=(double *) malloc(sizeInBytes);
  y=(double *) malloc(sizeInBytes);
  r=(double *) malloc(sizeInBytes);
  x2=(double *) malloc(sizeInBytes);
  vx2=(double *) malloc(sizeInBytes);
  vy2=(double *) malloc(sizeInBytes);
  y2=(double *) malloc(sizeInBytes);
  r2=(double *) malloc(sizeInBytes);
}


void main (void)
{
  //double vx[NMAX],vy[NMAX],x[NMAX],y[NMAX],r[NMAX];
  // double vx2[NMAX],vy2[NMAX],x2[NMAX],y2[NMAX];
 double vmax,vmag,vsum,vtx,vty;
 double xdiff,ydiff,vxdiff,vydiff;
 double oneau,kpc,t;
 double dx1,dy1,dx2,dy2;
 int  i,n1,n2,ik,ic;

char fname[14];

/* File pointers for data output */
 FILE *fp1,*fpE;
 initmem();
oneau = 149597871.0E3;
/* one kilo parsec*/
kpc = 2.06e8*oneau;

 // File Name Counter 

 for (ik=1281;ik<2482;ik++) {

   t=(double) ik;

     sprintf(fname,"data/vdat_%04d",ik);

     printf("%s\n",fname);
     fp1=fopen(fname,"r");
     fscanf(fp1,"%d %d\n",&n1,&n2);
     fscanf(fp1,"%le %le %le %le\n",&dx1,&dy1,&dx2,&dy2);
     printf("n1= %d n2= %d dx1= %e dy1= %e dx2= %e dy2= %e\n",n1,n2,dx1,dy1,dx2,dy2);

     for (i=0;i<n1;i++)
       {
	 fscanf(fp1,"%le %le %le %le\n",&x[i],&y[i],&vx[i],&vy[i]);
	 r[i]=sqrt((x[i]-dx1)*(x[i]-dx1) + (y[i]-dy1)*(y[i]-dy1));
	 //	 printf("%e %e %e %e %e\n",x[i],y[i],vx[i],vy[i],r[i]);
       }

     for (i=n1;i<n1+n2;i++)
       {
	 fscanf(fp1,"%le %le %le %le\n",&x[i],&y[i],&vx[i],&vy[i]);
	 r[i]=sqrt((x[i]-dx2)*(x[i]-dx2) + (y[i]-dy2)*(y[i]-dy2));
	 //	 printf("%e %e %e %e %e\n",x[i],y[i],vx[i],vy[i],r[i]);
       }
     fclose(fp1);

     printf("============================\n");

     sprintf(fname,"data/masC_%04d",ik);
     printf("%s\n",fname);
     fp1=fopen(fname,"w");
     fprintf(fp1,"$DATA=CURVE2D\n");
     // fprintf(fp1,"%%toplabel= \"WITH Dark Matter\"\n");
     //     fprintf(fp1,"%%toplabel= \"NO Dark Matter\"\n");
     fprintf(fp1,"%%toplabel= \" \"\n");
     fprintf(fp1,"%%subtitle= \"Time = %4.2f (myr)\"\n",t);
     fprintf(fp1,"%%xmin= -400\n");
     fprintf(fp1,"%%xmax= +400\n");
     fprintf(fp1,"%%ymin= -400\n");
     fprintf(fp1,"%%ymax= +400\n");
     fprintf(fp1,"%%equalscale= T\n");
     fprintf(fp1,"%%xlabel= \"x (kpc)\"\n");
     fprintf(fp1,"%%ylabel= \"y (kpc)\"\n");
     fprintf(fp1,"%%linetype= 0\n");

     for (i=0;i<n1;i++)
       {
//compute tangential component
	 vtx=vx[i]*(y[i]-dy1)/r[i];
	 vty=-vy[i]*(x[i]-dx1)/r[i];
// Add up vtx and vty to get total tangential velocity
// Divide by 1000 to get km/
	 vmag=fabs(vtx + vty)/1000.0;
	 //	 printf("vmag= %f\n",vmag);

	 //	 printf ("vmag= %f\n",vmag);
//Set color depending upon radial velocity
	 if (vmag > 250.0) {
 	   ic=4;}
	   else if (vmag > 200.0) {
	     ic=3;}
	   else if (vmag > 100.0) {
	     ic=5;}
	   else {
	     ic=0;}

//Fixed color 
	 ic=7;

         fprintf(fp1,"%%markertype= 1\n");
         fprintf(fp1,"%%markercolor= %d\n",ic);
	 fprintf(fp1,"%f %f\n",x[i],y[i]);
         fprintf(fp1,"\n");
       }

     for (i=n1;i<n1+n2;i++)
       {
//compute tangential component
	 vtx=vx[i]*(y[i]-dy2)/r[i];
	 vty=-vy[i]*(x[i]-dx2)/r[i];
// Add up vtx and vty to get total tangential velocity
// Divide by 1000 to get km/
	 vmag=fabs(vtx + vty)/1000.0;
	 //	 printf("vmag= %f\n",vmag);

	 //	 printf ("vmag= %f\n",vmag);
//Set color depending upon radial velocity
	 if (vmag > 250.0) {
 	   ic=4;}
	   else if (vmag > 200.0) {
	     ic=3;}
	   else if (vmag > 100.0) {
	     ic=5;}
	   else {
	     ic=0;}
//Fixed color 
	 ic=3;

         fprintf(fp1,"%%markertype= 1\n");
         fprintf(fp1,"%%markercolor= %d\n",ic);
	 fprintf(fp1,"%f %f\n",x[i],y[i]);
         fprintf(fp1,"\n");
       }


     fprintf(fp1,"\n");
     fclose(fp1);

 //Next time step
 }

     //END OF PROGRAM
}