/*
###############################################################################
#                                                                             #
#                 Dust V1.04 - Copyright ©1994 by A.Maschke                   #
#                           All rights reserved.                              #
#-----------------------------------------------------------------------------#
#                                                                             #
#              Tutorial 3: Create a particle-explosion using Dust             #
#                            -changes only particle-positions                 #
#                            -writes scriptfile to stdout                     #
#                                                                             #
# OPERATION: 1. write a script which creates a particle-object and outputs    #
#               the PPOS-array and OCOUNT, e.g.:                              #
#                "load(1,objects/s1)                                          #
#                 load(2,objects/c1)                                          #
#                 o2p(1,2,1,p)                                                #
#                 savep(1,PExample.obj)                                       #
#                 getocount(1)                                                #
#                 !copy T:Dust.output PExample.oCount                         #
#                 getppos(1)                                                  #
#                 !copy T:Dust.output PExample.PPOS                           #
#                 !delete T:Dust.output                                       #
#                 exit"                                                       #
#             2. execute this script using Dust                               #
#             3. execute this programm writing "PExample >pex.bat"            #
#             4. execute the created file "pex.bat" using Dust"               #
#             5. That's it (You got 12 particle objects).                     #
#                                                                             #
#-----------------------------------------------------------------------------#
#                                                                             #
#      Language: ANSI-C                                                       #
#      Compiler: MCPP                                                         #
# Last modified: 25 September 1994                                            #
#                                                                             #
#      Problems: GNU-C                                                        #
#                Other float-format ? Please help !                           #
#                                                                             #
###############################################################################
*/
/*
*******************************************************************************
* NOTE: This code is Public Domain. You can do what you want with it.         *
*******************************************************************************
*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

const ZERO=0.0001;
struct Vector{float x;float y;float z;};

int oCount;  /* number of particles (to be read) */
struct Vector* pos; /* particle positions (to be read) */
float bb;           /* needed by RND2() */
int i,j; 
float centreX,centreY,centreZ; /* centre of all particles */

/*
*******************************************************************************
* returns a factor between 0.9 and 1.0                                        *
*******************************************************************************
*/
float RND2()
{
 bb=(float) rand();
 bb=bb/(float) RAND_MAX*900.0;
 bb=bb+100.0;
 bb=bb/10000.0;
 if (rand()>16384) {return 1.0+bb;} else {return 1.0-bb;}
}
/*
*******************************************************************************
* generates a filename like "obj.0235"                                        * 
*******************************************************************************
*/
void GenFn(char* str1,int tt)
{
 printf("%s.",str1);
 if (tt<1000) printf("0");
 if (tt<100) printf("0");
 if (tt<10) printf("0");
 printf("%d",tt);
}
/*
*******************************************************************************
* Read (x,y,z)-array from EXFILE; this example: pos                           *
*******************************************************************************
*/
int ReadReal(char* fn,void* addr,int count)
{
 FILE* fh;
 long len;
 fh=fopen(fn,"r");
 if (fh==NULL) return 1;
 len=fread(addr,sizeof(struct Vector),count,fh);
 fclose(fh);
 if (len!=count) return 1;
 return 0;
}
/*
*******************************************************************************
* Read a integer from EXFILE; this example: oCount                            *
*******************************************************************************
*/
int ReadInt(char* fn,void* addr)
{
 FILE* fh;
 long len;
 fh=fopen(fn,"rb");
 if (fh==NULL) return 1;
 len=fread(addr,2,1,fh);
 fclose(fh);
 if (len!=1) {return 1;} else {return 0;}
}
/*
*******************************************************************************
* Call the SETPxxx-Procedure of Dust to change the vector vp^[ind];           *
* this example: SETPPOS(1,ind,pos^[ind].x,pos^[ind].x,pos^[ind].z)            *
*******************************************************************************
*/
void WritePXXX(float x,float y,float z,int ind,char* cmd)
{
 printf("\nSETP%s(1,%d,%f,%f,%f)",cmd,ind,x,y,z);
}
/*
*******************************************************************************
* Main procedure                                                              *
*  Arguments:                                                                 *
*    frames: number of objects (frames)                                       *
*        fn: base filename of the objects                                     *
*     timep: time to process                                                  *
*        gp: gravity constant, e.g. -10.0 (negative)                          *
*      etap: friction (Stokes), e.g. -0.00001 (negative)                      *
*      vv0p: speed at t=0 (positive)                                          *
*       fmt: save-format-string, e.g. "TDDD"                                  *
*                                                                             *
*******************************************************************************
*/
/*
*******************************************************************************
* calculates centre of all poarticles                                         *
*******************************************************************************
*/
void getCentre(float ggg)
{
 float xmin,xmax,ymin,ymax,zmin,zmax;
 xmin=pos[0].x;xmax=xmin;ymin=pos[0].y;ymax=ymin;zmin=pos[0].z;zmax=zmin;
 for(i=1;i<oCount;i++) {
  if (xmin>pos[i].x) xmin=pos[i].x;
  if (xmax<pos[i].x) xmax=pos[i].x;
  if (ymin>pos[i].y) ymin=pos[i].y;
  if (ymax<pos[i].y) ymax=pos[i].y;
  if (zmin>pos[i].z) zmin=pos[i].z;
  if (zmax<pos[i].z) zmax=pos[i].z;
 }
 centreX=xmin+(xmax-xmin)/2.0;
 centreY=ymin+(ymax-ymin)/2.0;
 centreZ=zmin+(zmax-zmin)/2.0;
 /* move the object into the positive z-space if necessary */
 if ((zmin<0) && (ggg!=0.0)) {
  xmin=1.1*zmin;
  for(i=0;i<oCount;i++) {
   pos[i].z=pos[i].z-xmin;
  }
  zmin=zmin-xmin;
  zmax=zmax-xmin;
  /* update centreZ */
  centreZ=zmin+(zmax-zmin)/2.0;
  /* change the original object */
  printf("\n;move the object into the positive z-space");
  for(i=0;i<oCount;i++) {
   WritePXXX(pos[i].x,pos[i].y,pos[i].z,i,"POS");
  }
 }
}

