/* BALL.C
   
   Interactive simulation of Free-Throw Space
   
   A project for Math 559, Chaotic Dynamical Systems,
   Professor Stephane Laederich

   For the Silicon Graphics GL package

   Copyright (C) Paul E. Debevec  1992-1998

   debevec@cs.berkeley.edu

   Release 1.0 April 1998

   Compile using the following command:

     cc ball.c -lc_s -lgl_s -lm -o ball

*/


#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include <gl/gl.h>
#include <gl/device.h>

#define BALL_RAD 0.395833333333    /* Ball Radius */
#define RIM_RATIO 0.027777777778   /* rim thickness over rim diameter */
#define XVIEW 1.25
#define YVIEW 0.5
#define ZVIEW 8
#define ZOOMFACTOR 4.0
#define MAXRES 40

/* Define Variable Modeling Parameters */
int courtmod = 1;
int netmod = 1;
int boardmod = 1;
int ballmod = 2;
int hoopmod = 3;
int vecmod = 1;
int postmod = 2;
int viewmod = 2;
int make_movie = 0;

int clongcuts =10;
int hlongcuts =10;
int hlatcuts =10;
int slongcuts =10;
int slatcuts =10;
int numweaves =10;
int numhooks =10;

int stepsperframe = 4;

/* Define pop-up menus */
int mainmenu;
int netmenu;
int viewmenu;
int hoopmenu;
int ballmenu;
int boardmenu;
int vecmenu;
int courtmenu;
int postmenu;
int savemenu;
int stepmenu;

Matrix Identity = {
    1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1
};

/************************/
/* Material Definitions */
/************************/

float vecmat[] = {             /* Vector indicator material */
    AMBIENT, 0.2,0.3,0.6,
    DIFFUSE,0.12,0.269,0.765,
    SPECULAR,0.2,0.2,0.2,
    SHININESS,10,
    LMNULL
};

float ballmat[] = {            /* Basketball material */
    AMBIENT, 0.2,0.1,0.0,
    DIFFUSE,0.72,0.269,0.065,
    SPECULAR,0.1,0.1,0.1,
    SHININESS,10,
    LMNULL
};

float rimmat[] = {             /* Hoop Rim material */
    AMBIENT, 0.2,0.1,0.1,
    DIFFUSE,0.82,0.209,0.065,
    SPECULAR,0.9,0.9,0.9,
    SHININESS,50,
    LMNULL
};

float bfloormat[] = {          /* Court floor material */
    AMBIENT, 0.75,0.4,0.2,
    DIFFUSE,0.2,0.2,0.2,
    SPECULAR,0.1,0.1,0.1,
    SHININESS,10,
    LMNULL
};

float boardmat[] = {           /* Backboard material */
    AMBIENT, 0.2,0.2,0.2,
    DIFFUSE,1,1,1,
    SPECULAR,0.1,0.1,0.1,
    SHININESS,10,
    LMNULL
};

static float lm[] = {         /* Defines the lighting model */
    AMBIENT,0.4,0.4,0.4,
    LOCALVIEWER,1,TWOSIDE,0,LMNULL
};


static float lt0 [] = {         /* defines the lighting source */
    LCOLOR,0.8,0.8,0.8,
    AMBIENT,0.3,0.3,0.3,
    POSITION,0.4,0.5,0.3,0,
    LMNULL
};

static float lt1 [] = {         /* defines the lighting source */
    LCOLOR,0.0,0.0,0.0,
    AMBIENT,0.1,0.1,0.1,
    POSITION,+0.5,1,0.5,0,
    LMNULL
};

/* Store initial positions of the valuators */
int d0init,d1init,d2init,d3init,d4init,d5init,d6init,d7init;

long mapwin,hoopwin;  /* Windows */
long curwin = 0;      /* ID of input focus window */
short val;      /* Used for event input */

/* The colormap and some colors */
static float colormap[256][3];
static float black[3] = {0.0,0.0,0.0};
static float white[3] = {1.0,1.0,1.0};
static float textc[3] = {0.2,0.8,1.0};   /* Used for text */

