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

void main(void);
int inlines(double x,double y, double z, double w);

/* These are constants used to jimmy the 
   basketball parameters. */

#define BBRADIUS     123.3

#define DOTRADIUS    (2.55/BBRADIUS)
#define DOTSPACING   (2.9/BBRADIUS)
#define NOSPOTWIDTH  (7.566666/BBRADIUS)
#define COLORWIDTH   (4.775/BBRADIUS)
#define CIRCPOS       66.5
#define CIRCRAD       47.0
#define PIPCLOSE      .99994652

#define OUTWIDTH 1024
#define OUTHEIGHT 512

void main(void)
{
 unsigned char g[OUTWIDTH*OUTHEIGHT];
 int t,i,p;
 double rad,x,y,z,phi,theta,r;
 double dotx[19000],doty[19000],dotz[19000]; 
     /* Ugly hardcoded arrays */
 int numspots;
 FILE *infile;
 
 infile=fopen("spherepoints","r"); 
        /* file of spherically distributed points */
 p=0;
 while (fscanf(infile,"%lg%lg%lg",&dotz[p],
               &doty[p],&dotx[p])!=EOF)
   {
     /* Normalize points to be at surface of the sphere */
     rad=1.0/sqrt(dotx[p]*dotx[p]+doty[p]*doty[p]+
                  dotz[p]*dotz[p]);
     dotx[p]*=rad;
     doty[p]*=rad;
     dotz[p]*=rad;
     /* If the points are not too close to one of the lines
        of the ball, add them to the list of pips. */
     if (!inlines(dotx[p],doty[p],dotz[p],
         NOSPOTWIDTH/2.0+DOTRADIUS)) p++; 
   }

 fclose(infile);
 numspots=p;

 for (i=0;i<OUTHEIGHT;i++) /* loop over output rows */
   {     
     phi=3.14159265*(i/(OUTHEIGHT-1.0)-0.5);
     z=sin(phi);
     printf("Rendering line %d\n",i);

     for (t=0;t<OUTWIDTH;t++) /* loop over output columns */
       {
         g[t+OUTWIDTH*i]=128;  
         /* The default output color is gray=128 */
         theta=t*2*3.141592654/OUTWIDTH; 
         r=cos(phi);
         x=r*cos(theta);
         y=r*sin(theta);
         /* The spherical coordinate transform! */

         /* if we're inside the lines of the ball,
            color this sample black. */
         if (inlines(x,y,z,COLORWIDTH/2.0))
           g[t+OUTWIDTH*i]=0;

         else for (p=0;p<numspots;p++)
           {
          /* Brute force, loop through the list of
             points. We use the dot product of the
             pip location and the surface to determine
             if the sample is "within" the pip. If it
             is, we'll color the surface  with the
             pip color. */

             if (x*dotx[p]+y*doty[p]+z*dotz[p]>PIPCLOSE)
               g[t+OUTWIDTH*i]=155;

           }
       }
   }

   /* A simple output function, writes the image to disk. */
   writeimage("bbout.raw",g,OUTWIDTH,OUTHEIGHT);
}


int inlines(double x,double y, double z, double w)
{

/* This routine decides whether an XYZ point on the surface
   of the sphere is within a certain distance of one of
   the lines on the basketball. */

  double r;


  /* easy tests... Z=0, Y=0 equator and meridian lines */

  if (fabs(z)<w) return(1);
  if (fabs(y)<w) return(1);

  /* 45 degree stripes outside of curve. Near the equator
     another case needs to be used, though  */

  if (fabs(y)>35.0/BBRADIUS && fabs(z)>35.0/BBRADIUS)
    {
      if (fabs(z-y)<1.4142*w) return(1);
      if (fabs(y+z)<1.4142*w) return(1);
      return(0);
    }

  /* Do a trivial rejection of areas far from the curved
     part of the lines */

  if (fabs(x)<80.0/BBRADIUS) return(0); 

  /* curved bit. We assume the curve is in the shape of a 
     circle defined by some (empirically determined)
     constants. We find the radius from this circle and
     see whether we are within w units of it. We use y and z
     symmetrically depending on which half (x positive or 
     negative) of the ball we are on.  */

  if (x>0.0)
    {
      if (y>0.0) r=sqrt(z*z+(y-CIRCPOS/BBRADIUS)
                        *(y-CIRCPOS/BBRADIUS));
      else r=sqrt(z*z+(y+CIRCPOS/BBRADIUS)
                       *(y+CIRCPOS/BBRADIUS));
      if (fabs(r-CIRCRAD/BBRADIUS)<w) return(1); 
    }
  
  else
    {
      if (z>0.0) r=sqrt(y*y+(z-CIRCPOS/BBRADIUS)
                        *(z-CIRCPOS/BBRADIUS));
      else r=sqrt(y*y+(z+CIRCPOS/BBRADIUS)
                        *(z+CIRCPOS/BBRADIUS));
      if (fabs(r-CIRCRAD/BBRADIUS)<w) return(1); 
    }

  return(0);
}