#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define N0 200 /*The default number of bodies in the 0-th group*/
#define N1 100 
#define N2 400 
#define N3 200
#define N4 600

#define X 0
#define Y 1
#define Z 2
#define Vx 3
#define Vy 4
#define Vz 5
#define MASS 6
#define DELTA 360.0
#define ALPHA 0.5
#define BETA 0.5
#define GAMMA 6.67e-11
#define POINTSIZE 1.0e7
#define MASSUNIT 1.0e23
#define XMAX 800.0
#define YMAX 800.0
#define ZMAX 800.0
#define BOXSIZE 120.0
#define EPS 0.001
#define CRITICAL_R 10.0
#define INTERVAL 1
#define MAXGROUPS 9 /*9*/
#define MAXBODIES 600 /*600*/

typedef double Triplet[3];
typedef struct{Triplet pos; Triplet v; double m;} Body;
typedef struct {double x, y, z, vx, vy, vz, m;} Stat;


int [host]M, /*The number of groups*/
    [host]N[MAXGROUPS];/*The numbers of bodies in groups*/
repl dM, dN[MAXGROUPS], di;
double [host]t;
Body (*[host]Galaxy[MAXGROUPS])[MAXBODIES];

nettype GalaxyNet(m, n[m]) {
  coord I=m;
  node {
    I>=0: fast*n[I] scalar;
  };
  /*
  link (J=m){
    J>0: length*(-1) [J]->[0];
    J>0: length*1    [I]->[J];
    };
    */
};

double [host]wtime;
int [host]Input(), [host]DrawGalaxy(void), [host]FirstDrawGalaxy(),
  [host]DrawBye(void);
int [*]MPC_Exit();
double MPC_Wtime();
void Merging(), SingleBody();