main()
{
  long xorigin,yorigin,xsize,ysize,xor,yor,ytext;
  float rx,ry,x,y,z;
  double vx,vy,x0,y0;
  short val;      /* Used for event input */

  long p[2];
  int i,j,m,bounces,replot;
  int imax,jmax;
  double zoom, vxcent, vycent, mvx,mvy;
  long key;
  char outstr[80];
  int xmouse,ymouse,oldx,oldy;
  int count,length,dir;
  int numzooms;
  int drawme;

  int pupval;

  void plot(double x,double y,double vx,double vy);
  int dynamic(double x0, double y0, double vx, double vy, int show);
  void cmapinterp(int i1,float r1,float g1,float b1,int i2,float r2,float g2,float b2);
  void initall();
  void menuinit();
  void winsave(char *filename);

  hoopwin = winopen("Free Throw Simulation");
  winposition(1270-800,1270,1000-500,1000);
  
  prefposition(5,260,5,400);
  mapwin = winopen("Free Throw Space");
  qdevice(ESCKEY);
  qdevice(LEFTMOUSE);
  qdevice(RIGHTMOUSE);
  qdevice(MIDDLEMOUSE);
  qdevice(MENUBUTTON);
  RGBmode();
  gconfig();
  clear();

  /* Set up the graphics windows */

  winset(hoopwin);
  getorigin(&xorigin,&yorigin);
  getsize(&xsize,&ysize);
  RGBmode();
  doublebuffer();
  gconfig();
  lsetdepth(getgdesc(GD_ZMIN),getgdesc(GD_ZMAX));
  zbuffer(1);
  mmode(MVIEWING);
  loadmatrix(Identity);
  perspective(450,xsize/(float)ysize,0.1,6000.0);
  lmdef(DEFLIGHT,1,0,lt0);
  lmdef(DEFLIGHT,2,0,lt1);
  lmdef(DEFLMODEL,1,0,lm);
  lmbind(LMODEL,1);
  lmbind(LIGHT0,1);
  lmbind(LIGHT1,2);

  d0init = getvaluator(DIAL0);
  d1init = getvaluator(DIAL1);
  d2init = getvaluator(DIAL2);
  d3init = getvaluator(DIAL3);
  d4init = getvaluator(DIAL4);
  d5init = getvaluator(DIAL5);
  d6init = getvaluator(DIAL6);
  d7init = getvaluator(DIAL7);

  menuinit();
  initall();

  x0=10.0;
  y0=(-3.0);
  vx=(-13);
  vy=15;

  plot(x0,y0,vx,vy);

  winset(mapwin);

  /* Create Color Map */
  cmapinterp(128,1,1,0,255,1,1,0);       /* Megabounce Makes */
  cmapinterp(136,0,0.7,0,154,1,1,1);     /* High-Bounce Makes */
  cmapinterp(129,0,0.15,0,136,0,0.7,0);  /* Low-Bounce Makes */
  
  cmapinterp(0,0,0.3,1,127,0,0.3,1);     /* Megabounce Misses */
  cmapinterp(120,1,0,0,127,0.2,0,0);     /* High-Bounce Misses */
  cmapinterp(107,0,0,1,120,1,0,0);       /* Low-Bounce Misses */

  imax=128;
  jmax=128;
  zoom=0.1;
  vxcent=(-15.0);
  vycent=15.0;
  rectzoom((float)3,(float)3);
  numzooms=0;

  drawme = 1;

  while (1) {

  replot=0;

  getorigin(&xor,&yor);
  count=0;
  length=0;
  dir=0;
  i=imax/2 -1;
  j=jmax/2 -1;
  oldx=oldy=0;
  c3f(white);
  sboxfi(0,256,31,256+31);

  /* Main plotting loop */

  for (m=0;(m<(imax*jmax)) && !replot;m++) {

    vx = vxcent + (2.0*(float)j/jmax -1) / zoom;
    vy = vycent + (2.0*(float)i/imax -1) / zoom;
    bounces=dynamic(x0,y0,vx,vy,0);
    c3f(colormap[bounces + 128]);

    sboxfi(2*j,2*i,2*j+1,2*i+1);

    xmouse=(getvaluator(MOUSEX)-xor)/2*2;
    ymouse=(getvaluator(MOUSEY)-yor)/2*2;

    if (drawme || (curwin==mapwin) && ((oldx != xmouse) || (oldy != ymouse))) {
      /*replot coordinates*/
      drawme = 0;
      rectcopy(xmouse-4,ymouse-4,xmouse+5,ymouse+5,1,257);

      mvx = vxcent + ((float)(xmouse/2) /(jmax/2) -1) / zoom;
      mvy = vycent + ((float)(ymouse/2) /(imax/2) -1) / zoom;

      c3f(black);
      sboxfi(32,256,255,300);

      c3f(textc);
      ytext=395;
      cmov2i(2,ytext-=12);
      sprintf(outstr,"A Chaotic Dynamic Simulation");
      charstr(outstr);
      cmov2i(2,ytext-=12);
      sprintf(outstr,"    of Free-Throw Space");
      charstr(outstr);
      cmov2i(2,ytext-=16);
      sprintf(outstr,"       Paul Debevec");
      charstr(outstr);
      cmov2i(2,ytext-=12);
      sprintf(outstr,"www.cs.berkeley.edu/~debevec");
      charstr(outstr);
      cmov2i(2,ytext-=18);
      sprintf(outstr," Left:Shoot Mid/Right: Zoom");
      charstr(outstr);
      cmov2i(40,ytext=284);
      sprintf(outstr,"Magnification %1.0f%%",zoom*1000);
      charstr(outstr);
      cmov2i(40,ytext-=12);
      sprintf(outstr,"x=%3.15f",mvx);
      charstr(outstr);
      cmov2i(40,ytext-=12);
      sprintf(outstr,"y= %3.15f",mvy);
      charstr(outstr);
      }

    oldx=xmouse;oldy=ymouse;

    if (getbutton(SW31)) {
      winset(hoopwin);
      plot(x0,y0,vx,vy);
      winset(mapwin);
      }

    /* Event Handler */

    if (qtest()) {
    key = qread(&val);

    switch (key) {
    case INPUTCHANGE:
      curwin=val;
      if (curwin==0) { /* Clear out little window */
	c3f(black);
	sboxfi(1,257,29,257+29);
      }
      break;
    case ESCKEY: exit(-1); break;
    case LEFTMOUSE:   /* Recreate the clicked-on shot */
      if (val!=0) break;
      if (curwin==mapwin) {
      winset(hoopwin);
      bounces=dynamic(x0,y0,mvx,mvy,1);
      winset(mapwin);
      if (bounces > 0)  /* Form output string */
        sprintf(outstr,"%d bounce%s, 2 pts.",bounces-1,(bounces==2)?"":"s");
        else if (bounces < 0)
        sprintf(outstr,"%d bounce%s, 0 pts.",-bounces-1,(bounces==(-2))?"":"s");
      c3f(black);
      sboxfi(32,301,255,320);
      cmov2i(40,302);
      c3f(textc);
      charstr(outstr);
      }
      break;
    case MIDDLEMOUSE: case RIGHTMOUSE:  /* Zoom in or out */
      if (curwin==mapwin && val==0) {
        vxcent = vxcent + ((float)(xmouse/2) /(jmax/2) -1) / zoom;
        vycent = vycent + ((float)(ymouse/2) /(imax/2) -1) / zoom;
        zoom*=(key == MIDDLEMOUSE)?ZOOMFACTOR:(double)1.0/ZOOMFACTOR;
        numzooms+=(key == MIDDLEMOUSE)?1:(-1);
        c3f(black);
        clear();
        replot=1;
        oldx=oldy=(-1);  /* force replot upon entry */
      } else {
        if ((key==RIGHTMOUSE) && (val!=0) && (curwin==hoopwin)) {
          pupval = dopup(mainmenu);
          winset(hoopwin);
          plot(x0,y0,vx,vy);
          winset(mapwin);
          }
      }
      break;
    case REDRAW:
      if (val==hoopwin) {
        winset(hoopwin);
        reshapeviewport(); 
        getsize(&xsize,&ysize);
        perspective(450,xsize/(float)ysize,0.1,6000.0);
        plot(x0,y0,vx,vy);
        winset(mapwin);
        }
      else if (val==mapwin) {
        reshapeviewport(); 
        c3f(black);
        clear();
        replot=1;
        }
      break;
    }
    }

    /* Calculate next point in spiral */
    switch(dir) {
      case 0:i++;break;
      case 1:j++;break;
      case 2:i--;break;
      case 3:j--;break;
      }
    count++;
    if (count > length/2) {length++;dir=(dir+1) % 4;count=0;}
    }
  }
}

