#include <stdio.h>
#include <strings.h>

#define  SUN
#define  U_CHAR   1
#define  U_INT    2
#define  U_FLOAT  3

char Revision[]={'1','.','1','1'};

main(argc, argv)
int argc;
char *argv[];
{
  char Rev[5];
  int i, j, nsp, nnp, np[10];
  float ftemp;
  short temp[10], temp1;
  FILE *DMPFile;

  float **x, **vx, **vy, **vz, t, sigma;
  float extq, extq_1, extq_2, extq_3, jwall[10];
   
  /*********************************************************/

  if(argc < 2)
  {
    puts("Must give 2 args.");
    exit(1);
  }
  if ((DMPFile = fopen(argv[1], "r+b")) == NULL)
  {
    puts("Conver: open failed on the first argument.");
    exit(1);
  }
  fread(Rev, sizeof(char), 4, DMPFile);
  for (i=0; i<4; i++)
    if (Rev[i]!=Revision[i])
    {
      puts("Incompatible dump file version");
      putchar(7);
      exit(1);
    }
  read_data(&t, 4, 1, DMPFile, U_FLOAT);
  read_data(&sigma, 4, 1, DMPFile, U_FLOAT);
  read_data(&extq, 4, 1, DMPFile, U_FLOAT);
  read_data(&extq_1, 4, 1, DMPFile, U_FLOAT);
  read_data(&extq_2, 4, 1, DMPFile, U_FLOAT);
  read_data(&extq_3, 4, 1, DMPFile, U_FLOAT);
  read_data(&nsp, 4, 1, DMPFile, U_INT);
  read_data(jwall, 4, nsp, DMPFile, U_FLOAT);
  read_data(np, 4, nsp, DMPFile, U_INT);

 /*********************************************************/
  
  x  = (float **)malloc(nsp*sizeof(float *));
  vx = (float **)malloc(nsp*sizeof(float *));
  vy = (float **)malloc(nsp*sizeof(float *));
  vz = (float **)malloc(nsp*sizeof(float *));
  
  for(i=0; i<nsp; i++)
  {
    x[i]  = (float *)malloc(2*np[i]*sizeof(float));
    vx[i] = (float *)malloc(2*np[i]*sizeof(float));
    vy[i] = (float *)malloc(2*np[i]*sizeof(float));
    vz[i] = (float *)malloc(2*np[i]*sizeof(float));
  }
  
  if (!x[nsp-1] || !vx[nsp-1] || !vy[nsp-1] || !vz[nsp-1])
  {
    puts("Convert: x, v malloc failed");
    exit(1);
  }

 /*********************************************************/

  for (j=0; j<nsp; j++)
  {
    for (i=0; i<np[j]; i++)
    {
      read_data(&ftemp, 4, 1, DMPFile, U_FLOAT);
      x[j][i] = ftemp;
      read_data(&ftemp, 4, 1, DMPFile, U_FLOAT);
      vx[j][i]= ftemp;
      read_data(&ftemp, 4, 1, DMPFile, U_FLOAT);
      vy[j][i]= ftemp;
      read_data(&ftemp, 4, 1, DMPFile, U_FLOAT);
      vz[j][i]= ftemp;
    }
  }
  fclose(DMPFile);

  /*******************************************************************/
  /* Doubling the number of particles by copying the list over again */
  
  for (j=0; j<nsp; j++)
  {
    for (i=0; i<np[j]; i++)
    {
      x[j][i+np[j]] = x[j][i];
      vx[j][i+np[j]]= vx[j][i];
      vy[j][i+np[j]]= vy[j][i];
      vz[j][i+np[j]]= vz[j][i];
    }
    np[j]*=2;
  }

  /*********************************************************/
  
  if ((DMPFile = fopen(argv[2], "w")) == NULL)
  {
    puts("Conver: open failed on the second argument.");
    return;
  }
  
  write_data(Revision, 1, 4, DMPFile, U_CHAR);
  write_data(&t, 4, 1, DMPFile, U_FLOAT);
  write_data(&sigma, 4, 1, DMPFile, U_FLOAT);
  write_data(&extq, 4, 1, DMPFile, U_FLOAT);
  write_data(&extq_1, 4, 1, DMPFile, U_FLOAT);
  write_data(&extq_2, 4, 1, DMPFile, U_FLOAT);
  write_data(&extq_3, 4, 1, DMPFile, U_FLOAT);
  write_data(&nsp, 4, 1, DMPFile, U_INT);
  write_data(jwall, 4, nsp, DMPFile, U_FLOAT);
  write_data(np, 4, nsp, DMPFile, U_INT);

  for (j=0; j<nsp; j++)
  {
    for (i=0; i<np[j]; i++)
    {
      ftemp = x[j][i];
      write_data(&ftemp, 4, 1, DMPFile, U_FLOAT);
      ftemp = vx[j][i];
      write_data(&ftemp, 4, 1, DMPFile, U_FLOAT);
      ftemp = vy[j][i];
      write_data(&ftemp, 4, 1, DMPFile, U_FLOAT);
      ftemp = vz[j][i];
      write_data(&ftemp, 4, 1, DMPFile, U_FLOAT);
    }
  }
  fclose(DMPFile);
}