void [*]main(int [host]argc, char **[host]argv)
{
  
  if(Input(argv[1])<0)     /*Initializing Galaxy, M and N*/
    {
      ([host]printf)("Error input.\n");
      MPC_Exit(-1);
    }
  
  dM=M;        /*Broadcasting the number of groups*/
  
  dN[]=N[];    /*Broadcasting the numbers of bodies in groups*/
  {
    net GalaxyNet(dM,dN) g;
    int [g]myN, [g]mycoord;
    Body [g]Group[MAXBODIES], [g]OldGroup[MAXBODIES];
    Triplet [g]Centers[MAXGROUPS], **[g]dV;
    double [g]Masses[MAXGROUPS];
    repl [g]group_count, [g]body_count, [g]k;
    void [net GalaxyNet(m, n[m])]MassesIntegrity(double (*)[MAXGROUPS]),
         [net GalaxyNet(m, n[m])]CentersIntegrity(Triplet (*)[MAXGROUPS]);
    int [g]i, [g]j, [g]dim;
    double [g]sigmax, [g]sigmay, [g]sigmaz,
           [g]gmax, [g]gmay, [g]gmaz, [g]F[3], [g]Fabs,
           [g]r, [g]r2, [g]dx, [g]dy, [g]dz;

  
    mycoord = I coordof body_count;
  
    myN = ([g]dN)[mycoord];
    dV=([g]calloc)(myN, [g](sizeof(Triplet *)));
      
    for(body_count=0; body_count<myN; body_count++)
      dV[body_count]=([g]calloc)(myN, [g](sizeof(Triplet)));

    for(group_count=0; group_count<[g]dM; group_count++)  /*Scattering groups*/
        [g:I==group_count]Group[] = (*Galaxy[group_count])[];

    ([g]Merging)(Group, myN);

    wtime=([host]MPC_Wtime)();
    
    if(FirstDrawGalaxy(argc, argv)<0)
      MPC_Exit(-1);
    
    for(body_count=0; body_count<myN; body_count++)
      OldGroup[body_count] = Group[body_count];

    for(body_count=0, Masses[mycoord]=0.0; body_count<myN; body_count++) 
      Masses[mycoord]+=Group[body_count].m;
    
    ([([g]dM,[g]dN)g])MassesIntegrity((void *)Masses);
    
    do {
      Centers[mycoord][]=0.0;
      for(body_count=0; body_count<myN; body_count++) 
        Centers[mycoord][] +=
          (OldGroup[body_count].m/Masses[mycoord])*(OldGroup[body_count].pos)[];
      
      ([([g]dM,[g]dN)g])CentersIntegrity((void *)Centers);
      
      
      for(j=0; j<myN; j++)
        if(OldGroup[j].m>0.0)
          {
            for(k=0, sigmax=0.0, sigmay=0.0, sigmaz=0.0; k<myN; k++)
              if(j<k&&OldGroup[k].m>0.0)
                {
                  Stat [g]jbody;
                  
                  dx=(OldGroup[k].pos)[X]-(OldGroup[j].pos)[X];
                  dy=(OldGroup[k].pos)[Y]-(OldGroup[j].pos)[Y];
                  dz=(OldGroup[k].pos)[Z]-(OldGroup[j].pos)[Z];
                  r2=dx*dx+dy*dy+dz*dz;
                  r=([g]sqrt)(r2);
                  
                  if(r/POINTSIZE<
                     ([g]sqrt)((OldGroup[k].m+OldGroup[j].m)
                               /MASSUNIT)+CRITICAL_R
                       ) {
                    jbody.x=dx;
                    jbody.y=dy;
                    jbody.z=dz;
                    jbody.vx=-(OldGroup[k].v)[X]+(OldGroup[j].v)[X];
                    jbody.vy=-(OldGroup[k].v)[Y]+(OldGroup[j].v)[Y];
                    jbody.vz=-(OldGroup[k].v)[Z]+(OldGroup[j].v)[Z];
                    jbody.m=OldGroup[k].m;
                    ([g]SingleBody)(&jbody);
                    dV[j][k][X]=jbody.vx;
                    dV[j][k][Y]=jbody.vy;
                    dV[j][k][Z]=jbody.vz;
                    sigmax+=dV[j][k][X]/DELTA;
                    sigmay+=dV[j][k][Y]/DELTA;
                    sigmaz+=dV[j][k][Z]/DELTA;
                    }
                    else {
                    Fabs=GAMMA*OldGroup[k].m/r2;
                    sigmax+=Fabs*(dx/r);
                    sigmay+=Fabs*(dy/r);
                    sigmaz+=Fabs*(dz/r);
                    }
              }
              else if(j>k&&OldGroup[k].m>0.0)
                {
                  dx=(OldGroup[k].pos)[X]-(OldGroup[j].pos)[X];
                  dy=(OldGroup[k].pos)[Y]-(OldGroup[j].pos)[Y];
                  dz=(OldGroup[k].pos)[Z]-(OldGroup[j].pos)[Z];
                  r2=dx*dx+dy*dy+dz*dz;
                  r=([g]sqrt)(r2);
                  if(r/POINTSIZE<
                     ([g]sqrt)((OldGroup[k].m+OldGroup[j].m)
                          /MASSUNIT)+CRITICAL_R
                     ) {
                    sigmax-=(OldGroup[k].m/
                             OldGroup[j].m)*dV[k][j][X]/DELTA;
                    sigmay-=(OldGroup[k].m/
                             OldGroup[j].m)*dV[k][j][Y]/DELTA;
                    sigmaz-=(OldGroup[k].m/
                             OldGroup[j].m)*dV[k][j][Z]/DELTA;
                  }
                  else {
                    Fabs=GAMMA*OldGroup[k].m/r2;
                    sigmax+=Fabs*(dx/r);
                    sigmay+=Fabs*(dy/r);
                    sigmaz+=Fabs*(dz/r);
                  }
                }
          for(k=0, gmax=0.0, gmay=0.0, gmaz=0.0; k<[g]dM; k++)
            if(k!=mycoord)
              {
                dx=Centers[k][X]-(OldGroup[j].pos)[X];
                dy=Centers[k][Y]-(OldGroup[j].pos)[Y];
                dz=Centers[k][Z]-(OldGroup[j].pos)[Z];
                r2=dx*dx+dy*dy+dz*dz;
                Fabs=GAMMA*Masses[k]/r2;
                r=([g]sqrt)(r2);
                gmax+=Fabs*(dx/r);
                gmay+=Fabs*(dy/r);
                gmaz+=Fabs*(dz/r);
            
              }
          sigmax=(sigmax+gmax);
          sigmay=(sigmay+gmay);
          sigmaz=(sigmaz+gmaz);
          (Group[j].v)[X]+=sigmax*DELTA;
          (Group[j].v)[Y]+=sigmay*DELTA;
          (Group[j].v)[Z]+=sigmaz*DELTA;
          (Group[j].pos)[X]+=(sigmax*DELTA/2.0+(ALPHA*(Group[j].v)[X]
                             +BETA*(OldGroup[j].v)[X]))*DELTA;
          (Group[j].pos)[Y]+=(sigmay*DELTA/2.0+(ALPHA*(Group[j].v)[Y]
                             +BETA*(OldGroup[j].v)[Y]))*DELTA;
          (Group[j].pos)[Z]+=(sigmaz*DELTA/2.0+(ALPHA*(Group[j].v)[Z]
                             +BETA*((OldGroup[j].v)[Z]-sigmaz*DELTA)))*DELTA;
            }
        ([g]Merging)(Group, myN);
        t+=DELTA;
    
        for(body_count=0; body_count<myN; body_count++)
          OldGroup[body_count] = Group[body_count];
        if(((int)(t/DELTA))%INTERVAL==0)
          {
            for(group_count=0; group_count<[g]dM; group_count++)  /*Gathering groups*/
              (*Galaxy[group_count])[]=[g:I==group_count]Group[];
            if(DrawGalaxy()<0)
              MPC_Exit(-1);
          }
      } while(1);
    DrawBye();
  }
}