void winsave(char *filename)
/* Saves the active RGB window to the specified file */
{
  FILE *f;
  unsigned char *buffer;
  unsigned int *lbuffer;
  unsigned int temp;
  long xsize,ysize;
  int numpix,size;
  register i,x,y;

  getsize(&xsize,&ysize);
  numpix=xsize*ysize;
  buffer = malloc(4*numpix);
  lrectread(0,0,xsize-1,ysize-1,(unsigned long *)buffer);

  lbuffer = (unsigned int *)buffer;

  /* Flip upside-down */
  for (y=0;y<ysize/2;y++) {
    for (x=0;x<xsize;x++) {
      temp = lbuffer[y*xsize + x];
      lbuffer[y*xsize + x] = lbuffer[(ysize-y-1)*xsize + x];
      lbuffer[(ysize-y-1)*xsize + x] = temp;
    }
  }
  /* Remap A,B,G,R to R,G,B */
  for (i=0;i<numpix;i++) {
    buffer[3*i] = buffer[4*i+3];
    buffer[3*i+1] = buffer[4*i+2];
    buffer[3*i+2] = buffer[4*i+1];
    }

  f=fopen(filename,"wb");
  fprintf(f,"P6\n# Generated by free-throw space program\n%d %d\n255\n",xsize,ysize);

  size = fwrite(buffer,1,3*xsize*ysize,f);
  printf("%d bytes written to %d x %d image %s.\n",size,xsize,ysize,filename);
  fclose(f);

  free(buffer);
}

