/*  This program will display FOUR dimensions on a TWO dimensional screen  */

#include <exec/types.h>
#include <exec/exec.h>
#include <intuition/intuition.h>
#include <intuition/intuitionbase.h>
#include <graphics/gfxbase.h>
#include <graphics/gfx.h>
#include <graphics/text.h>
#include <graphics/regions.h>
#include <graphics/copper.h>
#include <graphics/gels.h>
#include <devices/serial.h>
#include <devices/keymap.h>
#include <stdio.h>
#include <ctype.h>
#include <libraries/dos.h>
#include <math.h>
#ifdef LATTICE		/* Only Lattice currently has limits.h */
#include <limits.h>
#endif
#include <hardware/blit.h>

#define INTUITION_REV 1
#define GRAPHICS_REV  1
#define MILLION       1000000

struct IntuitionBase *IntuitionBase;
struct GfxBase       *GfxBase;
struct Screen        *Screen;
struct Window        *Window;
struct IntuiMessage  *Message;

struct NewScreen NewScreen =
 {
    0,1,640,400,1,1,2,LACE | HIRES, CUSTOMSCREEN,
    NULL,"FOUR DIMENSIONAL POINT MANIPULATION",NULL,NULL,
 };
struct NewWindow NewWindow =
 {
    0,00,639,399,1,1,
    CLOSEWINDOW | MOUSEMOVE | MOUSEBUTTONS,
    WINDOWCLOSE | WINDOWSIZING,
    NULL,NULL,NULL,NULL,NULL,
    0,0,0,0,CUSTOMSCREEN,
 };

/* This is the data all the functions need to see */

#define TOP 100 /* 100 4d points maximum, change at will       */
#define WID 5   /* Width of the 4d matrix x,y,z,t, perspective */
#define DEW 2   /* Width of display data  x,y                  */
#define PLA 6   /* Six Planes in 4 dimensions                  */
#define SIZ 100 /* Size of each line of capture at maximum     */

static double four [ WID ][ WID ];/* This matrix will contain the 4d formula */
static double deta [ TOP ][ WID ];/* This matrix will contain your data      */
          int disp [ TOP ][ WID ];/* Screen coords, note use of WID not DEW  */
          int conn [ 100 ][ 5   ];/* Point connection matrix(NOTE NOT PLANES)*/
static double rads [ PLA ]       ;/* Six rotation planes in four dimensions! */

/*********** RETRIEVE DATA FUNCTION ************/

retrieve (file)
char *file[];
{
     FILE *input;
     char buffer [ SIZ ],rec [ 20 ]; /* 80 chars on one line, 4 sets */
     int a,i,pos;
     pos = FALSE; /* if it stays false program will abort later */
                  /* else it will represent size of data        */

     if ((input = fopen(file,"r"))== NULL) goto HEAVEN;

     while ((fgets(buffer,SIZ,input)))
  {
          a = 0;
          for (i=0;i<WID;i++)
       {
               if (pos > TOP) goto HEAVEN;
               a = strcspn(buffer,",");
               strncpy (rec,buffer,a);
               deta [pos][i] = atof (rec);
               ++a;
               strcpy (buffer,buffer+a);
       }
     ++pos;
  }
HEAVEN:

     fclose (input);
     return (pos-1);
}

/*********** MULTIPLY DATA FUNCTION ************/