void [net GalaxyNet(m, n[m]) g] MassesIntegrity (double (*MassesOfGroups)[MAXGROUPS])
{
  double MassOfMyGroup;
  repl group_count, count;
  MassOfMyGroup=(*MassesOfGroups)[I coordof group_count];
  for(group_count=0; group_count<m; group_count++)
    for(count=0; count<m; count++)
      [g:I==group_count](*MassesOfGroups)[count] = [g:I==count]MassOfMyGroup;
}

void [net GalaxyNet(m, n[m]) g] CentersIntegrity(Triplet (*CentresOfMassOfGroups)[MAXGROUPS])
{
  Triplet CentreOfMassOfMyGroup;
  repl group_count, dim, count;
  CentreOfMassOfMyGroup[]=(*CentresOfMassOfGroups)[I coordof group_count][];
  for(group_count=0; group_count<m; group_count++)
    for(count=0; count<m; count++)
      [g:I==group_count](*CentresOfMassOfGroups)[count][] =
        [g:I==count]CentreOfMassOfMyGroup[];
}

void Merging(Body *rs, int n)
{
  int i, j;
  double dx, dy, dz, r2, r;
  for(i=0; i<n; i++)
    if(rs[i].m>0.0)
      for(j=i+1; j<n; j++)
        if(rs[j].m>0.0)
          {
            dx=(rs[j].pos)[X]-(rs[i].pos)[X];
            dy=(rs[j].pos)[Y]-(rs[i].pos)[Y];
            dz=(rs[j].pos)[Z]-(rs[i].pos)[Z];
            r2=dx*dx+dy*dy+dz*dz;
            r=sqrt(r2);
            if((int)(r/POINTSIZE)<=(int)(sqrt(rs[i].m/MASSUNIT)/2)+
               (int)(sqrt(rs[j].m/MASSUNIT)/2))
              {
                double m;
                m=rs[i].m+rs[j].m;
                (rs[i].pos)[X]=((rs[j].pos)[X]*rs[j].m+(rs[i].pos)[X]*rs[i].m)/m;
                (rs[i].pos)[Y]=((rs[j].pos)[Y]*rs[j].m+(rs[i].pos)[Y]*rs[i].m)/m;
                (rs[i].pos)[Z]=((rs[j].pos)[Z]*rs[j].m+(rs[i].pos)[Z]*rs[i].m)/m;
                (rs[i].v)[X]=((rs[j].v)[X]*rs[j].m+(rs[i].v)[X]*rs[i].m)/m;
                (rs[i].v)[Y]=((rs[j].v)[Y]*rs[j].m+(rs[i].v)[Y]*rs[i].m)/m;
                (rs[i].v)[Z]=((rs[j].v)[Z]*rs[j].m+(rs[i].v)[Z]*rs[i].m)/m;
                rs[i].m=m;
                rs[j].m=0.0;
                (rs[j].pos)[X]=0.0;
                (rs[j].pos)[Y]=0.0;
                (rs[j].pos)[Z]=0.0;
                (rs[j].v)[X]=0.0;
                (rs[j].v)[Y]=0.0;
                (rs[j].v)[Z]=0.0;
              }
          }
}