void plot(double x,double y,double vx,double vy)
/* Plots the current state of the basketball system */
{
  void drawtorus(double xc, double yc, double zc, double r1);
  void drawboard();
  void drawbfloor();
  void drawnet();
  void drawcyl(double x0, double y0, double z0, double x1, double y1, double z1, double r);
  void drawsphere(double xc, double yc, double zc, double r);
  void drawvec(double x0,double y0,double x1,double y1);

  czclear(0x404040, getgdesc(GD_ZMAX));
  pushmatrix();

  switch (viewmod) {
    case 0:  /* Slam-Cam */
      lookat(1.5,4,0,1.25,0,0,0);
      break;
    case 1:  /* Free-Throw Line */
      lookat(15,-4,0,1.25,0,0,0);
      break;
    case 2:  /* Stands */
      lookat(5,2,18,(x+1.25)/2,y/2,0,0);
      break;
    case 3:  /* Behind Ball */
      lookat(x+2,y,0,x,y,0,0);
      break;
    case 4:  /* Ball-Cam */
      lookat(x,y,0,1.25,0,0,0);
      break;
    case 5:  /* Dial-position view */
      /* Transform coordinates according to dial settings */
      translate((double)(getvaluator(DIAL0)-d0init)/200-XVIEW,
                (double)(getvaluator(DIAL2)-d2init) /200-YVIEW,
                (double)(getvaluator(DIAL4)-d4init) / 200-ZVIEW);
      rotate((double)(getvaluator(DIAL1)-d1init) / 2, 'y');
      rotate((double)(getvaluator(DIAL3)-d3init) / 2, 'z');
      break;
    case 6:  /* Lay-Up */
      lookat(2.63,-2.4,3.1,1,0,0,0);
      break;
    case 7:  /* Demo */
      lookat(6,-3,-4,3,0,0,0);
      break;
    default:
      translate(-XVIEW,-YVIEW,-ZVIEW);
      rotate(0, 'y');
      rotate(0, 'z');
      break;
    }

  lmdef(DEFMATERIAL,1,0,ballmat);
  lmbind(MATERIAL,1);
  if (ballmod != 0) drawsphere(x,y,0,BALL_RAD);

  lmdef(DEFMATERIAL,2,0,rimmat);
  lmbind(MATERIAL,2);
  if (hoopmod !=0) drawtorus(1.25,0,0,0.75);

  lmdef(DEFMATERIAL,3,0,boardmat);
  lmbind(MATERIAL,3);

  if (boardmod != 0) drawboard();
  
  lmbind(MATERIAL,3);
  if (postmod != 0) drawcyl(-4,-10,0,-4,2,0,0.25);

  lmbind(MATERIAL,3);
  if (netmod != 0) drawnet();

  lmdef(DEFMATERIAL,4,0,bfloormat);
  lmbind(MATERIAL,4);
  if (courtmod != 0) drawbfloor();

  lmdef(DEFMATERIAL,5,0,vecmat);
  lmbind(MATERIAL,5);
  if (vecmod !=0) drawvec(x,y,x+vx/10,y+vy/10);

  popmatrix();
  swapbuffers();

}

void cmapinterp(int i1,float r1,float g1,float b1,int i2,float r2,float g2,float b2)
/* User to create the ramped color map */
{
  int i;
  float v;
  float ctemp;

  if (i1==i2) {
    colormap[i1][0]=(r1+r2)/2;
    colormap[i1][2]=(g1+g2)/2;
    colormap[i1][3]=(b1+b2)/2;
    }

  else {

  if (i2<i1) {    /* Swap everything so it isn't backwards */
    i=i1;i1=i2;i2=i;
    ctemp=r1;r1=r2;r2=ctemp;
    ctemp=g1;g1=g2;g2=ctemp;
    ctemp=b1;b1=b2;b2=ctemp;
    }

  for (i=i1;i<=i2;i++) {
    v=(float)(i-i1)/(i2-i1);
    colormap[i][0]=r1 + v*(r2-r1);
    colormap[i][1]=g1 + v*(g2-g1);
    colormap[i][2]=b1 + v*(b2-b1);
    }
  }
}


static float net[MAXRES][MAXRES][3];
static float netn[MAXRES][MAXRES][3];
static float ball[MAXRES][MAXRES+1][3];
static float balln[MAXRES][MAXRES+1][3];
static float hoop[MAXRES][MAXRES+1][3];
static float hoopn[MAXRES][MAXRES+1][3];

/*********************************************************
  Used to specify the graphics primitives via the menu
*********************************************************/

setnet(n)
int n;
{
  void initnet(int model);

  netmod=n;
  initnet(netmod);
  return 8;
}

setball(n)
int n;
{
  void initball(int model);

  ballmod=n;
  initball(ballmod);
  return 8;
}

setboard(n)
int n;
{
  boardmod=n;
  return 8;
}

sethoop(n)
int n;
{
  void inithoop(int model);

  hoopmod=n;
  inithoop(hoopmod);
  return 8;
}

setpost(n)
int n;
{
  postmod=n;
  return 8;
}

setcourt(n)
int n;
{
  courtmod=n;
  return 7;
}

setvec(n)
int n;
{
  vecmod=n;
  return 7;
}

setview(n)
int n;
{
  viewmod=n;
  return 7;
}