void DoIt(int frames,char* fn,float timep,float gp,float etap,float vv0p,char* fmt)
{
 struct Vector *v0;
 float *oldX,*oldY;
 float eta,g,vv0,time,timestep,t,A,Ae,Aet,V0,X0,hh;
 float dx,dy,dz,vb,cx,cy,cz,newX,newY,newZ;
 float e1,e2,e3,e4,e5,e6;
 if (frames<1) return;
 /* allocate some arrays */
 v0=(struct Vector*)calloc(oCount,sizeof(struct Vector));
 if (v0==NULL) return;
 oldX=(float*)calloc(oCount,sizeof(float));
 if (oldX==NULL) {free(v0);return;}
 oldY=(float*)calloc(oCount,sizeof(float));
 if (oldY==NULL) {free(oldX);free(v0);return;}
 time=timep;
 g=gp;
 eta=etap;
 vv0=vv0p;
 if (g>0.0) {printf("\ng-error.\n");g=-10.0;}
 if (vv0<ZERO) {printf("\nvv0-error.\n");vv0=1.0;}
 if (eta>=0.0) {printf("\neta-error.\n");eta=-0.00001;}
 getCentre(g);

 /* compute speed at t=0 (vectors) */
 for(i=0;i<oCount;i++) {
  dx=pos[i].x-centreX;
  dy=pos[i].y-centreY;
  dz=pos[i].z-centreZ;
  vb=dx*dx;
  vb=vb+dy*dy;
  vb=vb+dz*dz;
  vb=sqrt(vb);
  e1=RND2();
  e2=RND2();
  e3=RND2();
  v0[i].x=dx/vb*vv0*e1;
  v0[i].y=dy/vb*vv0*e2;
  v0[i].z=dz/vb*vv0*e3;
 }

 /* generate filename and say "save it" (first and unchanged object) */
 if (!strcmp(fmt,"TDDD")) {printf("\nSavePTDDD(1,");}
 else {
  if (!strcmp(fmt,"VS")) {printf("\nSavePVS(1,");} else {printf("\nSaveP(1,");}
 }
 GenFn(fn,1);
 printf(")");

 /* main loop */
 timestep=time/(frames-1);
 for(i=1;i<frames;i++) {

  /* output frame as comment */
  printf("\n;frame %d",i);

  t=i*timestep; /* actual time */

  for(j=0;j<oCount;j++) {

   /* particle-size (friction) (to calculate this the array PSCL is needed) */
   e1=RND2();
   A=e1*10.0;
   Ae=A*eta;
   Aet=Ae*t;
   hh=exp(Aet);

   /* compute x-coordinate */
   V0=v0[j].x;
   X0=pos[j].x;
   cx=X0*Ae-V0;
   cx=cx+V0*hh;
   cx=cx/Ae;

   /* compute y-coordinate */
   V0=v0[j].y;
   X0=pos[j].y;
   cy=X0*Ae-V0;
   cy=cy+V0*hh;
   cy=cy/Ae;

   /* compute z-coordinate */
   V0=v0[j].z;
   X0=pos[j].z;
   cz=X0*Ae*Ae-g;
   cz=cz-V0*Ae;
   cz=cz+g*hh;
   cz=cz+V0*Ae*hh;  
   cz=cz-g*Aet;
   cz=cz/(Ae*Ae);
   /* motion must stop at z=0 */
   if (cz>0.0) {
    oldX[j]=cx;
    oldY[j]=cy;
   }
   else {
    cz=0.0;
    if (i>1) {
     cx=oldX[j];
     cy=oldY[j];
    }
    else {
     oldX[j]=cx;
     oldY[j]=cy;
    }
   }

   /* the final posion */
   newX=cx;
   newY=cy;
   newZ=cz;

   /* say: "apply them" */
   WritePXXX(newX,newY,newZ,j,"POS");

  }
  
  /* generate filename and say "save it" */
  if (!strcmp(fmt,"TDDD")) {printf("\nSavePTDDD(1,");}
  else {
   if (!strcmp(fmt,"VS")) {printf("\nSavePVS(1,");} else {printf("\nSaveP(1,");}
  }
  GenFn(fn,i+1);
  printf(")");

 }

 /* clean up */
 free(oldY);free(oldX);free(v0);
}

void main()
{
 /* load the ppos-array */
 if (ReadInt("shit:tmp/PEXAMPLE.OCOUNT",&oCount)==0) {
  oCount=oCount/65536;
  pos=(struct Vector*)calloc(oCount,sizeof(struct Vector));
  if (pos==NULL) {
   printf("\nCouldn't allocate PPOS-array.\n");
  }
  else
  { 
   if (ReadReal("shit:tmp/PEXAMPLE.PPOS",&pos[0],oCount)==0) {
    /* tell Dust to load the Source-Particle-Object */
    printf("\nloadp(1,PExample.obj)");
    /* do it !*/
    DoIt(12,"obj",12.8,-10.0,-0.001,52,"TDDD");
    /* important*/
    printf("\n");
    free(pos);
   }
   else
   {
    free(pos);
    printf("\nCouldn't read PPOS-file.\n");
   }
  }
 }
 else
 {
  printf("\nCouldn't read OCOUNT-file.\n");
 }
}