int [host]Input(char *fname)
{
  FILE *pf;
  int i, j;
  double x0, y0, z0;
  double r, phi, xi;
#define PI2 6.28
  
  if(fname!=NULL)
    {
      pf=fopen(fname, "r");
      if(pf==NULL)
        {
          printf("Can't open file '%s'\n", fname);
          return -1;
        }
      if((i=fscanf(pf, " Number of groups = %d", &M))<1 || M<=0 || M>9)
        {
          fclose(pf);
          if(i<1)
            printf("Problem with reading number of groups from file '%s'\n", fname);
          else
            printf("Number of groups = %d\n", M);
          return -1;
        }
      if((i=fscanf(pf, " Group sizes ="))<0)
        {
          fclose(pf);
          printf("Problem with reading groups' sizes from file '%s'\n", fname);
          return -1;
        }
      for(i=0; i<M; i++)
        if((j=fscanf(pf, " %d", N+i))<1 || N[i]<=0 || N[i]>MAXBODIES)
          {
            fclose(pf);
            if(j<1)
              printf("Problem with reading the size of group #%d from file '%s'\n",
                     i, fname);
            else
              printf("Size of group #%d = %d\n", i, N[i]);
            return -1;
          }
      fclose(pf);
      for(i=0; i<M; i++) {
        Galaxy[i]=calloc(MAXBODIES, sizeof(Body));
        if(Galaxy[i]==NULL)
          {
            printf("Can't allocate storage for Galaxy[%d]\n", i);
            return -1;
          }
      }
      
      for(i=0; i<M; i++)
        {
          switch(i) {
          case 0: x0=y0=BOXSIZE/2.0; break;
          case 1: x0=BOXSIZE*1.5; y0=YMAX-BOXSIZE/2.0; break;
          case 2: x0=XMAX/2.0; y0=YMAX/2.0; break;
          case 3: x0=XMAX-BOXSIZE/2.0; y0=BOXSIZE/2.0; break;
          case 4: x0=XMAX-BOXSIZE/2.0; y0=YMAX-BOXSIZE/2.0;
                  break;
          case 5: x0=BOXSIZE*1.5; y0=YMAX/2.0; break;
          case 6: x0=XMAX-BOXSIZE/2.0; y0=YMAX/2.0; break;
          case 7: x0=XMAX/2.0; y0=YMAX-BOXSIZE/2.0; break;
          case 8: x0=XMAX/2.0; y0=BOXSIZE/2.0; break;
          }
          for(j=0; j<N[i]; j++)
            {
              r=(BOXSIZE/2.0/RAND_MAX)*rand();
              phi=(PI2/RAND_MAX)*rand();
              ((*Galaxy[i])[j].pos)[X]=(x0+r*cos(phi))*POINTSIZE;
              ((*Galaxy[i])[j].pos)[Y]=(y0+r*sin(phi))*POINTSIZE;
              ((*Galaxy[i])[j].pos)[Z]=0.0;
              ((*Galaxy[i])[j].v)[X]=0.0;
              ((*Galaxy[i])[j].v)[Y]=0.0;
              ((*Galaxy[i])[j].v)[Z]=0.0;
              (*Galaxy[i])[j].m=(j==N[i]/2)?50.0*MASSUNIT:0.5*MASSUNIT;
            }
        }
    }
  else
    {
      M=5;
      N[0]=N0;
      N[1]=N1;
      N[2]=N2;
      N[3]=N3;
      N[4]=N4;
      for(i=0; i<M; i++)
        Galaxy[i]=calloc(MAXBODIES, sizeof(Body));
      
      for(i=0; i<M; i++)
        {
          switch(i) {
          case 0: x0=y0=BOXSIZE/2.0; break;
          case 1: x0=BOXSIZE*1.5; y0=YMAX-BOXSIZE/2.0; break;
          case 2: x0=XMAX/2.0; y0=YMAX/2.0; break;
          case 3: x0=XMAX-BOXSIZE/2.0; y0=BOXSIZE/2.0; break;
          case 4: x0=XMAX-BOXSIZE/2.0; y0=YMAX-BOXSIZE/2.0;
                  break;
          }
          for(j=0; j<N[i]; j++)
            {
              r=(BOXSIZE/2.0/RAND_MAX)*rand();
              phi=(PI2/RAND_MAX)*rand();
              ((*Galaxy[i])[j].pos)[X]=(x0+r*cos(phi))*POINTSIZE;
              ((*Galaxy[i])[j].pos)[Y]=(y0+r*sin(phi))*POINTSIZE;
              ((*Galaxy[i])[j].pos)[Z]=0.0;
              ((*Galaxy[i])[j].v)[X]=0.0;
              ((*Galaxy[i])[j].v)[Y]=0.0;
              ((*Galaxy[i])[j].v)[Z]=0.0;
              (*Galaxy[i])[j].m=(j==N[i]/2)?50.0*MASSUNIT:0.5*MASSUNIT;
            }
        }
    }
    return 1;
}