setstep(n)
int n;
{
  stepsperframe=n;
  return 7;
}

savescreen(n)
int n;
{
  void winsave(char *filename);

  switch (n) {
  case 0:
    winset(mapwin);
    winsave("map_screen.ppm");
    break;
  case 1:
    winset(hoopwin);
    winsave("hoop_screen.ppm");
    winset(mapwin);
    break;
  case 2:
    make_movie = !make_movie;
    printf("Movie!\n");
    break;
  }
}

void menuinit()
/* Set up the button menus */
{

  ballmenu = defpup("Balls %t %F|None %x0|Coarse %x1|Smooth %x2|Decorated %x3",setball);
  hoopmenu = defpup("Hoops %t %F|None %x0|Simple %x1|Complex %x2|Great %x3",sethoop);
  boardmenu = defpup("Boards %t %F|None %x0|Simple %x1|Complex %x2",setboard);
  netmenu = defpup("Nets %t %F|None %x0|Simple %x1|Complex %x2|Ornate %x3",setnet);
  postmenu = defpup("Posts %t %F|None %x0|Square %x1|Round %x2",setpost);
  stepmenu = defpup("Steps per Frame %t %F|One %x1|Two %x2|Three %x3|Four %x4|Five %x5|Seven %x7|Ten %x10|Fifteen %x15|Twenty-Five %x25",setstep);
  vecmenu = defpup("Vectors %t %F|None %x0|Simple %x1",setvec);
  savemenu = defpup("Save Screen %t %F|Map Screen %x0|Hoop Screen %x1|Make Movie of Next Shot %x2",savescreen);
  courtmenu = defpup("Courts %t %F|None %x0|Simple %x1|Decorated %x2",setcourt);
  viewmenu = defpup("Views %t %F|Slam-Cam %x0|Free-Throw %x1|Layup %x6|Stands %x2|Behind Ball %x3|From Ball %x4|From Dials %x5|Demo %x7",setview);
  
  mainmenu = defpup("Simulation %t|Ball %m|Hoop %m|Backboard %m|Net %m|Post %m|Vector %m|Court %m|View %m|Steps %m|Save %m",
    ballmenu,hoopmenu,boardmenu,netmenu,postmenu,vecmenu,courtmenu,viewmenu,stepmenu,savemenu);
}
/**********************************************************/


void drawcyl(double x0, double y0, double z0, double x1, double y1, double z1, double r)
/* Draws the cylinder graphics primitive */
{
    double phi,dphi = 2 * M_PI/clongcuts;
    float x,y,z;
    float n[3],v[3];
    int i,j;

pushmatrix();
        bgntmesh();
        for (j = 0, phi = 0; j <= clongcuts; j++, phi += dphi) {

	if (j == clongcuts) phi = 0;
	x = cos(phi);
	z = sin(phi);
	n[0] = x;
	n[1] = 0;
	n[2] = z;
	n3f(n);
	v[0] = r*x+x0;
	v[1] = y0;
	v[2] = r*z+z0;
	v3f(v);

	n3f(n);
	v[0] = r*x+x1;
	v[1] = y1;
	v[2] = r*z+z1;
	v3f(v);
    }
    endtmesh();
popmatrix();
}

void drawsphere(double xc, double yc, double zc, double r)
/* Draws the sphere graphics primitive */
{
  int i,j;

  pushmatrix();
  translate(xc,yc,zc);
  scale(r,r,r);
  for (i = 0; i < slatcuts; i++) {
    bgntmesh();
    for (j = 0; j <= slongcuts; j++) {
      n3f(balln[j%slongcuts][i]);
      v3f( ball[j%slongcuts][i]);
      n3f(balln[j%slongcuts][i+1]);
      v3f( ball[j%slongcuts][i+1]);
      }
    endtmesh();
    }
  popmatrix();
}


/* Initialize the coordinate arrays of all the objects */
void initall()
{
  void initball(int model);
  void initnet(int model);
  void inithoop(int model);

  initball(ballmod);
  inithoop(hoopmod);
  initnet(netmod);
}


void initnet(int model)
{
  int hook,weave,i,j;
  double x,y,z,r1,r2;
  double theta,theta2,dtheta,phi,dphi,length;

  numweaves=4;
  numhooks=6;

  if (model >= 2) {
    numweaves=9;
    numhooks=12;
    }

  /* Initialize Net Coordinates */
  dtheta = 2 * M_PI / numhooks;
  for (weave=0;weave < numweaves;weave++) {
    r1=0.75 - 0.3 * sqrt(sqrt(0.5*(double)weave));
    r2=0.012;
    y = (-1.7*(sqrt(weave)) / sqrt(numweaves));
    x = 1.25;
    z = 0;
    for (hook = 0, theta = 0 ; hook < numhooks; hook++, theta += dtheta) {
      theta2=theta + (weave % 2)*0.5*dtheta;
      net[weave][hook][0]=x+r1*cos(theta2);
      net[weave][hook][1]=y;
      net[weave][hook][2]=z+r1*sin(theta2);
      netn[weave][hook][0]=cos(theta2);
      netn[weave][hook][1]=0;
      netn[weave][hook][2]=sin(theta2);
      }
    }
}