/****************************************************************/


read_data(ptr, size, nitems, stream, type)
char *ptr;
FILE *stream;
int size, nitems, type;
{
#ifdef IBM
  sun_read(ptr, size, nitems, stream);
#endif

#ifdef DEC
  fread(ptr, size, nitems, stream);
#endif

#ifdef SUN
  sun_read(ptr, size, nitems, stream);
#endif

#ifdef UNICOS
  switch(type)
    {
    case U_FLOAT:
      unicos_read_float(ptr, nitems, stream);
      break;
      
    case U_INT:
      unicos_read_int(ptr, nitems, stream);
      break;
      
    case U_CHAR:
      fread(ptr, size, nitems, stream);
      break;
    }
#endif
}

/****************************************************************/

write_data(ptr, size, nitems, stream, type)
char *ptr;
FILE *stream;
int size, nitems, type;
{
#ifdef IBM
  sun_write(ptr, size, nitems, stream);
#endif
  
#ifdef DEC
  fwrite(ptr, size, nitems, stream);
#endif
  
#ifdef SUN
  sun_write(ptr, size, nitems, stream);
#endif
  
#ifdef UNICOS
  switch(type)
    {
    case U_FLOAT:
      unicos_write_float(ptr, nitems, stream);
      break;
      
    case U_INT:
      unicos_write_int(ptr, nitems, stream);
      break;
      
    case U_CHAR:
      fwrite(ptr, size, nitems, stream);
      break;
    }
#endif
}

/****************************************************************/

sun_read(ptr, size, nitems, stream)
char *ptr;
int size;
int nitems;
FILE *stream;
{
  int i, j;
  
  for(j= 0; j <nitems; j++)
    for(i= size-1; i >=0; i--)
      fread(ptr+j*size+i, 1, 1, stream);
}

/****************************************************************/

sun_write(ptr, size, nitems, stream)
char *ptr;
int size;
int nitems;
FILE *stream;
{
  int i, j;
  
  for(j= 0; j <nitems; j++)
    for(i= size-1; i >=0; i--)
      fwrite(ptr+j*size+i, 1, 1, stream);
}

/****************************************************************/

#ifdef UNICOS
unicos_read_float(ptr, nitems, stream)
float *ptr;
FILE *stream;
int nitems;
{
  unsigned long *iptr;
  unsigned char a,b,c,d,t;
  int j;

  for(j= 0; j <nitems; j++)
  {
    iptr = (unsigned long *)(ptr +j);
    
    fread(&d,1,1,stream);
    fread(&c,1,1,stream);
    fread(&b,1,1,stream);
    fread(&a,1,1,stream);

    if( a == 0 && b == 0 && c == 0 && d == 0) *iptr = 0;
    else 
    {
      t = ((a << 1) | (b >> 7));
      *iptr = ((a  >> 7) << 63) | ((t -126 + 040000) << 48);
      *iptr |= (b | 0x80) << 40;  
      *iptr |= c << 32;
      *iptr |= d << 24;
    }
  }
}
#endif