void SingleBody(Stat *body)
{
  double r2, r, a, ax, ay, az, dt=DELTA,
    t, xk1, yk1, zk1, xk2, yk2, zk2,
    vxk1, vyk1, vzk1, vxk2, vyk2, vzk2,
    axk, ayk, azk;
  int i=1;
  
  r2=body->x*body->x+body->y*body->y+body->z*body->z;
  a=GAMMA*body->m/r2;
  r=sqrt(r2);
  ax=a*(body->x/r);
  ay=a*(body->y/r);
  az=a*(body->z/r);
  vxk1=body->vx+ax*dt;
  vyk1=body->vy+ay*dt;
  vzk1=body->vz+az*dt;
  xk1=body->x+vxk1*dt/2.0+body->vx*dt;
  yk1=body->y+vyk1*dt/2.0+body->vy*dt;
  zk1=body->z+vzk1*dt/2.0+body->vz*dt;
  
l:
  for(dt/=2.0, t=0.0, vxk2=body->vx, vyk2=body->vy, vzk2=body->vz,
        xk2=body->x, yk2=body->y, zk2=body->z, axk=ax, ayk=ay, azk=az;
      t<DELTA;
      t+=dt
      ) {
    i++;
    vxk2+=axk*dt;
    xk2+=vxk2*dt/2.0+(ALPHA*vxk2+BETA*(vxk2-axk*dt))*dt;
    vyk2+=ayk*dt;
    yk2+=vyk2*dt/2.0+(ALPHA*vyk2+BETA*(vyk2-ayk*dt))*dt;
    vzk2+=azk*dt;
    zk2+=vzk2*dt/2.0+(ALPHA*vzk2+BETA*(vzk2-azk*dt))*dt;
    r2=xk2*xk2+yk2*yk2+zk2*zk2;
    r=sqrt(r2);
    a=GAMMA*body->m/r2;
    axk=a*(xk2/r);
    ayk=a*(yk2/r);
    azk=a*(zk2/r);
  }
  if(fabs(vxk2-vxk1)<=EPS*fabs(vxk2) &&
     fabs(vyk2-vyk1)<=EPS*fabs(vyk2) &&
     fabs(vzk2-vzk1)<=EPS*fabs(vzk2) &&
     fabs(xk2-xk1)<=EPS*fabs(xk2) &&
     fabs(yk2-yk1)<=EPS*fabs(yk2) &&
     fabs(zk2-zk1)<=EPS*fabs(zk2)
     )
        {
          body->x=xk2;
          body->y=yk2;
          body->z=zk2;
          body->vx=vxk2-body->vx;
          body->vy=vyk2-body->vy;
          body->vz=vzk2-body->vz;
          return;
        }
  else {
    vxk1=vxk2;
    vyk1=vyk2;
    vzk1=vzk2;
    xk1=xk2;
    yk1=yk2;
    zk1=zk2;
    goto l;
  }
}

