#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <assert.h>
#include "tipsydefs.h"  
/* needed by readdmstd */
#include <rpc/types.h>
#include <rpc/xdr.h>

static XDR xdrs; 
/**/

FILE *fp;

/* Prototypes */
void readdmtip_(int *ndark, int *nstar, int *nsph, float *xx, float *yy, float *zz,float *vx, float *vy, float *vz, float *mass, float *eps,
		float *phi,double *time, float *soft,float *metals,
                float *tform, float* rho,float* temp, float *pot ); 

void readdmstd_(int *ndark, int *nstar, int *nsph, float *xx, float *yy, float *zz, 
		float *vx, float *vy, float *vz, float *mass, float *eps,
		float *phi,double *time,float *soft,float *metals,
                float *tform, float* rho,float* temp, float *pot ); 



/* XDR handling functions */ 


int xdr_header()  
{
  int pad;
  
  if(xdr_double(&xdrs, &header.time) != TRUE)
    return 0; 
  if(xdr_int(&xdrs, &header.nbodies) != TRUE)
    return 0;
  if(xdr_int(&xdrs, &header.ndim) != TRUE)
    return 0;
  if(xdr_int(&xdrs, &header.nsph) != TRUE)
    return 0;
  if(xdr_int(&xdrs, &header.ndark) != TRUE)
    return 0;
  if(xdr_int(&xdrs, &header.nstar) != TRUE)
    return 0; 
  if(xdr_int(&xdrs, &pad) != TRUE)  
    return 0;  
  return 1;  
}  

int xdr_gas()  
{ 
  int i;  

  if(sizeof(Real) == sizeof(float)) 
    {  
      xdr_vector(&xdrs,(char *)gas_particles, 
	     header.nsph*(sizeof(*gas_particles)/sizeof(Real)), 
	     sizeof(Real),(xdrproc_t) xdr_float);   
    }   
  return 0;
}    

int xdr_dark() 
{
  int i;
  
  if(sizeof(Real) == sizeof(float)) 
    { 
      xdr_vector(&xdrs,(char *)dark_particles,
	     header.ndark*(sizeof(*dark_particles)/sizeof(Real)),
	     sizeof(Real),(xdrproc_t) xdr_float); 
    }
   return 0;
}  

int xdr_star()  
{
  int i;

  if(sizeof(Real) == sizeof(float))
    {
      xdr_vector(&xdrs,(char *)star_particles,
	     header.nstar*(sizeof(*star_particles)/sizeof(Real)),
	     sizeof(Real),(xdrproc_t) xdr_float);
    } 
   return 0;
}    



