#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define pi 3.14159265358979323844
#define xmargin 10
#define ymargin 10
/* This is the number of frames you wish to process */
#define nframes 100
/* This are the dimensions of the box to draw in */
#define xboxsize 240
#define yboxsize (int)(xboxsize*0.75)
#define psscale xboxsize/512
/* Choose imoldy=0 for BUU and 1 form mol. dynamics */
#define imoldy 0
/* These times are subtracted from times in file to set clock */
#define timeoffsetbuu 1.0
#define timeoffsetmoldy 0.0
/* Choose next 3 lines appropriately */
#define title1 "BUU"
#define title2 "Kr + Nb"
#define title3 "200x sampling"
//#define title1 "Mol. dynamics"
//#define title2 "A=93 + A=93"
//#define title3 "with Coulomb"
#define radius0 (1.2/(1.0+imoldy))
#define realboundary (20.0/(1.0+imoldy))
#define timescale (1.0+15.0*imoldy)
/* This chooses the point for the perspective */
#define realzperspective realboundary*1.5
#define realxperspective realboundary*.2
#define realyperspective realboundary*(-.1)
/* These are euler angles.  Choose theta=90 for sideview, 0 for down pipe */
#define thetaview 25.0
#define phi1view 0.0
#define phi2view 0.0
/* This determines the no. of circles to overlay to make sphere, 3 is good */
#define nspheredraw 2
/* Choose at least max. no. of particles + 1 */
#define npartmax 187
void psheader(double xmin,double ymin,double xmax,double ymax);
void rectdraw(double xa,double ya,double xb,double yb);
void rectdrawfill(double xa,double ya,double xb,double yb);
void linedraw(double xa,double ya,double xb, double yb);
void textwrite(double xx, double yy,double fontsize,double angle,char *s);
void spheredraw(double xx,double yy,double radius,int it);
void xyfind(double xx,double yy,double zz,double *x0,double *y0);
void rgb(double red,double green,double blue);
void shell_sort(int n);
void euler(int ip);
double x[npartmax],y[npartmax],z[npartmax],rp[npartmax];
double xperspective,yperspective,zperspective;
int itype[npartmax];
FILE *psfile,*datafile;
int main(){
  char psfilename[40];
  char textstring[40];
  int ifile;
  int xmin,xmax,ymin,ymax;
  double xa,ya,xb,yb;
  double xaa,yaa,xbb,ybb,xc,yc,xx,yy,zz,cc,ss;
  double boundary,scale;
  double xtext,ytext,fontsize;
  double radius,rp0;
  double epera,time;
  int npart,i,imax,ipart,np,ip;
  datafile=fopen("xyz.dat","r");
  fscanf(datafile,"%lf",&epera);
  if(imoldy==1)epera=epera*8.0;
  for(ifile=1;ifile<=nframes;ifile++){
    sprintf(psfilename,"xyz.%03d.eps",ifile);
    psfile=fopen(psfilename,"w+");
    /* A page caries a size of 792 x 612, though we will go for 512 by 384 */
    /* for the area in which we draw */
    ymin=ymargin;
    ymax=ymin+yboxsize-1;
    xmin=xmargin;
    xmax=xmin+xboxsize-1;
    psheader(xmin,ymin,xmax,ymax);
    rgb(0.0,0.0,0.7);
    rectdraw(xmin,ymin,xmax,ymax);
    xb=xmax-30*psscale;
    yb=ymax-60*psscale;
    ya=ymin+30*psscale;
    boundary=0.5*(yb-ya);
    xa=xb-2.0*boundary;
    rgb(0.9,0.7,0.7);
    rectdrawfill(xa,ya,xb,yb);
    xc=0.5*(xa+xb);
    yc=0.5*(ya+yb);
    scale=boundary/realboundary;
    xperspective=realxperspective*scale;
    yperspective=realyperspective*scale;
    zperspective=realzperspective*scale;
    rp0=sqrt(pow(xperspective,2)+pow(yperspective,2)+pow(zperspective,2));
    fscanf(datafile,"%d %lf",&npart,&time);
    time=time*timescale;
    if(imoldy==0) time=time-timeoffsetbuu;
    if(imoldy==1) time=time-timeoffsetmoldy;
    printf("%d  t= %6.1f\n",ifile,time);
    xyfind(-boundary,-boundary,-2.0*boundary,&xaa,&yaa);
    xaa=xaa+xc;yaa=yaa+yc;
    xyfind(boundary,boundary,-2.0*boundary,&xbb,&ybb);
    xbb=xbb+xc;ybb=ybb+yc;
    rgb(0.8,0.6,0.6);
    rectdrawfill(xaa,yaa,xbb,ybb);
    rgb(0.0,0.0,0.7);
    rectdraw(xa,ya,xb,yb);
    linedraw(xa,ya,xaa,yaa);
    linedraw(xb,ya,xbb,yaa);
    linedraw(xb,yb,xbb,ybb);
    linedraw(xa,yb,xaa,ybb);
    rectdraw(xaa,yaa,xbb,ybb);
    rgb(0.0,0.0,0.7);
    fontsize=24.0*psscale;
    xtext=xmin+6;
    ytext=ymax-50;
    textwrite(xtext,ytext,fontsize,0.0,(char*)title1);
    ytext=ytext-1.5*fontsize;
    sprintf(textstring,"E/A = %4.1f MeV",epera);
    textwrite(xtext,ytext,fontsize,0.0,textstring);
    ytext=ytext-1.5*fontsize;
    textwrite(xtext,ytext,fontsize,0.0,(char*)title2);
    ytext=ytext-1.5*fontsize;
    textwrite(xtext,ytext,fontsize,0.0,(char*)title3);
    ytext=ya+10.0;
    xtext=xc-1.5*fontsize;
    sprintf(textstring,"t =%4.0f fm/c",time);
    textwrite(xtext,ytext,fontsize,0.0,textstring);
    ipart=0;
    ip=1;
    np=0;
    while(ipart<npart){
      ipart=ipart+1;
      fscanf(datafile,
	     "%lf %lf %lf %d",&x[ip],&y[ip],&z[ip],&itype[ip]);
      if(fabs(x[ip])<realboundary&&fabs(y[ip])<realboundary
	&&fabs(z[ip])<realboundary){
	ip=ip+1;np=np+1;
	euler(np);
	x[np]=scale*x[np];
	y[np]=scale*y[np];
	z[np]=scale*(z[np]-realboundary);
	rp[np]=sqrt(pow(x[np]-xperspective,2)+pow(y[np]-yperspective,2)
		    +pow(z[np]-zperspective,2));
      }
    }
    /* Sort the particles in inverse order according to rp */
    shell_sort(np);
    for(ipart=1;ipart<=np;ipart++){
      radius=scale*radius0*rp0/rp[ipart];
      xyfind(x[ipart],y[ipart],z[ipart],&xx,&yy);
      spheredraw(xc+xx,yc+yy,radius,itype[ipart]);
    }
    fprintf(psfile,"showpage \n");
    fclose(psfile);
  }
  return(0);
}
void psheader(double xmin,double ymin,double xmax,double ymax){
  fprintf(psfile,"%%!PS-Adobe-2.0 EPSF-2.0 \n");
  fprintf(psfile,
	  "%%%%BoundingBox: %5.1f %5.1f %5.1f %5.1f \n",xmin,ymin,xmax,ymax);
  fprintf(psfile,"%%%%Creator: S. Pratt, M.S.U. \n");
  fprintf(psfile,"%%%%EndComments \n");
}
void xyfind(double xx,double yy,double zz,double *x0,double *y0){
  double factor;
  factor=-zperspective/(-zperspective+zz);
  *x0=factor*xx+(1.0-factor)*xperspective;
  *y0=factor*yy+(1.0-factor)*yperspective;
}
void rgb(double red,double green,double blue){
  fprintf(psfile,"%4.2f %4.2f %4.2f setrgbcolor \n",red,green,blue);
}
void linedraw(double xa,double ya,double xb,double yb){
  fprintf(psfile,"%6.2f %6.2f moveto %6.2f %6.2f lineto stroke \n",xa,ya,xb,yb);
}
void rectdrawfill(double xa,double ya,double xb,double yb){
  fprintf(psfile,"newpath \n%6.2f %6.2f moveto %6.2f %6.2f lineto \n",
	  xa,ya,xb,ya);
fprintf(psfile,"%6.2f %6.2f lineto %6.2f %6.2f lineto \nclosepath fill\n",
xb,yb,xa,yb);
}
void rectdraw(double xa,double ya,double xb,double yb){
  fprintf(psfile,"newpath %6.2f %6.2f moveto %6.2f %6.2f lineto \n",
	  xa,ya,xb,ya);
fprintf(psfile,"%6.2f %6.2f lineto %6.2f %6.2f lineto \nclosepath stroke \n",
xb,yb,xa,yb);
}
void spheredraw(double xx,double yy,double radius,int it){
  int i;
  double theta0,thetaf,graylevel,graylevel0;
  double r,r2,yx,xci,yci;
  theta0=0.0;
  thetaf=360.0;
  for(i=1;i<=nspheredraw;i++){
    if(it==1)graylevel0=.9;
    if(it==2)graylevel0=.7;
    if(it==3)graylevel0=.8;
    graylevel=graylevel0;
    if(nspheredraw>1)
      graylevel=graylevel-graylevel0*0.3*(double)(i-1)/(double)(nspheredraw-1);
    r2=pow(radius,2)*(1.0-(double)(i-1)/(double)(nspheredraw));
    r=sqrt(r2);
    fprintf(psfile,"newpath \n");
    xci=xx+((double)(i-1)/(double)(nspheredraw))*r*0.15;
    yci=yy-((double)(i-1)/(double)(nspheredraw))*r*0.1;
    fprintf(psfile," %6.2f %6.2f %6.2f %6.2f %6.2f arc closepath \n",
	    xci,yci,r,theta0,thetaf);
    if(it==1) rgb(graylevel,0.3*graylevel,0.3*graylevel);
    if(it==2) rgb(0.3*graylevel,graylevel,0.3*graylevel);
    if(it==3) rgb(0.3*graylevel,graylevel,0.3*graylevel);
    fprintf(psfile," fill stroke \n");
  }
}
void textwrite(double xtext,double ytext, double fontsize,double angle,char *textstring){
  int n;
  fprintf(psfile,"/Helvetica findfont %5.1f scalefont setfont ",fontsize);
  fprintf(psfile,"%6.2f %6.2f moveto \n",xtext,ytext);
  if(fabs(angle)>0.0001) fprintf(psfile,"%6.1f rotate ",angle);
  fprintf(psfile,"(%s)show \n",textstring);
  if(fabs(angle)>0.0001){
    n=(int)angle/360.0;
    angle=angle-360.0*n;
    fprintf(psfile,"%6.1f rotate (%s)show %6.1f rotate \n",
	    angle,textstring,-angle);
  }
  else{
    fprintf(psfile,"(%s)show \n",textstring);
  }
}
/* This the Shell method which works very well for n <= 50 */
/* We use it for sorting in reverse order (see below) */
void shell_sort(int n){
  int i,j,inc,ve;
  double va,vb,vc,vd;
  inc=1;
  do {
    inc *=3;
    inc++;
  } while (inc <= n);
  do {
    inc /= 3;
    for (i=inc+1;i<=n;i++) {
      va=rp[i];
      vb=x[i];
      vc=y[i];
      vd=z[i];
      ve=itype[i];
      j=i;
      /* while (x[j-inc]>va) {*/
      /* For sorting in reverse order use below statement */
      while (rp[j-inc]<va) { 
	rp[j]=rp[j-inc];
	x[j]=x[j-inc];
	y[j]=y[j-inc];
	z[j]=z[j-inc];
	itype[j]=itype[j-inc];
	j-=inc;
	if(j<=inc) break;
      }
      rp[j]=va;
      x[j]=vb;
      y[j]=vc;
      z[j]=vd;
      itype[j]=ve;
    }
  } while (inc>1);
}
void euler(int np){
  static double cc1,ss1,cc,ss,cc2,ss2;
  static int entry=0;
  double xx,yy,zz;
  entry=entry+1;
  if(entry==1){
    cc1=cos(phi1view*pi/180.0);ss1=sin(phi1view*pi/180.0);
    cc=cos(thetaview*pi/180.0);ss=sin(thetaview*pi/180.0);
    cc2=cos(phi2view*pi/180.0);ss2=sin(phi2view*pi/180.0);
  }
  xx=x[np]*cc1+y[np]*ss1;
  y[np]=y[np]*cc1-x[np]*ss1;
  x[np]=xx;
  zz=z[np]*cc+x[np]*ss;
  x[np]=x[np]*cc-z[np]*ss;
  z[np]=zz;
  xx=x[np]*cc2+y[np]*ss2;
  y[np]=y[np]*cc2-x[np]*ss2;
  x[np]=xx;
}