#pragma keywords ANSI
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>


#include <X11/bitmaps/icon>
#pragma keywords SHORT

#define WIDTH XMAX
#define HEIGHT YMAX

Display *[host]display;
int [host]screen;
Window [host]win;
GC [host]gc;
XEvent [host]report;
XFontStruct *[host]font_info;
char *[host]fontname="9x15";
void [host]DrawBody();
void [host]DrawText();
void [host]SayGoodBye();


int [host]FirstDrawGalaxy(int argc, char **argv)
{
  unsigned width=WIDTH, height=HEIGHT;
  int x=0, y=0;
  unsigned border_width=4;
  unsigned display_width, display_height;
  char *window_name="Galaxy Demo";
  char *icon_name="galaxy";
  Pixmap icon_pixmap;
  XSizeHints size_hints;
  unsigned long valuemask=0;
  XGCValues values;
  char *display_name=NULL;

  int i, j, k;

  if((display=XOpenDisplay(display_name))==NULL)
    {
      printf("FirstDrawGalaxy: can't connect to Xserver %s\n",
                     XDisplayName(display_name));
      return (-1);
    }
  screen=XDefaultScreen(display);

  win=XCreateSimpleWindow(display, XRootWindow(display,screen),
                          x, y, width, height, border_width,
                          XBlackPixel(display, screen),
                          XBlackPixel(display, screen)
                          );
  icon_pixmap=XCreateBitmapFromData(display, win, icon_bits,
                                    icon_width, icon_height);
  size_hints.flags=PPosition|PSize|PMinSize;
  size_hints.x=x;
  size_hints.y=y;
  size_hints.width=width;
  size_hints.height=height;
  size_hints.min_width=width;
  size_hints.min_height=height;

  XSetStandardProperties(display, win, window_name, icon_name,
                         icon_pixmap, argv, argc, &size_hints);
  XSelectInput(display, win,
               ExposureMask|KeyPressMask|ButtonPressMask|StructureNotifyMask);

  if((font_info=XLoadQueryFont(display, fontname))==NULL)
    {
      printf("FirstDrawGalaxy: can't open %s font\n", fontname);
      return -1;
    }
  
  gc=XCreateGC(display, win, valuemask, &values);
  XSetForeground(display, gc, XWhitePixel(display, screen));
  XSetLineAttributes(display, gc, 1, LineSolid, CapNotLast, JoinMiter);
  XSetFillStyle(display, gc, FillSolid);
  XSetFont(display, gc, font_info->fid);
  
  XMapWindow(display, win);
  while(1) {
    XNextEvent(display, &report);
    switch(report.type) {
    case Expose:
      while (XCheckTypedEvent(display, Expose, &report));
      for(i=0; i<M; i++) {
        for(j=0; j<N[i]; j++) {
          if((*Galaxy[i])[j].m>0.0)
            DrawBody(&gc,
                     (int)((*Galaxy[i])[j].m/MASSUNIT),
                     (int)(((*Galaxy[i])[j].pos)[X]/POINTSIZE),
                     (int)(((*Galaxy[i])[j].pos)[Y]/POINTSIZE),
                     (int)(((*Galaxy[i])[j].pos)[Z]/POINTSIZE));
        }
      }
      DrawText([host]t);
      XFlush(display);
      return 1;
    case ConfigureNotify: break;
    case ButtonPress:
    case KeyPress:
      return 1;
    default:
      break;
    }
  }
}