/* Function definitions */


   void readdmtip_(int *ndark, int *nstar, int *nsph, float *xx, float *yy, float *zz, 
		float *vx, float *vy, float *vz, float *mass, float *eps,
		float *phi,double *time, float *soft,float *metals,
                float *tform, float* rho,float* temp, float *pot )

 
{
  /*	 FILE *fp; */
	 struct dark_particle dp;
	 struct gas_particle gp;
         struct star_particle sp;  
	 struct dump head;
	 int i;
	 int j;

	 /*
	 ** reads header
	 */
	 printf(" read files\n");
	 fread(&head,sizeof(struct dump),1,fp);   
	 printf("  Particles in TIPSY file: %d gas, %d dark, %d stars\n",
	     head.nsph,head.ndark,head.nstar);
	 assert(head.ndim==3);
	 *time=head.time;
	 printf("  Time: %g\n",*time);

	 /*
	 ** reads   gas particles 
	 */ 
	 i= 0;                               
	 while (i<head.nsph) {
		 fread(&gp,sizeof(struct gas_particle),1,fp);
		 mass[i]=gp.mass;
		 xx[i]=gp.pos[0];
		 yy[i]=gp.pos[1];
		 zz[i]=gp.pos[2];
		 vx[i]=gp.vel[0];
		 vy[i]=gp.vel[1];
		 vz[i]=gp.vel[2];
		 rho[i]=gp.rho;
		 temp[i]=gp.temp;
		 soft[i]=gp.hsmooth;
		 metals[i]=gp.metals;
		 pot[i]=gp.phi;

		 ++i; 
	 }

	 /*
	 ** reads  dark particles
	 */
	 j = 0;
	 while (i<head.ndark) {
		 fread(&dp,sizeof(struct dark_particle),1,fp);
		 mass[i]=dp.mass;
		 xx[i] = dp.pos[0];
		 yy[i] = dp.pos[1];
		 zz[i] = dp.pos[2];
		 vx[i] = dp.vel[0];
		 vy[i] = dp.vel[1];
		 vz[i] = dp.vel[2];
		 eps[i] = dp.eps;
		 phi[i] = dp.phi;
		 ++i;
		 ++j;

	 }

/* reads  star particles */
   
	 j = 0;
	 while (j<head.nstar) {
	   fread(&sp,sizeof(struct star_particle),1,fp);
	   mass[i] = sp.mass;
	   xx[i] = sp.pos[0];
	   yy[i] = sp.pos[1];
	   zz[i] = sp.pos[2]; 
	   vx[i] = sp.vel[0];
	   vy[i] = sp.vel[1];
	   vz[i] = sp.vel[2];
	   metals[i]=sp.metals;
	   tform[i]=sp.tform;
	   soft[i]=sp.eps;
	   pot[i]=sp.phi; 
	   ++i; 
	   ++j;  
	 } 

	 
	 printf("  %d particles read from TIPSY file.\n",i);
	 	 *ndark = head.ndark;

 }



 void readdmstd_( int *ndark, int *nstar, int *nsph, float *xx, float *yy, float *zz, 
		float *vx, float *vy, float *vz, float *mass, float *eps,
		float *phi,double *time,float *soft,float *metals,
                float *tform, float* rho,float* temp, float *pot  ) 
{
  /*	 FILE *fp; */
  struct dark_particle dp;
  struct gas_particle gp;
  struct star_particle sp; 
  struct dump head;

  int i;
  int j;


  printf( " read files  1\n" );
  assert(sizeof(float)==sizeof(Real));

  
  /*
  **  read stream of data from file fp
  */
  xdrstdio_create(&xdrs, fp, XDR_DECODE);

  printf( " read files  2\n" );
  
  /*
  **  read Tipsy header
  */
  xdr_header(); 

  printf( " read files  3\n" );
  
  
  fprintf(stderr,"  Particles in TIPSY file: %d gas, %d dark, %d stars\n",
             header.nsph,header.ndark,header.nstar);
  assert(header.ndim==3);
  *time=header.time;
  printf("  Time: %g\n",*time);

  /*
  ** reads  any gas particles 
  */ 
  i=0;
  while (i<header.nsph) {
    /*    xdr_vector(&xdrs,(char *) &gp, (sizeof(gp)/sizeof(Real)),
	  sizeof(Real),(xdrproc_t) xdr_float);*/

    fread(&sp,sizeof(struct star_particle),1,fp);


    mass[i]=gp.mass; 
    xx[i]=gp.pos[0];
    yy[i]=gp.pos[1];
    zz[i]=gp.pos[2];
    vx[i]=gp.vel[0];
    vy[i]=gp.vel[1];
    vz[i]=gp.vel[2];
    rho[i]=gp.rho;
    temp[i]=gp.temp;
    soft[i]=gp.hsmooth;
    metals[i]=gp.metals;
    pot[i]=gp.phi; 

    ++i; 
  } 

  /*
  ** reads data for dark particles
  */
  j = 0;
  while (j<header.ndark) {
    xdr_vector(&xdrs,(char *) &dp, (sizeof(dp)/sizeof(Real)),
	       sizeof(Real),(xdrproc_t) xdr_float);
    mass[i]=dp.mass;
    xx[i] = dp.pos[0];
    yy[i] = dp.pos[1];
    zz[i] = dp.pos[2];
    vx[i] = dp.vel[0];
    vy[i] = dp.vel[1];
    vz[i] = dp.vel[2];
    eps[i] = dp.eps;
    phi[i] = dp.phi;

    /*    printf("  dark mass: %g\n",mass[i]);*/
    ++j;  
    ++i;} 

  j = 0;
  while (j<header.nstar) {

     xdr_vector(&xdrs,(char *) &sp, (sizeof(sp)/sizeof(Real)),
	       sizeof(Real),(xdrproc_t) xdr_float);

    mass[i] = sp.mass;
    xx[i] = sp.pos[0];
    yy[i] = sp.pos[1];
    zz[i] = sp.pos[2]; 
    vx[i] = sp.vel[0];
    vy[i] = sp.vel[1];
    vz[i] = sp.vel[2];
    metals[i]=sp.metals;
    tform[i]=sp.tform;
    soft[i]=sp.eps;
    pot[i]=sp.phi; 
    /*    printf("  stellar mass: %g\n",sp.mass); */
    ++i; 
    ++j;  
  }




  printf("  %d dark particles read from TIPSY file.\n",header.ndark);
     *ndark = header.ndark; 
  printf("  %d star particles read from TIPSY file.\n",header.nstar);
  *nstar = header.nstar;

   printf("  %d gas particles read from TIPSY file.\n",header.nsph); 
   *nsph = header.nsph; 

   
  xdr_destroy(&xdrs);
}