multiply (data,P1,P2)
int data;
double P1,P2;
{
    int c,r,i;
    float yx,zx,zy,tx,ty,tz;
    float tweedle;
    float Stz,Sty,Stx,Szy,Szx,Syx,Ctz,Cty,Ctx,Czy,Czx,Cyx;

/* pull desired rotation radians out of array */
    yx = rads [ 0 ];
    zx = rads [ 1 ];
    zy = rads [ 2 ];
    tx = rads [ 3 ];
    ty = rads [ 4 ];
    tz = rads [ 5 ];

/* Pre initialize Sin and Cosine Values, I evaluated and stripped negative Sins */
    Stz=sin(tz);Ctz=cos(tz);
    Sty=sin(ty);Cty=cos(ty);
    Stx=sin(tx);Ctx=cos(tx);

    Szy=sin(zy);Czy=cos(zy);
    Szx=sin(zx);Czx=cos(zx);

    Syx=sin(yx);Cyx=cos(yx);

/*  This is the 4d array, it does all the work
 *  See the program MATFINDER to see how this came into being
 *  At this point we multiply this throught to make a pure numeric -
 *  representation of the 4d matrix for this rotation,
 *  so we can cross multiply it later...
 *  Note again: I stripped all the symbols for negative sin off and -
 *  evaluated those expressions immediately, converting them to positive
 */
    four [0][0] = Ctz*Cty*Ctx;
    four [0][1] = (Stz*Czy-Ctz*Sty*Szy)*Czx-Ctz*Cty*Stx*Szx;
    four [0][2] = 0.0;
    four [0][3] = 0.0;
    four [0][4] = (Stz*Szy+Ctz*Sty*Czy)*Syx+((Stz*Czy-Ctz*Sty*Szy)*Szx
                  +Ctz*Cty*Stx*Czx)*Cyx*P2+(Stz*Szy+Ctz*Sty*Czy)*Cyx
                  -((Stz*Czy-Ctz*Sty*Szy)*Szx+Ctz*Cty*Stx*Czx)*Syx*P1;

    four [1][0] = -Stz*Cty*Ctx;
    four [1][1] = (Ctz*Czy+Stz*Sty*Szy)*Czx+Stz*Cty*Stx*Szx;
    four [1][2] = 0.0;
    four [1][3] = 0.0;
    four [1][4] = (Ctz*Szy-Stz*Sty*Czy)*Syx+((Ctz*Czy+Stz*Sty*Szy)*Szx
                  -Stz*Cty*Stx*Czx)*Cyx*P2+(Ctz*Szy-Stz*Sty*Czy)*Cyx-
                  ((Ctz*Czy+Stz*Sty*Szy)*Szx-Stz*Cty*Stx*Czx)*Syx*P1;

    four [2][0] = -Sty*Ctx;
    four [2][1] = Sty*Stx*Szx-Cty*Szy*Czx;
    four [2][2] = 0.0;
    four [2][3] = 0.0;
    four [2][4] = (Cty*Czy*Syx-(Cty*Szy*Szx+Sty*Stx*Czx)*Cyx)*P2
                  +(Cty*Czy*Cyx+(Cty*Szy*Szx+Sty*Stx*Czx)*Syx)*P1;

    four [3][0] = -Stx;
    four [3][1] = -Ctx*Szx;
    four [3][2] = 0.0;
    four [3][3] = 0.0;
    four [3][4] = Ctx*Czx*Cyx*P2-Ctx*Czx*Syx*P1;

    four [4][0] = 0.0;
    four [4][1] = 0.0;
    four [4][2] = 0.0;
    four [4][3] = 0.0;
    four [4][4] = 1.0;

/* Cross multipy transformation matrix agaist your points */

    for (r=0;r<data;r++)
 {
         for (c=0;c<WID;c++)
      {        
              tweedle = 0.0;
              for (i=0;i<WID;i++)
           {
                   tweedle += (deta [r][i]) * (four [i][c]);
           }
              disp [r][c] = ceil(tweedle);

      }
 }
return (TRUE);
}

main (argc,argv)
short argc;
char *argv[];
{

     double persx,persy;
     short forever = TRUE;
     int x,y,c,pos,top,i,bit,j,zap;
     double pi=3.141592654;

     short INTOPEN = FALSE,
           GFXOPEN = FALSE,
           SCNOPEN = FALSE,
           WINOPEN = FALSE,
           MENOPEN = FALSE;

     if (argc == 1){ printf ("Please supply an argument bitpattern also ex: four 63 (all planes active)\n");
                    goto HELL;}

     bit = atoi(argv[1]); /* get chosen rotations */

     if (!(IntuitionBase = (struct IntuitionBase *)OpenLibrary("intuition.library",INTUITION_REV)))
         goto HELL;INTOPEN = TRUE;

     if (!(GfxBase = (struct GfxBase *)OpenLibrary("graphics.library",GRAPHICS_REV)))
         goto HELL;GFXOPEN = TRUE;

     if (!(Screen = (struct Screen *)OpenScreen(&NewScreen)))
         goto HELL;SCNOPEN = TRUE;

     ShowTitle(Screen,TRUE);
     NewWindow.Screen = Screen;

     if (!(Window = (struct Window *)OpenWindow(&NewWindow)))
         goto HELL;WINOPEN = TRUE;

/* Start of Program */

    if (!(top = retrieve ("fourdata"))) goto HELL;

/* Main Program Loop */

    rads[0]=rads[1]=rads[2]=rads[3]=rads[4]=rads[5]=0.0;

    while (forever)
 {
        if (Message = (struct IntuiMessage *)GetMsg(Window->UserPort))
            if ((Message->Class == CLOSEWINDOW)) goto HELL;
    zap=1;
    for (i=0;i<PLA;i++)
 {
        if ((bit & zap))
      {    rads[i]+=.1;
           if (rads[i]>pi) rads[i]=-pi;
      }
        zap=1;
        for (j=0;j<i;j++) zap = zap*2;
 }

    persx = 1.0;persy = 1.0;

    if (!(multiply (top,persx,persy))) goto HELL;

    SetAPen  (Window->RPort,0);
    RectFill (Window->RPort,0,10,640,400);
    SetAPen  (Window->RPort,3);

    for (pos=0;pos<(top-1);pos++)
  {
        x = disp [ pos ][ 0 ];
        y = disp [ pos ][ 1 ];
        c = disp [ pos ][ 4 ];
        x += 320;y += 200;

        if (pos == 0)
             Move (Window->RPort,x,y);
        Draw (Window->RPort,x,y);
  }
 }
HELL:

      if (MENOPEN)  ClearMenuStrip (Window);
      if (WINOPEN)  CloseWindow    (Window);
      if (SCNOPEN)  CloseScreen    (Screen);
      if (GFXOPEN)  CloseLibrary   (GfxBase);
      if (INTOPEN)  CloseLibrary   (IntuitionBase);
      exit (TRUE);
}