void initball(int model)
{
  int i,j;
  double x,y,z,r1,r2;
  double theta,theta2,dtheta,phi,dphi,length;

  slatcuts  = 5;
  slongcuts = 8;

  if (model >= 2) {
    slatcuts=12;
    slongcuts=12;
    }

  /* Initialize Sphere Coordinates */
  dtheta = M_PI/slatcuts;
  dphi = 2 * M_PI/slongcuts;
  for (i = 0, theta = -M_PI/2 ; i <= slatcuts; i++, theta += dtheta) {
    for (j = 0, phi = 0; j < slongcuts; j++, phi += dphi) {
      ball[j][i][0] = balln[j][i][0] = cos(theta) * cos(phi);
      ball[j][i][1] = balln[j][i][1] = sin(theta);
      ball[j][i][2] = balln[j][i][2] = cos(theta) * sin(phi);
      }
    }
}

void inithoop(int model)
{
  int i,j;
  double x,y,z,r1,r2;
  double theta,theta2,dtheta,phi,dphi,length;
  double a,b,c,d,e,f;

  hlatcuts=6;
  hlongcuts=6;

  if (model == 2) {
    hlatcuts=12;
    hlongcuts=12;
    }

  if (model == 3) {
    hlatcuts=36;
    hlongcuts=24;
    }


  /* Initialize Torus Coordinates */
  dtheta = 2 * M_PI/hlatcuts;
  dphi = 2 * M_PI/hlongcuts;
  for (i = 0, theta = 0 ; i < hlatcuts; i++, theta += dtheta) {
    for (j = 0, phi = 0; j <= hlongcuts; j++, phi += dphi) {
	if (j == hlongcuts) phi = 0;
	hoop[j][i][0] = (1 + RIM_RATIO * cos(phi))* cos(theta);
        hoop[j][i][1] = RIM_RATIO * sin(phi);
	hoop[j][i][2] = (1 + RIM_RATIO * cos(phi))* sin(theta);
/*
	hoopn[j][i][0] = cos(theta) * cos(phi);
	hoopn[j][i][1] = sin(phi);
	hoopn[j][i][2] = sin(theta) * cos(phi);
*/
        d = (-sin(theta))*(1+RIM_RATIO*cos(phi));
        e = 0;
        f = cos(theta)*(1+RIM_RATIO*cos(phi));
        a = (-RIM_RATIO)*cos(theta)*sin(phi);
        b = RIM_RATIO*cos(phi);
        c = (-RIM_RATIO)*sin(theta)*sin(phi);
        
        x = b*f-c*e;  y = c*d-a*f;  z = a*e-b*d;

/*        x = RIM_RATIO*cos(phi)*cos(theta)*(1+RIM_RATIO*cos(phi));
        y = (-RIM_RATIO)*cos(theta)*cos(theta)*sin(phi)*(1+RIM_RATIO*cos(phi))
            + RIM_RATIO*sin(theta)*sin(theta)*sin(phi)*(1+RIM_RATIO*cos(phi));
        z = (-RIM_RATIO)*cos(phi)*sin(theta)*(1+RIM_RATIO*cos(phi));
*/
        length = sqrt(x*x + y*y + z*z);
	hoopn[j][i][0] = x / length;
	hoopn[j][i][1] = y / length;
	hoopn[j][i][2] = z / length;
      }
    }
}

void drawnet()
/* Draws the net */
{
  int hook,hook2,weave,strand;
  double radius=0.01;
  double h;
  int numballs;

  void drawsphere(double xc, double yc, double zc, double r);

  for (weave=0;weave < numweaves-1;weave++) {
    for (hook = 0; hook < numhooks; hook++) {
      hook2=(numhooks+hook+((weave%2 != 0)?1:(-1)))%numhooks;
      if (netmod == 3) {
        numballs=(int)(2.0*(net[weave][0][1]-net[weave+1][0][1])/radius);
        for (strand=0;strand<numballs;strand++) {
          h=(double)strand/numballs;
          drawsphere(h*net[weave][hook][0]+
                     (1-h)*net[weave+1][hook][0],
                     h*net[weave][hook][1]+
                     (1-h)*net[weave+1][hook][1],
                     h*net[weave][hook][2]+
                     (1-h)*net[weave+1][hook][2],radius);
          drawsphere(h*net[weave][hook][0]+
                     (1-h)*net[weave+1][hook2][0],
                     h*net[weave][hook][1]+
                     (1-h)*net[weave+1][hook2][1],
                     h*net[weave][hook][2]+
                     (1-h)*net[weave+1][hook2][2],radius);
          
          }
        }
      else {
        drawcyl(net[weave][hook][0],net[weave][hook][1],
                net[weave][hook][2],net[weave+1][hook][0],
                net[weave+1][hook][1],net[weave+1][hook][2],radius);
        drawcyl(net[weave][hook][0],net[weave][hook][1],
                net[weave][hook][2],net[weave+1][hook2][0],
                net[weave+1][hook2][1],net[weave+1][hook2][2],radius);
        }
/*      n3f(netn[weave][hook % numhooks]);
      v3f(net[weave][hook % numhooks]);
      n3f(netn[weave+1][hook % numhooks]);
      v3f(net[weave+1][hook % numhooks]);  */
      if (weave>0) drawsphere(net[weave][hook][0],net[weave][hook][1],net[weave][hook][2],1.75*radius);
      }
    }
}

