/* Galaxy colisions in 3D © 1989 AMIGA Plus
   Monday 19-Feb-90 09:50:27
   This program is written in C based on the BASIC
   program listed in the December, 1988 issue of "Astronomy"
   Compiled with Lattice 5.02
   Linked with a modified _main.o which opens a large CON: window
   if GC3D is run from the Workbench.
   This file (_main.o) is on this AMIGA Plus disk.  The original source
   is in the file umain.c in the lc/source directory of Lattice C
Compile: (Lc gc3d)
Link: (Blink lib:c.o,_main.o,gc3d.o to gc3d lib lib:lcm.lib,lc.lib,amiga.lib)
*/

#include <stdio.h>
#include <math.h>
#include <exec/types.h>
#include <exec/memory.h>
#include <intuition/intuition.h>

#define NUMSTARS 1000           /* Max number of stars */
#define SF2 2.0                 /* Softening factor */
#define ASPECT 0.44             /* Aspect ratio of a high res screen */
#define CXY_X 160               /* Center of X-Y plane, X */
#define CXY_Y 90                /* Center of X-Y plane, Y */
#define CXZ_X 480               /* Center of X-Z plane, X */
#define CXZ_Y 90                /* Center of X-Z plane, Y */
#define BLUP 0                  /* Color definitions */
#define WHTP 1
#define BLKP 2
#define ORGP 3

struct GfxBase *GfxBase;
struct IntuitionBase *IntuitionBase;
struct NewWindow nwindow = {
        0, 0,                   /* Left/Top Edges */
        640, 200,               /* Width/Height */
        ORGP, BLKP,             /* Detail Pen/Block Pen */
        MOUSEBUTTONS | RAWKEY | CLOSEWINDOW,
                                /* IDCMP Flags */
        SMART_REFRESH | ACTIVATE | BORDERLESS | WINDOWDEPTH | GIMMEZEROZERO |
        WINDOWCLOSE,            /* Flags */
        NULL,                   /* No Gadgets */
        NULL,                   /* Checkmark */
        " Galaxy Collision",    /* Title */
        NULL,                   /* Screen (dummy) */
        NULL,                   /* BitMap */
        100, 50, 640, 200,      /* MinWidth, MinHeight, MaxWidth, MaxHeight */
        WBENCHSCREEN            /* Screen Type */
};
struct Window *w;               /* Our Window pointer */
float x[NUMSTARS], y[NUMSTARS], z[NUMSTARS],
        vx[NUMSTARS], vy[NUMSTARS], vz[NUMSTARS];
int nrr;                        /* Number of rings of stars */
int nrs;                        /* Number of stars per ring */
int ns;                         /* Total number of stars */
float dr;                       /* Ring spacing */
float x1,y1,z1;                 /* Position of Galaxy #1 */
float x2,y2,z2;                 /* Position of Galaxy #2 */
float vx1,vy1,vz1;              /* Velocity of Galaxy #1 */
float vx2,vy2,vz2;              /* Velocity of Galaxy #2 */
float m1,m2;                    /* Masses of the Galaxies */
int time;                       /* Time */
short quit;                     /* Quit flag */
short pause_key = 0;            /* Pause flag */

void main(), drawstars(), done(), handle_messages();

void main(ac,av)

int ac;
char *av[];