/****************************************************************/

#ifdef UNICOS
unicos_read_int(ptr, nitems, stream)
int *ptr;
FILE *stream;
int nitems;
{
  unsigned long *iptr;
  unsigned char a,b,c,d,t;
  int j;

  for(j= 0; j <nitems; j++)
  {
    iptr = (unsigned long *)(ptr +j);
    
    fread(&a,1,1,stream);
    fread(&b,1,1,stream);
    fread(&c,1,1,stream);
    fread(&d,1,1,stream);

    *iptr = d << 24 | c << 16 | b << 8 | a;
    *iptr |= (d >> 7) ? 0xffffffff00000000 : 0x0;
  }
}
#endif

/****************************************************************/

#ifdef UNICOS
unicos_write_float(ptr, nitems, stream)
float *ptr;
FILE *stream;
int nitems;
{
  unsigned long *iptr;
  unsigned char a,b,c,d,t;
  int j;
  
  for(j= 0; j <nitems; j++)
  {
    iptr = (unsigned long *)(ptr +j);

    t = ((*iptr >> 48) & 0xff) - 040000 + 126;
    
    a = ((*iptr & (1 << 63)) ? 1 << 7 : 0) | t >> 1; 
    b = (t & 0x1) << 7 | ((*iptr >> 40) & 0x7f);
    c = *iptr >> 32;
    d = *iptr >> 24;
    
    fwrite(&d,1,1,stream);
    fwrite(&c,1,1,stream);
    fwrite(&b,1,1,stream);
    fwrite(&a,1,1,stream);
  }
}
#endif

/****************************************************************/

#ifdef UNICOS
unicos_write_int(ptr, nitems, stream)
int *ptr;
FILE *stream;
int nitems;
{
  unsigned long *iptr;
  unsigned char a,b,c,d,t;
  int j;
  
  for(j= 0; j <nitems; j++)
  {
    iptr = (unsigned long *)(ptr +j);
    
    a = *iptr & 0xff;
    b = *iptr >> 8 & 0xff;
    c = *iptr >> 16 & 0xff;
    d = *iptr >> 24 & 0xff;
    fwrite(&a,1,1,stream);
    fwrite(&b,1,1,stream);
    fwrite(&c,1,1,stream);
    fwrite(&d,1,1,stream);
  }
}
#endif

/****************************************************************/
/*  Heap sort  */

sort(n, ra, rb, rc, rd)
int n;
float *ra, *rb, *rc, *rd;
{
  int l, j, ir, i;
  float rra, rrb, rrc, rrd;

  l = (n >>1) +1;
  ir= n;
  
  for(;;)
  {
    if(l >1)
    {
      rra = ra[--l];
      rrb = rb[l];
      rrc = rc[l];
      rrd = rd[l];
    }
    else
    {
      rra   = ra[ir];
      rrb   = rb[ir];
      rrc   = rc[ir];
      rrd   = rd[ir];

      ra[ir]= ra[1];
      rb[ir]= rb[1];
      rc[ir]= rc[1];
      rd[ir]= rd[1];

      if(--ir == 1)
      {
	ra[1] = rra;
	rb[1] = rrb;
	rc[1] = rrc;
	rd[1] = rrd;
	return;
      }
    }
    i=l;
    j= l << 1;
    while(j <= ir)
    {
      if(j < ir && ra[j] < ra[j+1]) ++j;
      if(rra <ra[j])
      {
	ra[i] = ra[j];
	rb[i] = rb[j];
	rc[i] = rc[j];
	rd[i] = rd[j];
	j += (i=j);
      }
      else j=ir+1;
    }
    ra[i]= rra;
    rb[i]= rrb;
    rc[i]= rrc;
    rd[i]= rrd;
  }
}

/****************************************************************/