void drawtorus(double xc, double yc, double zc, double r1)
/* Draws the torus graphics primitive for the hoop */
{
  int i,j;

  pushmatrix();
  translate(xc,yc,zc);
  scale(r1,r1,r1);

  for (i = 0; i < hlatcuts; i++) {
    bgntmesh();
    for (j = 0; j <= hlongcuts; j++) {
      n3f(hoopn[j%hlongcuts][i]);
      v3f(hoop[j%hlongcuts][i]);
      n3f(hoopn[j%hlongcuts][(i+1)%hlatcuts]);
      v3f(hoop[j%hlongcuts][(i+1)%hlatcuts]);
      }
    endtmesh();
    }
  popmatrix();
}

static float vecn[3] = {0.0,0.0,1.0};
void drawvec(double x0,double y0,double x1,double y1)
{
  float vec[3];
  
  bgntmesh();
  n3f(vecn);
  vec[2] = 0;

  vec[0] = x1;
  vec[1] = y1;
  v3f(vec);

  vec[0] = x0-(y1-y0)/12;
  vec[1] = y0+(x1-x0)/12;
  v3f(vec);

  vec[0] = x0+(y1-y0)/12;
  vec[1] = y0-(x1-x0)/12;
  v3f(vec);
  endtmesh();
}

static float bfloorn[3] = {0.0,1.0,0.0};
static float bfloor[4][3] = { {-4.0 , -10.0 , 25.0},
                             {-4.0 , -10.0 , -25.0},
                             {50.0 , -10.0 , -25.0},
                             {50.0 , -10.0 , 25.0}};

static float boardn[3] = {1.0,0.0,0.0};
static float board[4][3] = { {0.0 , -0.5 , -3.0},
                             {0.0 , 4.0 , -3.0},
                             {0.0 , 4.0 ,  3.0},
                             {0.0 , -0.5 ,  3.0}};

static float connectn[3] = {0.0,-1.0,0.0};
static float connect[4][3] = { {0.0 , 0.0 , -0.25},
                               {0.5 , 0.0 , -0.25},
                               {0.5 , 0.0 ,  0.25},
                               {0.0 , 0.0 ,  0.25}};

void drawbfloor()
{

  bgntmesh();
  n3f(bfloorn);
  v3f(bfloor[0]);
  v3f(bfloor[1]);
  v3f(bfloor[3]);
  v3f(bfloor[2]);
  endtmesh();
}

void drawboard()
{

  bgnpolygon();
  /* n3f(boardn); */
  RGBcolor(255,255,255);
  v3f(board[0]);
  v3f(board[1]);
  v3f(board[2]);
  v3f(board[3]);
  endpolygon();
  
  lmdef(DEFMATERIAL,1,0,rimmat);
  lmbind(MATERIAL,1);
  
  bgnpolygon();
  n3f(connectn);
  v3f(connect[0]);
  v3f(connect[1]);
  v3f(connect[2]);
  v3f(connect[3]);
  endpolygon();
  
}

/* Beginning of the dynamic simulation code */

#define a 16.0837   /* g/2 in feet/sec */

#define xloc(t,vx,x0) ((x0)+((t)*(vx)))
#define yloc(t,vy,y0) ((y0)+((t)*(vy))-((a)*(t)*(t)))
#define cdist(x,y,xc,yc,r) (sqrt(((x)-(xc))*((x)-(xc)) + ((y)-(yc))*((y)-(yc))) - (r))

double bdist(double x,double y,double x0, double y0, double x1, double y1)
/* Distance of the ball from a box, ie the backboard */
{
  double maxyet;

  maxyet = x0-x;
  if ((x-x1) > maxyet) maxyet=(x-x1);
  if ((y0-y) > maxyet) maxyet=(y0-y);
  if ((y-y1) > maxyet) maxyet=(y-y1);

  return(maxyet);
}

double objdist(double x, double y, int whichobj)
/* Distance of the ball from a circle */
{
  double dist;

  double bdist(double x,double y,double x0, double y0, double x1, double y1);

  switch (whichobj) {
    case 0: dist = cdist(x,y,2,  0,BALL_RAD); break;
    case 1: dist = cdist(x,y,0.5,0,BALL_RAD); break;
    case 2: dist = cdist(x,y,0  ,4,BALL_RAD); break;
    case 3: dist = bdist(x,y,-BALL_RAD,0,BALL_RAD,4); break;
    }

  return(dist);
}