{
        char  str[200];
        short ir, it, i, j, k;
        char ans[100];
        float r, v, th, t, t1, r1, r2, ax, ay, az, rr;;

  if ((IntuitionBase=(struct IntuitionBase *)OpenLibrary("intuition.library",1))==NULL)
    done("Could not open intuition.library");
  if ((GfxBase=(struct GfxBase *)OpenLibrary("graphics.library",0))==NULL)
    done("Could not open graphics.library");

tryagain:

  printf("\nWelcome to GC3D, the 3D Galaxy Collision Simulator\n\n");
  printf("This program will ask you the number of star rings and the number of\n");
  printf("stars per ring. These two figures multiplied must be less than");
  printf(" %d.\n",NUMSTARS);
  printf("Pressing RETURN for each question selects the default values.\n");

  printf("\nInput the number of rings of stars in the TARGET galaxy (default 10) ? ");
        gets(str);
        if (sscanf(str," %d\n",&nrr) != 1) nrr=10;
  printf("Input the number of stars per ring in the TARGET galaxy (default 10) ? ");
        gets(str);
        if (sscanf(str," %d\n",&nrs) != 1) nrs=10;

        ns = nrr*nrs;
        if (ns > NUMSTARS) {
                printf("Too many stars specified, try again.");
                goto tryagain;
        }
        if (nrr > 1)
                dr = 20.0/(nrr-1);
        else
                dr = 20;
  printf("Input mass fraction of INTRUDER galaxy in units of the target\n");
  printf("galaxy mass (default 1) ? ");
        gets(str);
        if (sscanf(str," %f\n",&m2) != 1) m2=1;

  printf("\nThe TARGET galaxy is located at 0,0,0.\n");
  printf("Input the initial X,Y,Z coordinates of the INTRUDER galaxy.\n");
  printf("(e.g. 15 15 0) (Default 7.5 0 35) ? ");
        gets(str);
        if (sscanf(str," %f %f %f\n",&x2,&y2,&z2) != 3) {
          x2=7.5;
          y2=0;
          z2=35;
        }
  printf("\nThe TARGET galaxy is initially at rest:\n");
  printf("Input the initial X,Y,Z velocities of the INTRUDER galaxy.\n");
  printf("(e.g. -1 -1 0) (Default 0 0 -1) ? ");
        gets(str);
        if (sscanf(str," %f %f %f\n",&vx2,&vy2,&vz2) != 3) {
          vx2=0;
          vy2=0;
          vz2=-1;
        }

        printf("\nNumber of Rings of stars: %d\n",nrr);
        printf("Number of Stars per ring: %d\n",nrs);
        printf("Mass fraction of INTRUDER galaxy: %f\n",m2);
        printf("Initial location of INTRUDER galaxy: %f %f %f\n",x2,y2,z2);
        printf("Initial velocity of INTRUDER galaxy: %f %f %f\n",vx2,vy2,vz2);
        printf("\nAre these inputs correct (default Y) ? ");
        gets(str);
        if (sscanf(str," %s",ans) != 1) ans[0]='Y';
        if (toupper(ans[0]) != 'Y')
                goto tryagain;
        printf("Thinking. . .");

/* Now initialize the target galaxy's mass, position, and velocity */

        m1 = 5.0;
        x1 = 160.0;
        y1 = 100.0;
        vx1 = vy1 = vz1 = z1 = 0.0;

/* Scale the intruder galaxy's mass , postion, and velocities */

        time = 0;
        m2 *= m1;
        x2 += x1;
        y2 += y1;
        z2 += z1;

/* Place stars in concentric rings around TARGET galaxy */

        for (ir=i=0; ir < nrr; ir++) {
                r =10.0+ir*dr;
                v = sqrt(m1/r);
                th = (0.5*v/r)*(180.0/PI);
                if (r==10.0)
                        v *= 0.9;
                for (it=0; it <nrs; it++,i++) {
                        t = it*360.0/nrs;
                        t1 = PI*(t-th)/180.0;

/* Initialize the star positions and velocities */

                        x[i] = r*cos(t/57.2958)+x1;
                        y[i] = r*sin(t/57.2958)+y1;
                        z[i] = vz[i] = 0.0;
                        vx[i] = -v*sin(t1);
                        vy[i] = v*cos(t1);
                }
        }

/* Now open the window and print the initial positions out */

        if ((w=(struct Window *)OpenWindow(&nwindow))==NULL)
                done("Could not open window.");
        drawstars();

/* Now loop */

        printf("\n");
        for (k=0,quit=0; !quit; k++) {  /* Loop Until closed */
                for (j=0; j < 1; j++) {
                        for (i=0; i < ns; i++) {

/* Determine the force on a star due to the galactic centers.  SF2 below
   called the softening factor, is used to prevent overflows during force
   calculation.  SF2 is assigned a value above.
*/

                                r1 = m1/pow(pow(x[i]-x1,2.0)+pow(y[i]-y1,2.0)+pow(z[i]-z1,2.0)+SF2,1.5);
                                r2 = m2/pow(pow(x[i]-x2,2.0)+pow(y[i]-y2,2.0)+pow(z[i]-z2,2.0)+SF2,1.5);
                                ax = r1*(x1-x[i])+r2*(x2-x[i]);
                                ay = r1*(y1-y[i])+r2*(y2-y[i]);
                                az = r1*(z1-z[i])+r2*(z2-z[i]);
                                vx[i] += ax;
                                vy[i] += ay;
                                vz[i] += az;
                                x[i] += vx[i];
                                y[i] += vy[i];
                                z[i] += vz[i];
                        }

/* Update positions and velocities of galactic centers */

                        rr = pow(pow(x1-x2,2.0)+pow(y1-y2,2.0)+pow(z1-z2,2.0)+SF2,1.5);
                        ax = (x2-x1)/rr;
                        ay = (y2-y1)/rr;
                        az = (z2-z1)/rr;
                        vx1 += m2*ax;
                        vy1 += m2*ay;
                        vz1 += m2*az;
                        vx2 -= m1*ax;
                        vy2 -= m1*ay;
                        vz2 -= m1*az;
                        x1 += vx1;
                        y1 += vy1;
                        z1 += vz1;
                        x2 += vx2;
                        y2 += vy2;
                        z2 += vz2;
                        time++;
                }
                drawstars();
                handle_messages();
                while (!quit && pause_key) {
                        Wait(1<<w->UserPort->mp_SigBit);
                        handle_messages();
                }
        }
        CloseWindow(w);
        done(NULL);
}