int [host]DrawGalaxy(void)
{
  int i, j, k, size;
  XEvent event;

  event.type=Expose;
  XSendEvent(display, win, False, 0L, &event);
  while(1) {
    XNextEvent(display, &report);
    switch(report.type) {
    case Expose:
      while (XCheckTypedEvent(display, Expose, &report));
      XClearWindow(display, win);
      for(i=0; i<M; i++) {
        for(j=0; j<N[i]; j++) {
          if((*Galaxy[i])[j].m>0.0)
            DrawBody(&gc,
                   (int)((*Galaxy[i])[j].m/MASSUNIT),
                   (int)(((*Galaxy[i])[j].pos)[X]/POINTSIZE),
                   (int)(((*Galaxy[i])[j].pos)[Y]/POINTSIZE),
                   (int)(((*Galaxy[i])[j].pos)[Z]/POINTSIZE));
        }
      }
      DrawText([host]t);
      XFlush(display);
      XSendEvent(display, win, False, 0L, &event);
      return 1;
    case ConfigureNotify: break;
    case ButtonPress:
    case KeyPress:
      XSendEvent(display, win, False, 0L, &event);
      DrawBye();
      return -1;
    default:
      break;
    }
  }
}

int [host]DrawBye(void)
{
  int i, j;
  while(1) {
    XNextEvent(display, &report);
    switch(report.type) {
    case Expose:
      while (XCheckTypedEvent(display, Expose, &report));
      XClearWindow(display, win);
      for(i=0; i<M; i++) {
        for(j=0; j<N[i]; j++)
          if((*Galaxy[i])[j].m>0.0)
            DrawBody(&gc,
                   (int)((*Galaxy[i])[j].m/MASSUNIT),
                   (int)(((*Galaxy[i])[j].pos)[X]/POINTSIZE),
                   (int)(((*Galaxy[i])[j].pos)[Y]/POINTSIZE),
                   (int)(((*Galaxy[i])[j].pos)[Z]/POINTSIZE));
        DrawText([host]t);
        SayGoodBye();
        XFlush(display);
      }
      break;
    case ConfigureNotify: break;
    case ButtonPress:
    case KeyPress:
      XFreeGC(display, gc);
      XCloseDisplay(display);
      return 1;
    default:
      break;
    }
  }
}

void [host]DrawBody(GC *gc, int m, int x, int y, int z)
{
  int xx, yy, r;
  
  r=(int)(sqrt((double)m)/2);
  for(xx=x-r; xx<=x+r; xx++)
    for(yy=y-r; yy<=y+r; yy++)
      if((xx-x)*(xx-x)+(yy-y)*(yy-y)<=r*r)
        XDrawPoint(display, win, *gc, xx, yy);
}

void [host]DrawText(double t)
{
  char text[25];
  int len, wid;
  double gt, wt;
  
  sprintf(text, "Galaxy time = %.1f hours", (gt=t/3600.0));
  len=strlen(text);
  wid=XTextWidth(font_info, text, len);
  XDrawString(display, win, gc, (int)((WIDTH-wid)/2), 20, text, len);
  sprintf(text, "Wall time = %.1f seconds", (wt=MPC_Wtime()-wtime));
  len=strlen(text);
  wid=XTextWidth(font_info, text, len);
  XDrawString(display, win, gc, (int)((WIDTH-wid)/2), 40, text, len);
  sprintf(text, "Rate = %.1f sec/h", gt?wt/gt:0.0);
  len=strlen(text);
  wid=XTextWidth(font_info, text, len);
  XDrawString(display, win, gc, (int)((WIDTH-wid)/2), 60, text, len);
}

void [host]SayGoodBye(void)
{
  char text[50];
  int len, wid;
  
  sprintf(text, "Good bye! Hope to see you soon.");
  len=strlen(text);
  wid=XTextWidth(font_info, text, len);
  XDrawString(display, win, gc, (int)((WIDTH-wid)/2), (int)(HEIGHT/3),
                      text, len);
  XFlush(display);
}