double Dist(double x, double y, int *closest)
/* Distance of the ball to the hoop & board structure */
{
  double objdist(double x, double y, int whichobj);

  int i;
  double dist,mindist;
  
  mindist = objdist(x,y,0);
  *closest = 0;

  for (i=1;i<4;i++) {
    dist = objdist(x,y,i);
    if (dist < mindist) {
      mindist=dist;
      *closest=i;
      }
    }
  return(mindist);
}

void Bounce(double x,double y,double *vx,double *vy,int closest)
/* Bounces the ball off the board or rim */
{
  double xc,yc,vxt,vyt,r,sinth,costh;

  if (closest == 3) *vx = (- (*vx));
  else {
    switch (closest) {
      case 0: xc=2;yc= 0; break;
      case 1: xc = 0.5;yc=0; break;
      case 2: xc = 0;yc=4; break;
      }

    r = sqrt( (x-xc)*(x-xc) + (y-yc)*(y-yc) );

    sinth = (y-yc) / r;
    costh = (x-xc) / r;

    vxt = (*vx) * costh + (*vy) * sinth;
    vyt = (-(*vx)) * sinth + (*vy) * costh;

    vxt = (-vxt);

    *vx = vxt * costh - vyt * sinth;
    *vy = vxt * sinth + vyt * costh;
    }
}

#define terminated(show,x,y,vx,vy) (((y)<endy) && ((vy)<0))

/* The loop for dynamic simulation

   Returns whether it went in and how many bounces
*/
int dynamic(double x0, double y0, double vx, double vy, int show)
{
  void plot(double x,double y,double vx,double vy);
  void Bounce(double x,double y,double *vx,double *vy,int closest);
  void winsave(char *filename);

  double x,y,t,ta,tb,cvy,begint;
  double oldy;
  double dt = 0.005;
  double firststep,lastplot;
  double endy,dif,olddif,DistA,DistB;
  short val;
  char filename[256];

  int i;
  int iter;
  int closest,dummy;
  int firstone;
  int numbounces;
  int madeit;
  int framenum=0;


  x=x0;
  oldy=y=y0;
  t=0;
  firststep=dt;

  i=0;
  firstone=1;
  numbounces=0;
  madeit=0;
  endy=(show)?(-10+BALL_RAD):(-1);

  /* Advance ball to within hoop range */
  if (!show) {
    begint=(2.6-x0)/vx;
    if (begint>0) {
      while (begint>(t+dt)) t+=dt;
      x=xloc(t,vx,x0);
      y=yloc(t,vy,y0);
      }
    }

  /* Begin main loop */

  while (!terminated(show,x,y,vx,vy-2*a*t)) {

  cvy=vy-2*a*t;

  /***********/
  /* Move ball until intersection detected */
  /***********/
  
  while ((!terminated(show,x,y,vx,cvy)) && (firstone || (Dist(x,y,&closest) > 0))) {
    if (((i++)%stepsperframe == 0) && show) {
      plot(x,y,vx,cvy); /* Display state of system */
      if (make_movie) {
	printf("Making movie frame...\n");
	sprintf(filename,"Movie/movie%04d.ppm",framenum++);
	winsave(filename);
      }
    }
    
    t+=(firstone)?firststep:dt; 
    x=xloc(t,vx,x0);
    y=yloc(t,vy,y0);
    cvy=vy-2*a*t;
    madeit=madeit || ((x>0.5) && (x<2) && (y<0) && (oldy>=0) && (cvy<0));
    firstone=0;
    oldy=y;
    }

  if (!terminated(show,x,y,vx,cvy)) {

  /***********/
  /* Find exact time of intersection with bisection method*/
  /***********/

  ta=t-dt;
  lastplot=tb=t;
  iter=0;

  olddif=0;
  while ((dif=tb-ta) !=olddif) {
    olddif=dif;
    iter++;
    t = 0.5 * (tb+ta);
    if (Dist(xloc(tb,vx,x0),yloc(tb,vy,y0),&dummy) *
        Dist(xloc(t ,vx,x0),yloc(t ,vy,y0),&closest) < 0.0)
         ta = t;
    else tb = t;
    }
  
  x=xloc(t,vx,x0);
  y=yloc(t,vy,y0);
  madeit=madeit || ((x>0) && (x<2) && (y<0) && (oldy>=0) && (vy-2*a*t<0));

  /***********/
  /* Bounce the ball off the object */
  /***********/

  x0=x;
  y0=y;
  vy = vy - 2*a*t;
  firststep = lastplot + dt - t;
  t=0;

  Bounce(x,y,&vx,&vy,closest);
  numbounces++;

  firstone=1;
  }
  }  /* Loop */

  if (madeit) return(numbounces+1);
  else return(-(numbounces+1));
}