/* Handle messages from user port */

void handle_messages()

{
        struct IntuiMessage *msg;
        short cls, cde, msey;

        while (msg=(struct IntuiMessage *)GetMsg(w->UserPort)) {
                cls = msg->Class;
                cde = msg->Code;
                msey = msg->MouseY;
                switch (cls) {
                        case CLOSEWINDOW:
                                quit++;
                        case MOUSEBUTTONS:
                                if ((cde == SELECTDOWN) && ((msg->MouseY+
                                        msg->MouseX) < 10))
                                        quit++;
                        case RAWKEY:
                                if (msg->Code == 0x19)  /* 'P'ause key */
                                        pause_key ^= 1;
                }
                ReplyMsg(msg);
        }
}

/* Come here to draw the stars on the screen */

void drawstars()

{
        float xc, yc, zc, xx1, yy1, zz1, xx2, yy2, zz2;
        short i, xi, yi, zi;
        char tbuf[100];

/* Clear the screen and set up the lines */

        SetRast(w->RPort,BLKP);
        SetDrMd(w->RPort,JAM1);
        SetAPen(w->RPort,BLUP);
        Move(w->RPort,0,180);
        Draw(w->RPort,639,180);
        Move(w->RPort,320,0);
        Draw(w->RPort,320,180);
        Move(w->RPort,150,189);
        sprintf(tbuf,"X - Y           Time = %4d             X - Z",time);
        Text(w->RPort,tbuf,strlen(tbuf));

/* Calculate the center of mass of the system */

        xc = (m1*x1+m2*x2)/(m1+m2);
        yc = (m1*y1+m2*y2)/(m1+m2);
        zc = (m1*z1+m2*z2)/(m1+m2);

/* Calculate position of galactic centers and display on screen */

        xx1 = x1-xc;
        yy1 = y1-yc;
        zz1 = z1-zc;
        xx2 = x2-xc;
        yy2 = y2-yc;
        zz2 = z2-zc;
        xi = (short)(xx2/ASPECT);
        yi = (short)yy2;
        zi = (short)zz2;
        SetAPen(w->RPort,ORGP);
        if (xi>-160 && xi<160 && yi>-90 && yi<90)
                DrawEllipse(w->RPort,CXY_X+xi,CXY_Y-yi,5,2);
        if (xi>-160 && xi<160 && yi>-90 && yi<90)
                DrawEllipse(w->RPort,CXZ_X+xi,CXZ_Y-zi,5,2);
        xi = (short)(xx1/ASPECT);
        yi = (short)yy1;
        zi = (short)zz1;
        SetAPen(w->RPort,WHTP);
        if (xi>-160 && xi<160 && zi>-90 && zi<90)
                DrawEllipse(w->RPort,CXY_X+xi,CXY_Y-yi,9,4);
        if (xi>-160 && xi<160 && zi>-90 && zi<90)
                DrawEllipse(w->RPort,CXZ_X+xi,CXZ_Y-zi,9,4);

/* Now set all the stars in onto the display*/

        for (i=0; i<ns; i++) {
                xi = (short)((x[i]-xc)/ASPECT);
                yi = (short)(y[i]-yc);
                zi = (short)(z[i]-zc);
                if (xi>-160 && xi<160 && yi>-90 && yi<90)
                        WritePixel(w->RPort,CXY_X+xi,CXY_Y-yi);
                if (xi>-160 && xi<160 && zi>-90 && zi<90)
                        WritePixel(w->RPort,CXZ_X+xi,CXZ_Y-zi);
        }
}

/* Come here to exit */

void done(s)

char *s;

{
        if (s)
                printf("%s\n");
        exit(1);
}